View Javadoc

1   package org.textensor.stochdiff.disc;
2   
3   import java.util.ArrayList;
4   
5   import org.textensor.stochdiff.numeric.morph.TreePoint;
6   
7   /* This will remove points that divide segments into too many
8    * sections, subject to criteria on the allowed change in radius and
9    * the maximum length.
10   */
11  
12  
13  
14  public class SegmentMerger {
15  
16      private TreePoint[] merge(TreePoint[] pt, double maxdr, double maxlen, double tol) {
17  
18  
19          for (int i = 0; i < pt.length; i++) {
20              pt[i].setWork(2);
21          }
22          TreePoint p0 = pt[0];
23          recMerge(p0, maxdr, maxlen, tol);
24  
25          // remove dead points.
26          int nl = 0;
27          for (int i = 0; i < pt.length; i++) {
28              if (pt[i].getWork() > 0)
29                  nl++;
30          }
31          TreePoint[] pr = new TreePoint[nl];
32          nl = 0;
33          for (int i = 0; i < pt.length; i++) {
34              if (pt[i].getWork() > 0)
35                  pr[nl++] = pt[i];
36          }
37          return pr;
38      }
39  
40  
41      private void recMerge(TreePoint cp, double maxdr, double maxlen, double tol) {
42          // may have already been done
43          cp.setWork(1);
44  
45          for (int j = 0; j < cp.nnbr; j++) {
46              TreePoint cq = cp.nbr[j];
47              if (cq.getWork() == 2 && cq.nnbr == 2
48                      && Math.abs((cq.r - cp.r) / (cq.r + cp.r)) < 0.5 * maxdr
49                      && cp.distanceTo(cq) < maxlen) {
50                  // if nnbr is not two, either it is a terminal, or a branch point:
51                  // in either case there is nothing to do;
52  
53                  TreePoint cprev = cp;
54                  ArrayList<TreePoint> vpt = new ArrayList<TreePoint>();
55                  double ltot = 0.;
56                  double ldtot = 0.;
57                  while (cq.nnbr == 2 && cprev.distanceTo(cq) < maxlen
58                          && Math.abs((cq.r - cp.r) / (cq.r + cp.r)) < 0.5 * maxdr) {
59                      vpt.add(cq);
60                      double dl = cprev.distanceTo(cq);
61                      ltot += dl;
62                      ldtot += dl * (cprev.r + cq.r);
63  
64                      TreePoint cnxt = (cq.nbr[0] == cprev ? cq.nbr[1] : cq.nbr[0]);
65                      cprev = cq;
66                      cq = cnxt;
67                  }
68                  double dl = cprev.distanceTo(cq);
69                  ltot += dl;
70                  ldtot += dl * (cprev.r + cq.r);
71  
72  
73                  double lab = cp.distanceTo(cq);
74                  double ldab = lab * (cp.r + cq.r);
75                  int nadd = (int)(ltot / maxlen);
76                  if (nadd > vpt.size())
77                      nadd = vpt.size(); // cant happen at present;
78                  // recycle nadd points;
79  
80                  boolean cor = (Math.abs((lab - ltot) / (lab + ltot)) > 0.5 * tol || Math.abs((ldab - ldtot)
81                                 / (ldab + ldtot)) > 0.5 * tol);
82                  if (cor && nadd == 0)
83                      nadd = 1;
84  
85                  if (nadd == 0) {
86                      cp.replaceNeighbor(vpt.get(0), cq);
87                      cq.replaceNeighbor(vpt.get(vpt.size() - 1), cp);
88                  } else {
89                      for (int i = 0; i < nadd; i++) {
90                          TreePoint cm = vpt.get(i);
91                          cm.setWork(1);
92                          cm.locateBetween(cp, cq, (1. + i) / (1. + nadd));
93                          if (i == nadd - 1 && nadd < vpt.size()) {
94                              cm.replaceNeighbor(vpt.get(nadd), cq);
95                              cq.replaceNeighbor(vpt.get(vpt.size() - 1), cm);
96                          }
97                      }
98                  }
99  
100 
101                 // just kill the rest;
102                 for (int jd = nadd; jd < vpt.size(); jd++) {
103                     TreePoint cd = vpt.get(jd);
104                     cd.nnbr = 0;
105                     cd.setWork(0);
106                 }
107 
108 
109                 if (cor) {
110                     double dpar = lab / (nadd + 1);
111                     double fl = ltot / lab;
112                     double dperp = Math.sqrt((fl * fl - 1) * dpar * dpar);
113                     double fr = ldtot / (fl * ldab);
114                     for (int i = 0; i < nadd; i++) {
115                         TreePoint cm = (vpt.get(i));
116                         cm.r *= fr;
117                         if (i % 2 == 0) {
118                             cm.movePerp(cp, cq, ((i / 2) % 2 == 0 ? dperp : -dperp));
119                         }
120                     }
121 
122                 }
123 
124 
125 
126             }
127             if (cq.getWork() == 2)
128                 recMerge(cq, maxdr, maxlen, tol);
129         }
130     }
131 
132 
133 }