1    /*
2    Copyright 2000 by Ralph Hartley
3    This software is licenced under the terms of the
4    Gnu Public Licence
5    */
6    import java.awt.*;
7    import java.util.*;
8    import java.awt.geom.*;
9    import java.awt.image.*;
10   import javax.swing.*;
11   
12   public class Mapping implements Dragable,java.io.Serializable{
13   
14     private static final long serialVersionUID = Version.getSUID();
15   
16     public int stacnt  = 0;
17     public Vertex[] spoint;
18     float[][] im;
19     int stacntplus = 0;
20     public boolean degenerate = false;
21   
22     static public double maxslope = 20;
23   
24     public boolean isVertical = false;
25   
26     double[] xcoif,ycoif;
27     boolean solved = false;
28     float[][] mapx;
29     transient int selected = -1;
30     transient int[] dragpos;
31     transient int oldi;
32     transient int oldj;
33     transient java.awt.Graphics draggraphics;
34     transient double dscale=1;
35   
36     static public double beta;
37     static public double theta;
38     static public double sideoff;
39     static public double mindist;
40     static public double linethresh;
41   
42     HashSet nearverticals = new HashSet();
43   
44     int wide1 = 0;
45     int wide2 = 0;
46     AffineTransform overlaytrans = null;
47     boolean tchanged = false;
48   
49     double maxdist = 0;
50   
51     int width = -1;
52     int height = -1;
53   
54     static double MAXERR = 0.0001;
55     static int MAXINTS = 15;
56   
57     transient QuadTree tree = null;
58   
59   //  transient SoftReference scoords = new SoftReference(null);
60   //  transient SoftReference spixels = new SoftReference(null);
61   
62     transient public SegEditor editor = null;
63   
64     public Mapping(int snum,boolean isVertical){
65       this.isVertical = isVertical;
66       spoint = new Vertex[snum+1];
67       im = new float[snum+1][2];
68     }
69   
70     public boolean used(Vertex station) {
71       for (int i=0;i<stacnt;i++)
72         if (spoint[i]==station) return(true);
73       return(false);
74     }
75   
76     public void remove(Vertex station) {
77       for (int i=0;i<stacnt;i++)
78         if (spoint[i]==station)
79           unmark(i);
80     }
81   
82     public void setScale(double scale) {
83       setOverlayTrans();
84       if (overlaytrans==null||isVertical)
85         overlaytrans = AffineTransform.getScaleInstance(scale,-scale);
86     }
87   
88     public double getScale() {
89       if (overlaytrans==null) return(Segment.defscale);
90       return(overlaytrans.getScaleX());
91     }
92   
93     void adjustCount() {
94       if (stacnt>(spoint.length-10)) {
95         solved = false;
96         Vertex[] nsp = new Vertex[spoint.length * 2];
97         float[][] nim = new float[spoint.length * 2][];
98         for (int i=0;i<stacnt;i++) {
99           nsp[i] = spoint[i];
100          nim[i] = im[i];
101        }
102        for (int i=stacnt; i<spoint.length*2; i++)
103          nim[i] = new float[2];
104        spoint = nsp;
105        im = nim;
106      }
107      else if (stacnt<(spoint.length/3) && spoint.length>30) {
108        solved=false;
109        Vertex[] nsp = new Vertex[spoint.length/2];
110        float[][] nim = new float[spoint.length/2][];
111        for (int i=0;i<spoint.length/2;i++) {
112          nsp[i] = spoint[i];
113          nim[i] = im[i];
114        }
115        spoint = nsp;
116        im = nim;
117      }
118    }
119  
120    public void reattachStations(Survey survey) {
121      for (int i=0;i<stacnt;i++) {
122        spoint[i] = survey.getStation(spoint[i]);
123      }
124      solved = false;
125    }
126  
127    public void btrans(Point2D p,Point2D res) {
128      if (isVertical) {
129        overlaytrans.transform(p,res);
130        return;
131      }
132      double x = (float)p.getX();
133      double y = (float)p.getY();
134  
135      int i;
136      double u = xcoif[stacntplus] + xcoif[stacntplus+1] * x + xcoif[stacntplus+2] *y;
137      double v = ycoif[stacntplus] + ycoif[stacntplus+1] * x + ycoif[stacntplus+2] *y;
138  
139      for (i=0;i<stacntplus;i++) {
140        double localh = h(x-mapx[i][0],y-mapx[i][1]);
141        u += xcoif[i] * localh;
142        v += ycoif[i] * localh;
143      }
144  
145      res.setLocation(u,v);
146    }
147  
148  /*                 
149    public double[] getMorphData(Point2D p) {
150  
151      double x = (float)p.getX();
152      double y = (float)p.getY();
153  
154      double u = xcoif[stacntplus] + xcoif[stacntplus+1] * x + xcoif[stacntplus+2] *y;
155      double v = ycoif[stacntplus] + ycoif[stacntplus+1] * x + ycoif[stacntplus+2] *y;
156  
157      for (int i=0;i<stacntplus;i++) {
158        double localh = h(x-mapx[i][0],y-mapx[i][1]);
159        u += xcoif[i] * localh;
160        v += ycoif[i] * localh;
161      }
162  
163      double[][] grad = new double[2][2];
164  
165      getGrad(grad,u,v);
166      System.out.println("grad = "+grad);
167  
168      double[] res = new double[4];
169  
170      res[0] = u;
171      res[1] = v;
172      res[2] = Math.atan2(grad[0][0]+grad[1][1],grad[0][1]-grad[1][0]);
173      res[3] = Math.sqrt(grad[0][1]*grad[1][0] - grad[1][1]*grad[0][0]);
174  
175      System.out.println(" det = "+(grad[0][0]*grad[1][1]-grad[0][1]*grad[1][0]));
176  
177      System.out.println("res x = "+res[0]+" y = "+res[1]+
178        " theta = "+ res[2]*180/Math.PI + " scale = "+ res[3]);
179  
180      return(res);
181    }
182  */               
183  /**
184   * Finish the drag perminantly fixing the selected item at the
185   * location of the last move.
186   * Return true if succesfull, otherwize item is still selected.
187   * @return boolean
188   */
189    public void fix() {
190      if (draggraphics==null) return;
191      im[selected][0] = oldi;
192      im[selected][1] = oldj;
193  //    if (editor!=null) editor.select(spoint[selected]);    (it is already selected)
194      setOverlayTrans();
195  //    System.out.println("trans = "+overlaytrans);
196    }
197  
198    public void delete() {
199      if (draggraphics==null) return;
200      if (stacnt!=0) {
201        unmark(selected);
202        setOverlayTrans();
203      }
204    }
205  
206    public void unmark(int p) {
207      solved=false;
208      for (int i=p;i<stacnt-1;i++) {
209        spoint[i] = spoint[i+1];
210        im[i][0] = im[i+1][0];
211        im[i][1] = im[i+1][1];
212      }
213      stacnt--;
214      adjustCount();
215    }
216  
217    public void setSelection(Vertex v) {
218      selected = -1;
219      if (v==null) return;
220      for (int i=0;i<stacnt;i++)
221        if (v==spoint[i]) selected = i;
222  //    if (editor!=null) editor.select(spoint[selected]);
223    }
224  
225    public void locate(Object what,int x,int y){
226      if (isVertical) {
227        ErrorLog.log("This segment is marked as 'vertical',\n"+
228          	   "meaning that it is a cross section.\n"+
229          	   "Point based morphing is not implemented\nfor cross sections.");
230        return;
231      }
232      Vertex v = null;
233      if (what instanceof Vertex) v = (Vertex)what;
234      else return;
235      for (int i=0;i<stacnt;i++) {
236        if (v==spoint[i]) return;
237  //      if (checkVertical(v,spoint)) return;;
238      }
239      solved = false;
240      spoint[stacnt] = v;
241      im[stacnt][0] = x;
242      im[stacnt][1] = y;
243      stacnt++;
244  
245      adjustCount();
246  
247      setOverlayTrans();
248  
249      selected = stacnt-1;
250  
251  //    if (editor!=null) editor.select(spoint[selected]);
252    }
253  
254    public AffineTransform setOverlayTrans() {
255      if (isVertical) {
256        return(overlaytrans);
257      }
258  
259      int i,j;
260      maxdist=0;
261      double dist=0;
262  
263      for (i=0;i<stacnt;i++)
264        for (j=i+1;j<stacnt;j++) {
265          dist = spoint[i].position.hdist(spoint[j].position);
266          double vdist = spoint[i].position.vdist(spoint[j].position);
267          if (dist==0) unmark(j--);
268          else if (vdist>dist*maxslope &&
269            !nearverticals.contains(new Line(spoint[i],spoint[j]))) {
270            String options[] = {"unmark "+spoint[i].label,"unmark "+spoint[j].label,
271          		    "keep both","unmark both"};
272            int res = JOptionPane.showOptionDialog(
273              null,"Near vertical shots may give bad morphing results",
274              "Vertical",0,JOptionPane.QUESTION_MESSAGE,
275              null,options,null);
276            if (res==0) {
277              unmark(i--);
278              break;
279            }
280            if (res==1) unmark(j--);
281            if (res==2)
282              nearverticals.add(new Line(spoint[i],spoint[j]));
283            if (res==3) {
284              unmark(j);
285              unmark(i--);
286              break;
287            }
288          }
289          if (dist>maxdist) {
290            maxdist = dist;
291            wide1=i;
292            wide2=j;
293          }
294        }
295  //    System.out.println("wide1 = spoint["+wide1+"] ="+spoint[wide1].label+" at "+spoint[wide1].position + "  " +
296  //        	       "wide1 = spoint["+wide2+"] ="+spoint[wide2].label+" at "+spoint[wide2].position);
297  //    System.out.println("dist = "+maxdist);
298  //    System.out.println("image1 = ("+im[wide1][0]+","+im[wide1][1]+") "+
299  //        	       "image2 = ("+im[wide2][0]+","+im[wide2][1]+")");
300  
301      if (maxdist==0) return(overlaytrans = null);
302  
303      double a = ((spoint[wide2].position.X[0]-spoint[wide1].position.X[0])*(im[wide2][0]-im[wide1][0]) -
304          	(spoint[wide2].position.X[1]-spoint[wide1].position.X[1])*(im[wide2][1]-im[wide1][1]))/
305        (maxdist*maxdist);
306      double b = ((spoint[wide2].position.X[1]-spoint[wide1].position.X[1])*(im[wide2][0]-im[wide1][0]) +
307          	(spoint[wide2].position.X[0]-spoint[wide1].position.X[0])*(im[wide2][1]-im[wide1][1]))/
308        (maxdist*maxdist);
309  
310  //    System.out.println("a = "+a+" b = "+b);
311  
312      AffineTransform newoverlaytrans = 
313        new AffineTransform(a,b,b,-a,
314          		  im[wide1][0]-a*spoint[wide1].position.X[0]-b*spoint[wide1].position.X[1],
315          		  im[wide1][1]-b*spoint[wide1].position.X[0]+a*spoint[wide1].position.X[1]);
316  
317      if (overlaytrans==null || !overlaytrans.equals(newoverlaytrans)) {
318        overlaytrans = newoverlaytrans;
319        tchanged = true;
320      }
321      return(overlaytrans);
322    }
323  
324    public void getGrad(double[][] A, double x, double y) {
325  
326      A[0][0] = xcoif[stacntplus+1];
327      A[0][1] = xcoif[stacntplus+2];
328      A[1][0] = ycoif[stacntplus+1];
329      A[1][1] = ycoif[stacntplus+2];
330  
331  
332      for (int i=0;i<stacntplus;i++) {
333        double[] grad = gradh(x-mapx[i][0],y-mapx[i][1]);
334        A[0][0] += xcoif[i] * grad[0];
335        A[0][1] += xcoif[i] * grad[1];
336        A[1][0] += ycoif[i] * grad[0];
337        A[1][1] += ycoif[i] * grad[1];
338      }
339  
340    }
341    
342  // ftrans is the tricky one. The mapping we use does not give
343  // this directly. So we do it by approximation (Newton).
344    public Point2D ftrans(Point2D p,Point2D.Double guess) {
345      if (isVertical)
346        try {
347          return(overlaytrans.inverseTransform(p,guess));
348        } catch (NoninvertibleTransformException ex) {ErrorLog.exception(ex,"bad transform (internal error)");}
349  
350      // as a starting guess use the location of the closest known point
351      double min = (p.getX()-im[0][0])*(p.getX()-im[0][0])+(p.getY()-im[0][1])*(p.getY()-im[0][1]);
352      int mini = 0;
353      double dist;
354      int i;
355      for (i=1;i<stacntplus;i++) {
356        dist = (p.getX()-im[i][0])*(p.getX()-im[i][0])+(p.getY()-im[i][1])*(p.getY()-im[i][1]);
357        if (dist<min) {
358          min=dist;
359          mini = i;
360        }        	
361      }        		
362      guess.setLocation(mapx[mini][0],mapx[mini][1]);
363      Linset jacobian = new Linset(2);
364      Point2D q = new Point2D.Double();
365  //    float[] e = new float[2];
366      double[] grad = new double[2];
367      float error;
368  //System.out.println("ftrans");
369      for (int j=0;j<MAXINTS;j++) {
370        btrans(guess,q);
371        jacobian.b[0] = p.getX()-q.getX();
372        jacobian.b[1] = p.getY()-q.getY();
373        error = (float)(jacobian.b[0]*jacobian.b[0]+jacobian.b[1]*jacobian.b[1]);
374  //        System.out.println("error = "+error);
375        if (error<MAXERR) break;
376        getGrad(jacobian.A,guess.getX(),guess.getY());
377  
378        jacobian.ludcmp();
379        jacobian.lubksb();
380        guess.x += jacobian.b[0];
381        guess.y += jacobian.b[1];
382   
383  //      if (j==MAXINTS-1) System.out.println("Max ints");
384     }
385      return guess;
386    }
387  
388    Rectangle2D btransBounds(Rectangle2D box,AffineTransform trans, double margin) {
389      Point2D pos = new Point2D.Double(box.getX(),box.getY());
390      ttrans(pos,pos,trans);
391      Rectangle2D res = new Rectangle2D.Double(pos.getX(),pos.getY(),0,0);
392      pos.setLocation(box.getX(),box.getMaxY());
393      ttrans(pos,pos,trans);
394      res.add(pos);
395      pos.setLocation(box.getMaxX(),box.getY());
396      ttrans(pos,pos,trans);
397      res.add(pos);
398      pos.setLocation(box.getMaxX(),box.getMaxY());
399      ttrans(pos,pos,trans);
400      res.add(pos);
401      res.add(res.getX()-margin,res.getY()-margin);
402      res.add(res.getMaxX()+margin,res.getMaxY()+margin);
403      return(res);
404    }
405  
406    public void ttrans(Point2D pos,Point2D res,AffineTransform trans) {
407      try {
408        trans.inverseTransform(pos,res);
409      } catch (java.awt.geom.NoninvertibleTransformException ex) {
410        ErrorLog.exception(ex);
411      }
412      btrans(res,res);
413    }
414  
415    public BufferedImage getImage(Image raw,AffineTransform trans,
416          			Rectangle2D bounds,Polygon boundary,
417          			String str,
418          			boolean transparent) throws CookKillException {
419  
420      if (raw==null) return(null);
421  
422      adjustCount();
423  
424  //    solve();
425  //    message.setText("Building quadtree "+str);
426      CartoFrame.message.setText("Reading "+str);
427  
428      int pixels[] = null;
429  //    int width=-1,height=-1;
430  
431      try {
432        for (/*width= -1*/;width<=0;width = raw.getWidth(null)) Thread.sleep(50);
433        for (/*height= -1*/;height<=0;height = raw.getHeight(null)) Thread.sleep(50);
434      } catch(InterruptedException e){ErrorLog.exception(e);};
435      pixels = new int[width*height];
436      PixelGrabber pg = new PixelGrabber(raw,0,0,width,height,pixels,0,width);
437      try {
438        pg.grabPixels();
439      } catch(InterruptedException e){ErrorLog.exception(e);};
440      pixels[0] = 0;
441  
442      CartoFrame.message.setText("Morphing "+str);
443  
444      int w = (int)bounds.getWidth() + 1;
445      int h = (int)bounds.getHeight() + 1;
446  
447      int[] tpixels = new int[w*h];
448      int x,y;
449      int progress = -1;
450  
451      if (tree==null) tree = new QuadTree(this,trans,bounds,boundary);
452      int[][][] coords = tree.getCoords();
453      tree = null;
454      int[] pos=null;
455  
456      for (x=0;x<w;x++) {
457        if (x*20/w > progress) {
458          if (Carto.kill==Carto.COOKKILL) {
459            throw new CookKillException();
460          }
461          progress = x*20/w;
462          CartoFrame.message.setText(str + " %" + (5*progress));
463        }
464        for (y=0;y<h;y++) {
465          int p = x+w*y;
466          pos = coords[x][y];
467          if (pos[0]>=0 && pos[0]<width-1 &&   //check for karma - a bad boundary could fail.
468              pos[1]>=0 && pos[1]<height-1) {
469            tpixels[p] = pixels[pos[0]+pos[1]*width];
470            if (transparent &&
471                (0==((~tpixels[p])&((1<<24)-1))))
472              tpixels[p]=0;
473          }
474          else tpixels[p]=0;
475        }
476      }
477  
478  //System.out.println("mapping done by getimage");
479  //    message.setText("Finishing "+str);
480      pixels = null;
481      coords = null;
482      BufferedImage res = new BufferedImage(w,h,BufferedImage.TYPE_INT_ARGB);
483      res.getRaster().setDataElements(0,0,w,h,tpixels);
484      tpixels = null;
485      return(res);
486    }
487  
488    double[] gradh(double x, double y) {
489      double[] res = new double[2];
490      if (x==0 && y==0) { res[0]=res[1]=0.0; return res;}
491      double mag = theta*(1+Math.log(x*x+y*y));
492      res[0]= x*mag;
493      res[1]= y*mag;
494      return res;
495    }
496  
497    double h(double x, double y) {
498      return (theta*(x*x+y*y)*Math.log(x*x+y*y+mindist));
499    }
500  
501    double hsq(double xsq) {
502      return (theta*xsq*Math.log(xsq+mindist));
503    }
504  
505    double hrange(float[] mapx, Rectangle2D zone) {
506  
507      int outcode = zone.outcode(mapx[0],mapx[1]);
508      double mindistsq = 0;
509      if (0!=(outcode & Rectangle2D.OUT_LEFT))
510        mindistsq += (mapx[0]-zone.getX())*(mapx[0]-zone.getX());
511      else if (0!=(outcode & Rectangle2D.OUT_RIGHT))
512        mindistsq += (mapx[0]-zone.getMaxX())*(mapx[0]-zone.getMaxX());
513      if (0!=(outcode & Rectangle2D.OUT_TOP))
514        mindistsq += (mapx[1]-zone.getY())*(mapx[1]-zone.getY());
515      else if (0!=(outcode & Rectangle2D.OUT_BOTTOM))
516        mindistsq += (mapx[1]-zone.getMaxY())*(mapx[1]-zone.getMaxY());
517  
518      double maxdistsq;
519      if (mapx[0]>zone.getX()+zone.getWidth()/2.0)
520        maxdistsq = (mapx[0] - zone.getX()) * (mapx[0] - zone.getX());
521      else 
522        maxdistsq = (mapx[0] - zone.getMaxX()) * (mapx[0] - zone.getMaxX());
523  
524      if (mapx[1]>zone.getY()+zone.getHeight()/2.0)
525        maxdistsq += (mapx[1] - zone.getY()) * (mapx[1] - zone.getY());
526      else 
527        maxdistsq += (mapx[1] - zone.getMaxY()) * (mapx[1] - zone.getMaxY());
528  
529      double minh = hsq(mindistsq);
530      double maxh = hsq(maxdistsq);
531  
532      System.out.println("mind = "+Math.sqrt(mindistsq)+" h ="+minh+"maxd = "+Math.sqrt(maxdistsq)+" h ="+maxh);
533  
534      if (mindistsq<1/Math.E && maxdistsq<1/Math.E) return(minh-maxh);
535      if (mindistsq<1/Math.E && maxdistsq>1/Math.E) {
536        if (minh<maxh) return(maxh-theta*2.0/Math.E);
537        else return(minh-theta*2.0/Math.E);
538      }
539      return(maxh-minh);
540    }
541  
542  /**
543   * provisionally move the selected item.
544   * Return true if the new position is OK
545   * @param x float
546   * @param y float
547   */
548  public void move(int x,int y)
549  {
550    if (draggraphics==null) return;
551    showStation(draggraphics,oldi/dscale,oldj/dscale);
552          
553    oldi = (int)im[selected][0]+x-dragpos[0];
554    oldj = (int)im[selected][1]+y-dragpos[1];
555  
556    showStation(draggraphics,oldi/dscale,oldj/dscale);
557  
558    solved=false;
559  
560    if (editor!=null) editor.select(spoint[selected]);
561  }
562          
563  /**
564   * select the item closest to the position given. Return true if an
565   * item was selected. g is a graphics object on which to display the process.
566   * @return boolean
567   * @param e java.awt.Event
568   */
569  public boolean select(int x,int y, java.awt.Graphics g)
570  {
571          float dist;
572          
573          dragpos=new int[2];
574          dragpos[0] = x;
575          dragpos[1] = y;
576  
577          if (stacnt==0) return false;
578  
579  //        int oldselect = selected;
580          selected = 0;
581          float min = (dragpos[0]-im[0][0])*(dragpos[0]-im[0][0])+(dragpos[1]-im[0][1])*(dragpos[1]-im[0][1]);
582  
583          for (int i=1;i<stacnt;i++) {
584          	dist = (dragpos[0]-im[i][0])*(dragpos[0]-im[i][0])+(dragpos[1]-im[i][1])*(dragpos[1]-im[i][1]);
585          	if (dist<min) {
586          		min=dist;
587          		selected = i;
588          	}		
589          }
590  
591  //        if (editor!=null) editor.select(spoint[selected]);
592  
593          draggraphics = g;
594          g.setXORMode(java.awt.Color.white);
595          oldi = (int)im[selected][0];
596          oldj = (int)im[selected][1];
597  
598                 showStation(g,oldi/dscale,oldj/dscale);
599          return true;
600  }        
601  
602    public void showStation(Graphics g,double x,double y){
603      x *= dscale;
604      y *= dscale;
605      g.drawLine((int)(x+8),(int)(y+4),(int)(x-8),(int)(y+4));
606      g.drawLine((int)(x+8),(int)(y+4),(int)x,(int)(y-8));
607      g.drawLine((int)(x-8),(int)(y+4),(int)x,(int)(y-8));
608    }
609  
610    public void paint(Graphics g,double scale) {
611      dscale = scale;
612      for (int i=0;i<stacnt;i++) {
613        if (i==selected) g.setColor(Segment.selectc);
614        showStation(g,im[i][0],im[i][1]);
615        g.drawString(spoint[i].label,(int)(im[i][0]*dscale+5),(int)(im[i][1]*dscale+5));
616        if (i==selected) g.setColor(Segment.stationc);
617      }
618    }
619  
620    public void solve() {
621      if (solved) return;
622  
623      tree = null;
624  
625      if (isVertical) {
626        if (overlaytrans==null) overlaytrans = AffineTransform.getScaleInstance(1,-1);
627        return;
628      }
629      solved = true;
630  
631      int i,j;
632  
633      double dist = 0;
634  
635  //    System.out.println("Debug data:");
636  //    System.out.println("p1 = "+spoint[p1].label+" at "+spoint[p1].position + "  " +
637  //        	       "p2 = "+spoint[p2].label+" at "+spoint[p2].position);
638  //    System.out.println("dist = "+maxdist);
639  //    System.out.println("image1 = ("+im[p1][0]+","+im[p1][1]+")"+" image2 = ("+im[p2][0]+","+im[p2][1]+")");
640  
641      degenerate = true;
642      double x0 = spoint[wide1].position.X[0];
643      double x1 = spoint[wide2].position.X[0];
644      double y0 = spoint[wide1].position.X[1];
645      double y1 = spoint[wide2].position.X[1];
646      double co = x0*(y1-y0) - y0*(x1-x0);
647  //    double den = Math.sqrt((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1));
648      for (i=0;i<stacnt;i++) {
649        dist = (spoint[i].position.X[0]*(y0-y1) + spoint[i].position.X[1]*(x1-x0) + co)/maxdist;
650  //    System.out.println("spoint["+i+"] = "+spoint[i].label+" at "+spoint[i].position +
651  //        	       " dist = "+dist);
652        if (dist<0) dist = -dist;
653        if (dist>maxdist*linethresh) degenerate = false;
654      }
655  
656      if (degenerate) {
657        spoint[stacnt] = new Vertex(new Vect((x0+x1)/2+(y0-y1),
658          				   (y0+y1)/2+(x1-x0),
659          				   0));
660        im[stacnt][0] = (im[wide1][0]+im[wide2][0])/2 + im[wide2][1] -im[wide1][1];
661        im[stacnt][1] = (im[wide1][1]+im[wide2][1])/2 + im[wide1][0] -im[wide2][0];
662        stacntplus = stacnt+1;
663  //      System.out.println("map point = "+spoint[stacnt].position+" imgpoint = ("+im[stacnt][0]+","+im[stacnt][1]+")");
664      }
665      else stacntplus = stacnt;
666  
667       mapx = new float[stacntplus][2];
668       for (i=0;i<stacntplus;i++) {
669         mapx[i][0] = (float)spoint[i].position.X[0];
670         mapx[i][1] = (float)spoint[i].position.X[1];
671       }
672  
673       Linset set = new Linset(stacntplus+3);
674  //System.out.println("solve setting up matrix");
675       for (i=0;i<stacntplus;i++) {
676         set.A[i][stacntplus] = 1.0;
677         set.A[stacntplus][i] = 1.0;
678         set.A[i][stacntplus+1] =
679             set.A[stacntplus+1][i] = mapx[i][0];
680         set.A[i][stacntplus+2] = 
681             set.A[stacntplus+2][i] = mapx[i][1];
682  
683         for (j=0;j<stacntplus;j++) {
684           set.A[i][j] = h(mapx[j][0]-mapx[i][0],mapx[j][1]-mapx[i][1]);
685           if (i==j) set.A[i][j] = 1.0/beta;
686         }
687       }
688  
689       for (i=stacntplus;i<stacntplus+3;i++)
690         for (j=stacntplus;j<stacntplus+3;j++)
691           set.A[i][j]=0;
692       set.ludcmp();
693       for (i=0;i<stacntplus;i++)
694         set.b[i] = im[i][0];
695       for (i=stacntplus;i<stacntplus+3;i++)
696         set.b[i] = 0;
697  
698       set.lubksb();
699  
700       xcoif = new double[stacntplus+3];
701  
702       for (i=0;i<stacntplus+3;i++)
703         xcoif[i] = set.b[i];
704       for (i=0;i<stacntplus;i++)
705         set.b[i] = im[i][1];
706       for (i=stacntplus;i<stacntplus+3;i++)
707         set.b[i] = 0;
708       
709       set.lubksb();
710  
711       ycoif = new double[stacntplus+3];
712  
713       for (i=0;i<stacntplus+3;i++)
714         ycoif[i] = set.b[i];
715  
716  //System.out.println("solve done");
717  //for (i=0;i<stacnt+3;i++)
718  //System.out.println("xcoif["+i+"] = "+xcoif[i]+" ycoif["+i+"] = "+ycoif[i]);
719    }
720  }
721