View Javadoc
1   package com.github.celldynamics.quimp.geom.filters;
2   
3   import java.util.ArrayList;
4   import java.util.Collections;
5   import java.util.Iterator;
6   import java.util.List;
7   
8   import org.scijava.vecmath.Point2d;
9   import org.slf4j.Logger;
10  import org.slf4j.LoggerFactory;
11  
12  import com.github.celldynamics.quimp.Outline;
13  import com.github.celldynamics.quimp.PointsList;
14  import com.github.celldynamics.quimp.Shape;
15  import com.github.celldynamics.quimp.Snake;
16  import com.github.celldynamics.quimp.Vert;
17  import com.github.celldynamics.quimp.geom.ExtendedVector2d;
18  import com.github.celldynamics.quimp.plugin.ana.ANA_;
19  
20  /**
21   * Support algorithms for processing outlines using their specific properties.
22   * 
23   * @author p.baniukiewicz
24   * @param <T> Outline or Snake class
25   * @see ANA_
26   * @see PointListProcessor
27   */
28  public class OutlineProcessor<T extends Shape<?>> {
29    /**
30     * The Constant LOGGER.
31     */
32    static final Logger LOGGER = LoggerFactory.getLogger(OutlineProcessor.class.getName());
33  
34    private T outline;
35  
36    /**
37     * Assign outline to be processed to object.
38     * 
39     * @param outline Reference to Outline to be processed
40     */
41    public OutlineProcessor(T outline) {
42      this.outline = outline;
43    }
44  
45    /**
46     * Apply running mean filter to Shape.
47     * 
48     * <p>Do not create new Shape but modify nodes of existing one. Compute
49     * {@link Shape#calcCentroid()}, {@link Shape#updateNormals(boolean)} and
50     * {@link Shape#setPositions()}. Normals are updated inwards.
51     * 
52     * @param window size of filter window, must be uneven
53     * @param iters number of smoothing iterations
54     * @return reference to this object. Allows chaining
55     */
56    public OutlineProcessor<T> runningMean(int window, int iters) {
57      PointListProcessorgeom/filters/PointListProcessor.html#PointListProcessor">PointListProcessor pp = new PointListProcessor(outline.asList());
58      pp.smoothMean(window, iters);
59      List<Point2d> p = pp.getList();
60      Iterator<?> it = outline.iterator();
61      Iterator<Point2d> itl = p.iterator();
62      while (it.hasNext() && itl.hasNext()) {
63        PointsList<?> n = (PointsList<?>) it.next();
64        ExtendedVector2d v = n.getPoint();
65        Point2d pr = itl.next();
66        v.x = pr.x;
67        v.y = pr.y;
68      }
69      outline.calcCentroid();
70      outline.updateNormals(true);
71      outline.setPositions();
72      return this;
73    }
74  
75    /**
76     * Apply running median filter to Shape.
77     * 
78     * <p>Do not create new Shape but modify nodes of existing one. Compute
79     * {@link Shape#calcCentroid()}, {@link Shape#updateNormals(boolean)} and
80     * {@link Shape#setPositions()}. Normals are updated inwards.
81     * 
82     * @param window size of filter window, must be uneven
83     * @param iters number of smoothing iterations
84     * @return reference to this object. Allows chaining
85     */
86    public OutlineProcessor<T> runningMedian(int window, int iters) {
87      PointListProcessorgeom/filters/PointListProcessor.html#PointListProcessor">PointListProcessor pp = new PointListProcessor(outline.asList());
88      pp.smoothMedian(window, iters);
89      List<Point2d> p = pp.getList();
90      Iterator<?> it = outline.iterator();
91      Iterator<Point2d> itl = p.iterator();
92      while (it.hasNext() && itl.hasNext()) {
93        PointsList<?> n = (PointsList<?>) it.next();
94        ExtendedVector2d v = n.getPoint();
95        Point2d pr = itl.next();
96        v.x = pr.x;
97        v.y = pr.y;
98      }
99      outline.calcCentroid();
100     outline.updateNormals(true);
101     outline.setPositions();
102     return this;
103   }
104 
105   /**
106    * Compute average curvature of {@link Outline}.
107    * 
108    * <p>Does not apply to {@link Snake}. Require {@link Outline#updateCurvature()} to be called
109    * first to fill {@link Vert#curvatureLocal} field.
110    * 
111    * @param averdistance Average distance in pixels to compute curvature. Nodes in half of distance
112    *        on left and right from current node will be included in average curvature for this node.
113    * @return instance of this object.
114    */
115   public OutlineProcessor<T> averageCurvature(double averdistance) {
116     Vert v;
117     Vert tmpV;
118     double totalCur;
119     double distance;
120     int count = 0;
121     if (!(outline instanceof Outline)) {
122       throw new IllegalArgumentException("This method applies to Outline only");
123     }
124     // Average over curvatures
125     if (averdistance > 0) {
126       // System.out.println("new outline");
127       v = (Vert) outline.getHead();
128       do {
129         // System.out.println("\tnew vert");
130         totalCur = v.curvatureLocal; // reset
131         count = 1;
132 
133         // add up curvatures of prev nodes
134         // System.out.println("\t prev nodes");
135         tmpV = v.getPrev();
136         distance = 0;
137         do {
138           distance += ExtendedVector2d.lengthP2P(tmpV.getNext().getPoint(), tmpV.getPoint());
139           totalCur += tmpV.curvatureLocal;
140           count++;
141           tmpV = tmpV.getPrev();
142         } while (distance < averdistance / 2);
143 
144         // add up curvatures of next nodes
145         distance = 0;
146         tmpV = v.getNext();
147         do {
148           distance += ExtendedVector2d.lengthP2P(tmpV.getPrev().getPoint(), tmpV.getPoint());
149           totalCur += tmpV.curvatureLocal;
150           count++;
151           tmpV = tmpV.getNext();
152         } while (distance < averdistance / 2);
153 
154         v.curvatureSmoothed = totalCur / count;
155 
156         v = v.getNext();
157       } while (!v.isHead());
158     }
159     LOGGER.trace("Average curvature over:" + averdistance + " nodes number:" + count);
160     return this;
161   }
162 
163   /**
164    * Sum smoothed curvature over a region of the membrane. For {@link Outline} only.
165    * 
166    * <p><p>Does not apply to {@link Snake}. Require {@link #averageCurvature(double)} to be called
167    * first to fill {@link Vert#curvatureSmoothed} field.
168    * 
169    * @param averdistance Average distance in pixels to compute curvature. Nodes in half of distance
170    *        on left and right from current node will be included in average curvature for this node.
171    * @return instance of this object.
172    */
173   public OutlineProcessor<T> sumCurvature(double averdistance) {
174     Vert v;
175     Vert tmpV;
176     double totalCur;
177     double distance;
178     if (!(outline instanceof Outline)) {
179       throw new IllegalArgumentException("This method applies to Outline only");
180     }
181     // Average over curvatures
182     if (averdistance > 0) {
183       v = (Vert) outline.getHead();
184       do {
185         // System.out.println("\tnew vert");
186         totalCur = v.curvatureSmoothed; // reset
187         // add up curvatures of prev nodes
188         // System.out.println("\t prev nodes");
189         tmpV = v.getPrev();
190         distance = 0;
191         do {
192           distance += ExtendedVector2d.lengthP2P(tmpV.getNext().getPoint(), tmpV.getPoint());
193           totalCur += tmpV.curvatureSmoothed;
194           tmpV = tmpV.getPrev();
195         } while (distance < averdistance / 2);
196 
197         // add up curvatures of next nodes
198         distance = 0;
199         tmpV = v.getNext();
200         do {
201           distance += ExtendedVector2d.lengthP2P(tmpV.getPrev().getPoint(), tmpV.getPoint());
202           totalCur += tmpV.curvatureSmoothed;
203           tmpV = tmpV.getNext();
204         } while (distance < averdistance / 2);
205 
206         v.curvatureSum = totalCur;
207 
208         v = v.getNext();
209       } while (!v.isHead());
210     }
211     return this;
212   }
213 
214   /**
215    * Shrink outline linearly.
216    * 
217    * @param steps steps
218    * @param stepRes stepRes
219    * @param angleTh angleTh
220    * @param freezeTh freezeTh
221    * @see #shrinknl(double, double, double, double, double, double, double, double)
222    */
223   public void shrinkLin(double steps, double stepRes, double angleTh, double freezeTh) {
224     final double defSigma = 0.3;
225     final double defAmpl = 1.0; // linear shrink
226     final double defCurvAverDist = 1.0;
227     final double defNormalsDist = 0; // off
228     shrinknl(steps, stepRes, angleTh, freezeTh, defSigma, defAmpl, defCurvAverDist, defNormalsDist);
229   }
230 
231   /**
232    * Shrink the outline nonlinearly.
233    * 
234    * <p>This method uses {@link Vert#curvatureSum} computed by
235    * {@link OutlineProcessor#sumCurvature(double)} to produce nonlinear scaling. Scaling factor
236    * depends on {@link Vert#curvatureSum} and it is computed by
237    * {@link #amplificationFactor(double, double, double)} method. Nodes are shifted inwards
238    * according to normals ({@link Vert#getNormal()}). If parameter <i>averageNormalsDist</i> is
239    * greater than zero, nodes within this distance to current node will be analysed to find lowest
240    * (must be negative) {@link Vert#curvatureSum} and then normal vector of this node will be copied
241    * to all other nodes
242    * in current window.
243    * 
244    * <p>{@link #amplificationFactor(double, double, double)} increase step size used by
245    * {@link #shrinkLin(double, double, double, double)} maximally by factor of <t>magn</t>
246    * ({@link #amplificationFactor(double, double, double)} option).
247    * 
248    * <p>For linear shrinking use <t>magn</t>=1 and <t>averageNormalsDist</t>=0.
249    * 
250    * @param steps number of steps - integer
251    * @param stepRes length of the step
252    * @param angleTh angle threshold
253    * @param freezeTh freeze threshold
254    * @param sigma sigma of Gaussian
255    * @param magn maximum amplification (for curv<<0)
256    * @param averageCurvDist number of neighbouring nodes to average curvature over. This is
257    *        multiplier of average node distance. Always at least three nodes are averaged
258    *        (n-1,n,n+1).
259    * @param averageNormalsDist number of neighbouring nodes to set normals the same (set 0 to
260    *        disable).
261    * @see Outline#scaleOutline(double, double, double, double)
262    * @see OutlineProcessor#amplificationFactor(double, double, double)
263    */
264   public void shrinknl(double steps, double stepRes, double angleTh, double freezeTh, double sigma,
265           double magn, double averageCurvDist, double averageNormalsDist) {
266     if (!(outline instanceof Outline)) {
267       throw new IllegalArgumentException("This method applies to Outline only");
268     }
269     LOGGER.debug("Steps: " + steps);
270     LOGGER.debug("Original res: " + outline.getNumPoints());
271 
272     Vert n;
273     int max = 10000;
274     double d = outline.getLength() / outline.getNumPoints();
275     LOGGER.debug("steps:" + steps + " stepRes:" + stepRes + " angleTh:" + angleTh + " freezeTh:"
276             + freezeTh + " sigma:" + sigma + " magn:" + magn + " averageCurvDist:" + averageCurvDist
277             + " averageNormalsDist" + averageNormalsDist + " averDistanceBefore:" + d);
278     ((Outline) outline).correctDensity(d, d / 2);
279     ((Outline) outline).updateNormals(true);
280     ((Outline) outline).updateCurvature();
281     double d1 = outline.getLength() / outline.getNumPoints();
282     averageCurvature(d1 * averageCurvDist).sumCurvature(d1 * averageCurvDist);
283     constrainNormals(d1 * averageNormalsDist);
284 
285     LOGGER.debug("Res after 1: " + outline.getNumPoints());
286 
287     for (int j = 0; j < steps; j++) {
288       n = (Vert) outline.getHead();
289       do {
290         if (!n.isFrozen()) {
291           n.setX(n.getX() - stepRes * amplificationFactor(n.curvatureSum, sigma, magn)
292                   * n.getNormal().getX());
293           n.setY(n.getY() - stepRes * amplificationFactor(n.curvatureSum, sigma, magn)
294                   * n.getNormal().getY());
295         }
296         n = n.getNext();
297       } while (!n.isHead());
298 
299       ((Outline) outline).removeProx(1.5, 1.5);
300       ((Outline) outline).freezeProx(angleTh, freezeTh);
301       ((Outline) outline).correctDensity(d, d / 2);
302       // ((Outline) outline).updateNormals(true); // sometimes it is better
303       ((Outline) outline).updateCurvature();
304       d1 = outline.getLength() / outline.getNumPoints();
305       averageCurvature(d1 * averageCurvDist).sumCurvature(d1 * averageCurvDist);
306       // equalNormales(d1 * 5); // better if original normals stay. interpolate can add vertex but
307       // ith will have updated normales
308 
309       // do not shrink if there are 4 nodes or less
310       if (outline.getNumPoints() <= 4) {
311         LOGGER.debug("Stopped iterations");
312         break;
313       }
314 
315       if (j > max) {
316         LOGGER.warn("shrink (336) hit max iterations!");
317         break;
318       }
319     }
320 
321     if (outline.getNumPoints() < 3) {
322       LOGGER.info("ANA 377_NODES LESS THAN 3 BEFORE CUTS");
323     }
324 
325     if (((Outline) outline).cutSelfIntersects()) {
326       LOGGER.debug("ANA_(382)...fixed ana intersects");
327     }
328 
329     if (outline.getNumPoints() < 3) {
330       LOGGER.info("ANA 377_NODES LESS THAN 3");
331     }
332   }
333 
334   /**
335    * Set normals to the same direction according to local minimum of sumCurvature.
336    * 
337    * @param averdistance distance to look along. Set 0 to turn off any modification. Any value
338    *        smaller than distance will cause that at least 3 nodes are taken.
339    */
340   private void constrainNormals(double averdistance) {
341     Vert v;
342     Vert tmpV;
343     double distance;
344     if (!(outline instanceof Outline)) {
345       throw new IllegalArgumentException("This method applies to Outline only");
346     }
347     LOGGER.debug("constrainNormals over:" + averdistance);
348     if (averdistance == 0) {
349       return;
350     }
351     List<Double> curvSums = new ArrayList<>(); // will keep local curvature sums in +-distance
352     List<ExtendedVector2d> norms = new ArrayList<>(); // will keep normals in +-distance
353 
354     if (averdistance > 0) {
355       v = (Vert) outline.getHead();
356       do { // main loop
357         curvSums.clear(); // new vertex
358         norms.clear();
359         curvSums.add(v.curvatureSum); // remember it
360         norms.add(v.getNormal());
361         tmpV = v.getPrev();
362         distance = 0;
363         do { // after current vert in +distance
364           distance += ExtendedVector2d.lengthP2P(tmpV.getNext().getPoint(), tmpV.getPoint());
365           curvSums.add(tmpV.curvatureSum); // Remember all verts in this window
366           norms.add(tmpV.getNormal());
367           tmpV = tmpV.getPrev();
368         } while (distance < averdistance / 2);
369         distance = 0;
370         tmpV = v.getNext();
371         do { // before current vert
372           distance += ExtendedVector2d.lengthP2P(tmpV.getPrev().getPoint(), tmpV.getPoint());
373           curvSums.add(tmpV.curvatureSum); // Remember all verts in this window
374           norms.add(tmpV.getNormal());
375           tmpV = tmpV.getNext();
376         } while (distance < averdistance / 2);
377         // v is current vert
378         // find min curv
379         Double minCurv = Collections.min(curvSums);
380         int indmin = curvSums.indexOf(minCurv); // min curvature index
381         if (minCurv < 0) { // propagate only in curvature is negative
382           ExtendedVector2d normmin = norms.get(indmin); // get its normale
383           v.setNormal(normmin.getX(), normmin.getY()); // and propagate to current
384         }
385         v = v.getNext(); // go to next vert
386       } while (!v.isHead());
387     }
388   }
389 
390   /**
391    * Return amplification factor dependent on curvature.
392    * 
393    * <p>For positive curvature it returns 1, for negative a reversed Gaussian curve scaled between
394    * 1.0 and specified magnitude is returned.
395    * 
396    * <pre>
397    * <code>
398    * x=linspace(-2,0,100);
399    * sig = 0.3;
400    * u = 0;
401    * a = 4
402    * y = (a*(1-exp(-0.5*((x-u)/sig).^2)))/(a/(a-1))+1;
403    * plot(x,y)
404    * </code>
405    * </pre>
406    * 
407    * @param curv curvature value to compute magnification coefficient for
408    * @param sigma sigma of Gaussian
409    * @param magn maximum amplification (for curv<<0). Set to 1.0 for constant output of 1.0
410    * @return amplification coefficient
411    */
412   public double amplificationFactor(double curv, double sigma, double magn) {
413     double ampl;
414     double u = 0;
415     if (magn == 1) {
416       return 1.0;
417     }
418     if (curv >= 0) {
419       ampl = 1.0;
420     } else {
421       ampl = (magn * (1 - Math.exp(-0.5 * (Math.pow((curv - u) / sigma, 2))))) / (magn / (magn - 1))
422               + 1;
423     }
424     return ampl;
425   }
426 
427   /**
428    * Outline getter.
429    * 
430    * @return the outline
431    */
432   public T getO() {
433     return outline;
434   }
435 
436 }