Package com.nr.cg

Source Code of com.nr.cg.Delaunay$Hash1

package com.nr.cg;

import static com.nr.cg.Triel.*;
import com.nr.ran.Hash;
import com.nr.ran.Ranhash;

import static com.nr.NRUtil.*;

/**
* Delaunay triangulation
* Copyright (C) Numerical Recipes Software 1986-2007
* Java translation Copyright (C) Huang Wen Hui 2012
*
* @author hwh
*
*/
public class Delaunay {
  static final int DIM=2;
  private static final double fuzz= 1.0e-6, bigscale=1000.0;
 
  /**
   * Random number counter
   */
  private static int jran = 14921620;
  public int npts,ntri,ntree,ntreemax,opt;
  public double delx,dely;
  public Point[] pts;
  public Triel[] thelist;
 
  /**
   * Create the hash memories with null hash function.
   */
  public Hash1 linehash;
  public Hash1 trihash;
  public int[] perm;
 
  /**
   * The raw hash function
   */
  public Ranhash hashfn = new Ranhash();
 
  public class Hash1 extends Hash<Long,Integer>{
    public Hash1(final int nh, final int nm){
      super(nh, nm);
    }
   
    public long fn(final Long  k) {
      return k & 0x7FFFFFFFFFFFFFFFL;
    }
  }

  public int dim(){return DIM;}

  public Delaunay(final Point[] pvec){
    this(pvec,0);
  }
 
  /**
   * Construct Delaunay triangulation from a vector of points pvec. If bit 0 in options is nonzero,
   * hash memories used in the construction are deleted. (Some applications may want to use them
   * and will set options to 1.)
   *
   * @param pvec
   * @param options
   */
  public Delaunay(final Point[] pvec, final int options) {
    for(int i=0;i<pvec.length;i++)
      if(pvec[i].dim() != DIM)
        throw new IllegalArgumentException("Need same dim!");
   
    npts = pvec.length;
    ntri=0;
    ntree=0;
    ntreemax = 10*npts+1000;
    opt=options;
    pts = new Point[npts+3];
    thelist = new Triel[ntreemax];
    for(int i=0;i<ntreemax;i++)
      thelist[i] = new Triel();
   
    int j;
    double xl,xh,yl,yh;
    linehash = new Hash1(6*npts+12,6*npts+12);
    trihash = new Hash1(2*npts+6,2*npts+6);
    perm = new int[npts];
    xl = xh = pvec[0].x[0];
    yl = yh = pvec[0].x[1];
    for (j=0; j<npts; j++) {
      pts[j] = pvec[j];
      perm[j] = j;
      if (pvec[j].x[0] < xl) xl = pvec[j].x[0];
      if (pvec[j].x[0] > xh) xh = pvec[j].x[0];
      if (pvec[j].x[1] < yl) yl = pvec[j].x[1];
      if (pvec[j].x[1] > yh) yh = pvec[j].x[1];
    }
    delx = xh - xl;
    dely = yh - yl;
    pts[npts] = new Point(0.5*(xl + xh), yh + bigscale*dely);
    pts[npts+1] = new Point(xl - 0.5*bigscale*delx,yl - 0.5*bigscale*dely);
    pts[npts+2] = new Point(xh + 0.5*bigscale*delx,yl - 0.5*bigscale*dely);
    storetriangle(npts,npts+1,npts+2);
    for (j=npts; j>0; j--) {
      //SWAP(perm[j-1],perm[hashfn.int64p(jran++) % j]);
      int nnn = (int)(hashfn.int64p(jran++) % j);
      swap(perm, j-1, nnn);
    }
    for (j=0; j<npts; j++) insertapoint(perm[j]);
    for (j=0; j<ntree; j++) {
      if (thelist[j].stat > 0) {
        if (thelist[j].p[0] >= npts || thelist[j].p[1] >= npts ||
        thelist[j].p[2] >= npts) {
          thelist[j].stat = -1;
          ntri--;
        }
      }
    }
    /*
    if (!(opt & 1)) {
      delete [] perm;
      delete trihash;
      delete linehash;
    }
    */
  }
 
