1    import java.awt.geom.*;
2    import java.awt.*;
3    
4    public class QuadTree implements java.io.Serializable {
5    
6      private static final long serialVersionUID = Version.getSUID();
7    
8      Polygon boundary;
9    
10     int[] data = new int[100];
11     int nodes;
12     transient Line2D[] blines;
13   
14   //  public static Point2D emptypos = new Point2D.Double(0,0);
15     transient public Rectangle2D bounds = null;
16   
17     public int width;
18     public int height;
19     public double basex;
20     public double basey;
21     public double basescale;
22   
23     static double iscale = 1000;
24   //  static int treesamps = 2;
25     static double MAXERR = 0.2;
26     static double ERRORFLOOR = 0.01;
27     static double MINBOX = 1;
28   
29   //  static int fmask = 1<<31;
30   //  static int cmask = -1 & ~fmask;
31   
32     static int IN = 1<<16;
33     static int OUT = 2<<16;
34     static int CORNER = 4<<16;
35   
36     static int LMASK = (1<<16)-1;
37   
38   //  transient Rectangle pbounds = null;
39   
40     void pack() {
41       if (nodes*4<data.length) {
42         int[] newdata = new int[nodes*4];
43         for (int i=0;i<nodes*4;i++)
44           newdata[i] = data[i];
45         data = newdata;
46       }
47     }
48   
49     void expand(int newsize) {
50       if (newsize*4>data.length) {
51         int[] newdata = new int[newsize*6];
52         for (int i=0;i<nodes*4;i++)
53           newdata[i] = data[i];
54         data = newdata;
55       }
56     }
57   
58   //  private void readObject(java.io.ObjectInputStream stream)
59   //      throws java.io.IOException,java.lang.ClassNotFoundException {
60   
61   //    bounds = new Rectangle2D.Double(basex,basey,basescale,basescale);
62   //  }
63   
64     AffineTransform linearFit(int node) {
65   
66       int[] x = {data[4*data[node*4]],data[4*data[node*4+2]],
67           	data[4*data[node*4+1]],data[4*data[node*4+3]]};
68       int[] y = {data[4*data[node*4]+1],data[4*data[node*4+2]+1],
69           	data[4*data[node*4+1]+1],data[4*data[node*4+3]+1]};
70   
71       double x0 = (   x[0]+x[1]+x[2]+x[3] )/(4.0*iscale);
72       double y0 = (   y[0]+y[1]+y[2]+y[3] )/(4.0*iscale);
73   
74       double dxx = ( -x[0]+x[1]-x[2]+x[3] )/(2.0*iscale);
75       double dyy = ( -y[0]+y[2]-y[1]+y[3] )/(2.0*iscale);
76   
77       double dyx = ( -y[0]+y[1]-y[2]+y[3] )/(2.0*iscale);
78       double dxy = ( -x[0]+x[2]-x[1]+x[3] )/(2.0*iscale);
79   
80       return(new AffineTransform(dxx,dyx,dxy,dyy,x0,y0));
81     }
82   
83     AffineTransform getFit(Point2D[][] p) {
84   
85       double[] x = {p[0][0].getX(),p[1][0].getX(),
86                  p[0][1].getX(),p[1][1].getX()};
87       double[] y = {p[0][0].getY(),p[1][0].getY(),
88                  p[0][1].getY(),p[1][1].getY()};
89   
90       double x0 = (   x[0]+x[1]+x[2]+x[3] )/4.0;
91       double y0 = (   y[0]+y[1]+y[2]+y[3] )/4.0;
92   
93       double dxx = ( -x[0]+x[1]-x[2]+x[3] )/2.0;
94       double dyy = ( -y[0]+y[2]-y[1]+y[3] )/2.0;
95   
96       double dyx = ( -y[0]+y[1]-y[2]+y[3] )/2.0;
97       double dxy = ( -x[0]+x[2]-x[1]+x[3] )/2.0;
98   
99       return(new AffineTransform(dxx,dyx,dxy,dyy,x0,y0));
100    }
101  
102  
103    Point2D makeNode(int node, Point2D mpos, Mapping map, AffineTransform trans) {
104      Point2D.Double pos = new Point2D.Double();
105      map.ttrans(mpos,pos,trans);
106      data[node*4] = (int)(iscale*pos.getX());
107      data[node*4+1] = (int)(iscale*pos.getY());
108      data[node*4+2]=0;
109      return(pos);
110    }
111  
112    double getError(AffineTransform fit,double x,double y,Point2D tpos) {
113      Point2D mypos= new Point2D.Double(x,y);
114      fit.transform(mypos,mypos);
115      double dif = mypos.getX() - tpos.getX();
116      double err = dif*dif;
117      dif = mypos.getY() - tpos.getY();
118      return(err + dif*dif);
119    }
120  
121    void subDivide(int node, Rectangle2D box,
122      Mapping map, int[] lines, double error, AffineTransform trans, int code) {
123  
124  //    message.setText("width ="+box.getWidth()+" n = "+node+" of "+nodes+" "+str);
125  
126      int[][][] sublines = new int[2][2][];
127      Rectangle2D[][] subrect = new Rectangle2D[2][2];
128      Point2D[][] pos = new Point2D[2][2];
129      Point2D[][] mpos = new Point2D[2][2];
130      Point2D mytpos = new Point2D.Double(box.getX()+box.getWidth()/2,
131          			       box.getY()+box.getHeight()/2);
132      boolean dividechildren = false;
133      double myerror=0;
134  //    int mask = 0;
135  
136      map.ttrans(mytpos,mytpos,trans);
137  
138      expand(nodes+4);
139  
140      for (int i=0;i<=1;i++)
141        for (int j=0;j<=1;j++) {
142            data[4*node+2*i+j]=nodes++;
143            subrect[i][j] = new Rectangle2D.Double(box.getX()+i*box.getWidth()/2,
144          					 box.getY()+j*box.getHeight()/2,
145          					 box.getWidth()/2,box.getHeight()/2);
146            mpos[i][j] = new Point2D.Double(subrect[i][j].getX()+subrect[i][j].getWidth()/2,
147          				  subrect[i][j].getY()+subrect[i][j].getHeight()/2);
148            pos[i][j] = makeNode(data[4*node+2*i+j],mpos[i][j],map,trans);
149          }
150  
151      if (box.getWidth()>MINBOX || box.getHeight()>MINBOX) {
152        AffineTransform fit = linearFit(node);
153        myerror = getError(fit,0,0,mytpos);
154        for (int i=0;i<=1;i++)
155          for (int j=0;j<=1;j++)
156            myerror += getError(fit,i-0.5,j-0.5,pos[i][j]);
157        myerror = Math.sqrt(myerror);
158  //      System.out.println("code = "+Integer.toString(code,4)+" width ="+box.getWidth()+" n = "+node+ " e = "+myerror);//+" l = "+lines.length);
159  
160        if ((myerror > MAXERR || myerror > error) &&
161            myerror>ERRORFLOOR) 
162          dividechildren = true;
163  
164        if (Double.isInfinite(myerror) || Double.isNaN(myerror)) {
165          System.out.println("bad approximation "+myerror);
166          ErrorLog.log("bad approximation "+myerror);
167          dividechildren = false;
168        }
169  
170      }
171  
172      for (int i=0;i<=1;i++)
173        for (int j=0;j<=1;j++) {
174            Rectangle2D ibounds = map.btransBounds(subrect[i][j],trans,3);
175            int[] sub = new int[lines.length];
176            int count = 0;
177            for (int k=0;k<lines.length;k++)
178              if (blines[lines[k]].intersects(ibounds)) {
179                sub[count]=lines[k];
180  //              if ((count>0 && sub[0]!=lines[k-1]) &&
181  //        	  !(lines[k]==boundary.npoints-1 && sub[0]==0))
182  //        	dividechildren = true;
183                count++;
184              }
185            if (count>2 ||
186                (count==2 && sub[0]+1!=sub[1] &&
187                 (sub[0]!=0 || sub[1]!=boundary.npoints-1)))
188              dividechildren = true;
189  
190            sublines[i][j] = new int[count];
191            for (int k=0;k<count;k++) sublines[i][j][k] = sub[k];
192        }
193  
194  //    dividechildren = !(box.getWidth()<MINBOX && box.getHeight()<MINBOX);
195      if (box.getWidth()<MINBOX && box.getHeight()<MINBOX)
196        dividechildren = false;
197  
198      for (int i=0;i<=1;i++)
199        for (int j=0;j<=1;j++) {
200          if (dividechildren) {
201            subDivide(data[4*node+2*i+j], subrect[i][j],
202          	    map, sublines[i][j],myerror,trans,(code<<2)+2*i+j);
203          }
204          else if (sublines[i][j].length==0) {
205            if (boundary.contains(pos[i][j]))
206              data[data[4*node+2*i+j]*4+3]=IN;
207            else
208              data[data[4*node+2*i+j]*4+3]=OUT;
209          }
210          else if (sublines[i][j].length==1)
211            data[data[4*node+2*i+j]*4+3] = sublines[i][j][0];
212          else if (sublines[i][j][0]!=0 || sublines[i][j][1]!=boundary.npoints-1)
213            data[data[4*node+2*i+j]*4+3] = sublines[i][j][0] | CORNER;
214          else
215            data[data[4*node+2*i+j]*4+3] = sublines[i][j][1] | CORNER;
216  
217  //        if ((data[4*data[4*node+2*i+j]+2]&fmask)!=0) mask|=1;
218  //        else if ((data[4*data[4*node+2*i+j]+3]&fmask)!=0) mask|=2;
219  //        else mask |= 4;
220        }
221  //    if (mask==1) data[4*node+2] |= fmask;
222  //    if (mask==2) data[4*node+3] |= fmask;    
223    }
224  
225    public QuadTree(Mapping map,AffineTransform trans,
226          	  Rectangle2D bounds,Polygon boundary
227  //        	  ,javax.swing.JLabel message, String str
228      ) {
229      this.boundary = boundary;
230  
231      width = (int)bounds.getWidth()+1;
232      height = (int)bounds.getHeight()+1;
233      basescale = (bounds.getWidth()>bounds.getHeight() ? bounds.getWidth() : bounds.getHeight());
234      basex = bounds.getX();
235      basey = bounds.getY();
236  
237      this.bounds = new Rectangle2D.Double(basex,basey,basescale,basescale);
238  
239      int[] lines = new int[boundary.npoints]; 
240      
241      blines = new Line2D[boundary.npoints]; 
242  
243      for (int i=0;i<boundary.npoints;i++) {
244        blines[i] = new Line2D.Double(boundary.xpoints[i],boundary.ypoints[i],
245          			 boundary.xpoints[(i+1)%boundary.npoints],
246          			 boundary.ypoints[(i+1)%boundary.npoints]);
247        lines[i] = i;
248      }
249  
250      nodes = 1;
251      subDivide(0,this.bounds,
252                map,lines,0,trans,1);
253  
254      blines = null;
255      preCheck();
256  
257  //    System.out.println("quadtree has "+nodes+" nodes");
258  //    System.out.println("bounds = "+this.bounds+" size = ("+width+","+height+") base ="+basescale);
259  //    System.out.println("data.length = "+data.length);
260    }
261  
262    public static double delta = 1;
263  
264    double[] nx=null;
265    double[] ny=null;
266    double[] off=null;
267    boolean[] convex=null;
268  
269    void preCheck() {
270      nx = new double[boundary.npoints];
271      ny = new double[boundary.npoints];
272      off = new double[boundary.npoints];
273      convex = new boolean[boundary.npoints];
274      for (int p1=0;p1<boundary.npoints;p1++) {
275        int p2 = (p1+1)%boundary.npoints;
276        nx[p1] = boundary.ypoints[p2] - boundary.ypoints[p1];
277        ny[p1] = boundary.xpoints[p1] - boundary.xpoints[p2];
278  
279        double len = Math.sqrt(nx[p1]*nx[p1]+ny[p1]*ny[p1]);
280        if (len>0) {
281          nx[p1] /= len;
282          ny[p1] /= len;
283        }
284  
285        off[p1] = nx[p1]*boundary.xpoints[p1] + ny[p1]*boundary.ypoints[p1];
286  
287        if (boundary.contains((boundary.xpoints[p1]+boundary.xpoints[p2])*0.5+nx[p1]*delta,
288          		    (boundary.ypoints[p1]+boundary.ypoints[p2])*0.5+ny[p1]*delta)){
289          nx[p1]= -nx[p1];
290          ny[p1]= -ny[p1];
291          off[p1]= -off[p1];
292        }
293  
294        if (len>0 &&
295            boundary.contains(-delta/len*boundary.xpoints[p1]+(1+delta/len)*boundary.xpoints[p2],
296          		    -delta/len*boundary.ypoints[p1]+(1+delta/len)*boundary.ypoints[p2]))
297            convex[p1]=false;
298            else convex[p1]=true;
299  //      System.out.println("p1 = "+p1+"n = ("+nx[p1]+","+ny[p1]+") off = "+off[p1]+" convex = "+convex[p1]);
300      }
301    }
302  
303    boolean lcheck(double x,double y,int pos) {
304      if (pos==IN) return(true);
305      int p = pos&LMASK;
306  //    if ((pos&CORNER)!=0) {
307  //      System.out.println("("+x+","+y+") p = "+p+" in[p] "+(x*nx[p]+y*ny[p]<off[p])+
308  //        		 " in[p+] "+(x*nx[(p+1)%boundary.npoints]+y*ny[(p+1)%boundary.npoints]<off[(p+1)%boundary.npoints])+
309  //        		 " convex = "+convex[p]);
310  //    }
311      if (x*nx[p]+y*ny[p]<off[p]) {
312        if ((pos&CORNER)==0 || !convex[p]) return(true);
313      }
314      else
315        if ((pos&CORNER)==0 || convex[p]) return(false);
316  
317      p  = (p+1)%boundary.npoints;
318      return(x*nx[p]+y*ny[p]<off[p]);
319    }
320  
321    Point2D fillCoords(int[][][] coords,int node,Rectangle2D box,int code) {
322  //    System.out.println("code "+Integer.toString(code,4)+" size = "+box.getWidth()+" n = "+node+
323  //        	       " at ("+(int)box.getX()+","+(int)box.getY()+")");
324      int expanded = 0;
325      Point2D[][] children = new Point2D[2][2];
326      Point2D pos = new Point2D.Double();
327  
328      for (int i=0;i<=1;i++)
329        for (int j=0;j<=1;j++) {
330          int child = data[node*4+2*i+j];
331          if ((data[4*child+2])==0)
332            children[i][j] = new Point2D.Double(data[child*4]/iscale,data[child*4+1]/iscale);
333          else {
334            children[i][j] = fillCoords(coords,data[4*node+2*i+j],
335          			      new Rectangle2D.Double(box.getX()+i*box.getWidth()/2,
336          						     box.getY()+j*box.getHeight()/2,
337          						     box.getWidth()/2,box.getHeight()/2),
338          			      (code<<2)+2*i+j);
339            expanded++;
340          }
341        }
342  
343      if (expanded<4) {
344  
345        int imin = (int)box.getX();
346        double imax = box.getWidth();
347        int jmin = (int)box.getY();
348        double jmax = box.getHeight();
349  
350        AffineTransform fit = getFit(children);
351  
352        for (int n=0;n<=1;n++)
353          for (int m=0;m<=1;m++) {
354            int child = data[node*4+ n*2+ m];
355            if (data[4*child+2]==0) {
356  
357              if (data[4*child+3]==OUT)
358                for (int i=(int)(n*imax/2.0);i<imax/(2.0-n) && i<width-imin;i++)
359          	for (int j=(int)(m*jmax/2.0);j<jmax/(2.0-m) && j<height-jmin;j++)
360          	  coords[i+imin][j+jmin][0] = -1; //&~(1<<20);
361              
362              else
363                for (int i=(int)(n*imax/2.0);i<imax/(2.0-n) && i<width-imin;i++)
364          	for (int j=(int)(m*jmax/2.0);j<jmax/(2.0-m) && j<height-jmin;j++) {
365  
366          	  pos.setLocation(2*(i+imin-box.getX())/box.getWidth()-1,
367          			  2*(j+jmin-box.getY())/box.getHeight()-1);
368          	  fit.transform(pos,pos);
369  
370          	  if (lcheck(pos.getX(),pos.getY(),data[4*child+3])) {
371          	    coords[i+imin][j+jmin][0] = (int)(pos.getX()+0.5);
372          	    coords[i+imin][j+jmin][1] = (int)(pos.getY()+0.5);
373  //        	    if (data[4*child+3]==IN) coords[i+imin][j+jmin][0] |= (2<<20);
374          	  }
375          	  else
376          	    coords[i+imin][j+jmin][0] = -1; //&~(2<<20);
377          	}
378            }
379          }
380      }
381  
382      return(new Point2D.Double((children[0][0].getX()+children[0][0].getX()+
383          		       children[0][0].getX()+children[0][0].getX())/4,
384          		      (children[0][0].getY()+children[0][0].getY()+
385          		       children[0][0].getY()+children[0][0].getY())/4));
386      }
387  
388    public int[][][] getCoords() {
389  //    pbounds = boundary.getBounds();
390  //    System.out.println("getCoords started");
391      bounds = new Rectangle2D.Double(basex,basey,basescale,basescale);
392      preCheck();
393      int[][][] coords = new int[width][height][2];
394   
395      fillCoords(coords,0,new Rectangle2D.Double(0,0,basescale,basescale),1);
396  
397      return(coords);
398    }
399  
400  }
401