  /**
   * Add the point with index r incrementally to the Delaunay triangulation
   *
   * @param r
   */
  public void insertapoint(final int r) {
    int i,j,k,l,s,tno=0,ntask,d0,d1,d2;
    long key;
    int[] tasks = new int[50], taski = new int[50], taskj = new int[50];
    for (j=0; j<3; j++) {
      tno = whichcontainspt(pts[r],1);
      if (tno >= 0) break;
      pts[r].x[0] += fuzz * delx * (hashfn.doub(jran++)-0.5);
      pts[r].x[1] += fuzz * dely * (hashfn.doub(jran++)-0.5);
    }
    if (j == 3) throw new IllegalArgumentException("points degenerate even after fuzzing");
    ntask = 0;
    i = thelist[tno].p[0]; j = thelist[tno].p[1]; k = thelist[tno].p[2];
    if ((opt & 2)!=0 && i < npts && j < npts && k < npts) return;
    d0 =storetriangle(r,i,j);
    tasks[++ntask] = r; taski[ntask] = i; taskj[ntask] = j;
    d1 = storetriangle(r,j,k);
    tasks[++ntask] = r; taski[ntask] = j; taskj[ntask] = k;
    d2 = storetriangle(r,k,i);
    tasks[++ntask] = r; taski[ntask] = k; taskj[ntask] = i;
    erasetriangle(i,j,k,d0,d1,d2);
    while (ntask!=0) {
      s=tasks[ntask]; i=taski[ntask]; j=taskj[ntask--];
      key = hashfn.int64(j) - hashfn.int64(i);
      Integer[] l_w = new Integer[1];
      if (linehash.get(key,l_w, 0)==0 ) continue;
      l = l_w[0];
      if (incircle(pts[l],pts[j],pts[s],pts[i]) > 0.0){
        d0 = storetriangle(s,l,j);
        d1 = storetriangle(s,i,l);
        erasetriangle(s,i,j,d0,d1,-1);
        erasetriangle(l,j,i,d0,d1,-1);
        key = hashfn.int64(i)-hashfn.int64(j);
        linehash.erase(key);
        key = 0 - key;
        linehash.erase(key);
        tasks[++ntask] = s; taski[ntask] = l; taskj[ntask] = j;
        tasks[++ntask] = s; taski[ntask] = i; taskj[ntask] = l;
      }
    }
  }
 
  public int whichcontainspt(final Point p){
    return whichcontainspt(p,0);
  }
 
  /**
   * Given point p, return index in thelist of the triangle in the
   * triangulation that contains it, or return 1 for failure. If
   * strict is nonzero, require strict containment, otherwise allow
   * the point to lie on an edge.
   *
   * @param p
   * @param strict
   * @return
   */
  public int whichcontainspt(final Point p, final int strict) {
    int i,j=0,k=0;
    while (thelist[k].stat <= 0) {
      for (i=0; i<3; i++) {
        if ((j = thelist[k].d[i]) < 0) continue;
        if (strict!=0) {
          if (thelist[j].contains(p) > 0) break;
        } else {
          if (thelist[j].contains(p) >= 0) break;
        }
      }
      if (i == 3) return -1;
      k = j;
    }
    return k;
  }

  /**
   * Erase triangle abc in trihash and inactivate it in thelist after
   * setting its daughters.
   *
   * @param a
   * @param b
   * @param c
   * @param d0
   * @param d1
   * @param d2
   */
  public void erasetriangle(final int a, final int b, final int c, final int d0, final int d1, final int d2) {
    long key;
    int j;
    key = hashfn.int64(a) ^ hashfn.int64(b) ^ hashfn.int64(c);
    Integer[] j_w = new Integer[1];
    if (trihash.get(key,j_w,0) == 0) throw new IllegalArgumentException("nonexistent triangle");
    j = j_w[0];
    trihash.erase(key);
    thelist[j].d[0] = d0; thelist[j].d[1] = d1; thelist[j].d[2] = d2;
    thelist[j].stat = 0;
    ntri--;
  }

  /**
   * Store a triangle with vertices a, b, c in trihash. Store its points in
   * linehash under keys to opposite sides. Add it to thelist, returning its
   * index there.
   *
   * @param a
   * @param b
   * @param c
   * @return
   */
  public int storetriangle(final int a, final int b, final int c) {
    long key;
    thelist[ntree].setme(a,b,c,pts);
    key = hashfn.int64(a) ^ hashfn.int64(b) ^ hashfn.int64(c);
    trihash.set(key,ntree);
    key = hashfn.int64(b)-hashfn.int64(c);
    linehash.set(key,a);
    key = hashfn.int64(c)-hashfn.int64(a);
    linehash.set(key,b);
    key = hashfn.int64(a)-hashfn.int64(b);
    linehash.set(key,c);
    if (++ntree == ntreemax) throw new IllegalArgumentException("thelist is sized too small");
    ntri++;
    return (ntree-1);
  }
 
  public double interpolate(final Point p, final double[] fnvals, final double defaultval) {
    int n,i,j,k;
    double[] wgts = new double[3];
    int[] ipts = new int[3];
    double sum, ans = 0.0;
    n = whichcontainspt(p);
    if (n < 0) return defaultval;
    for (i=0; i<3; i++) ipts[i] = thelist[n].p[i];
    for (i=0,j=1,k=2; i<3; i++,j++,k++) {
      if (j == 3) j=0;
      if (k == 3) k=0;
      wgts[k]=(pts[ipts[j]].x[0]-pts[ipts[i]].x[0])*(p.x[1]-pts[ipts[i]].x[1])
        - (pts[ipts[j]].x[1]-pts[ipts[i]].x[1])*(p.x[0]-pts[ipts[i]].x[0]);
    }
    sum = wgts[0] + wgts[1] + wgts[2];
    if (sum == 0) throw new IllegalArgumentException("degenerate triangle");
    for (i=0; i<3; i++) ans += wgts[i]*fnvals[ipts[i]]/sum;
    return ans;
  }
}
TOP

Related Classes of com.nr.cg.Delaunay$Hash1

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.