?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Use internal SNAPSHOT couplings again
[trakem2.git] / TrakEM2_ / src / main / java / ini / trakem2 / utils / M.java
blob95abd5a67834bd4253ab693b12579202a63750cd
1 package ini.trakem2.utils;
2
3 import ij.gui.PolygonRoi;
4 import ij.gui.Roi;
5 import ij.gui.ShapeRoi;
6 import ij.process.FloatPolygon;
7 import ini.trakem2.display.VectorDataTransform;
8 import ini.trakem2.persistence.Loader;
9 import ini.trakem2.vector.VectorString2D;
10
11 import java.awt.Color;
12 import java.awt.Graphics2D;
13 import java.awt.Polygon;
14 import java.awt.Rectangle;
15 import java.awt.Shape;
16 import java.awt.geom.AffineTransform;
17 import java.awt.geom.Area;
18 import java.awt.geom.GeneralPath;
19 import java.awt.geom.NoninvertibleTransformException;
20 import java.awt.geom.PathIterator;
21 import java.awt.geom.Point2D;
22 import java.awt.image.BufferedImage;
23 import java.awt.image.DataBufferByte;
24 import java.util.ArrayList;
25 import java.util.Collection;
26 import java.util.List;
27
28 import javax.vecmath.Point3f;
29 import javax.vecmath.Tuple3d;
30 import javax.vecmath.Vector3d;
31 import javax.vecmath.Vector3f;
32
33 /** TrakEM2's mathematician. */
34 public final class M {
35
36 /*============== 3D =============*/
37
38 static public final double distancePointToSegment
39 (final double px, final double py, final double pz,
40 final double lx1, final double ly1, final double lz1,
41 final double lx2, final double ly2, final double lz2)
42 {
43 return distancePointToSegment(new Vector3d(px, py, pz),
44 new Vector3d(lx1, ly1, lz1),
45 new Vector3d(lx2, ly2, lz2));
46 }
47 static public final double distancePointToSegmentSq
48 (final double px, final double py, final double pz,
49 final double lx1, final double ly1, final double lz1,
50 final double lx2, final double ly2, final double lz2)
51 {
52 return distancePointToSegmentSq(new Vector3d(px, py, pz),
53 new Vector3d(lx1, ly1, lz1),
54 new Vector3d(lx2, ly2, lz2));
55 }
56
57 /** Minimum distance between point v0 and a line segment defined by points v1 and v2. */
58 static public final double distancePointToSegment(final Vector3d p, final Vector3d v1, final Vector3d v2) {
59 final Vector3d v = new Vector3d();
60 v.sub(v2, v1);
61 final Vector3d w = new Vector3d();
62 w.sub(p, v1);
63
64 final double c1 = w.dot(v);
65 if (c1 <= 0) return distance(p, v1);
66
67 final double c2 = v.dot(v);
68 if (c2 <= c1) return distance(p, v2);
69
70 final double b = c1 / c2;
71
72 final Vector3d pb = new Vector3d(v);
73 pb.scale(b);
74 pb.add(v1);
75 return distance(p, pb);
76 }
77
78 static public final double distancePointToSegmentSq(final Vector3d p, final Vector3d v1, final Vector3d v2) {
79 final Vector3d v = new Vector3d();
80 v.sub(v2, v1);
81 final Vector3d w = new Vector3d();
82 w.sub(p, v1);
83
84 final double c1 = w.dot(v);
85 if (c1 <= 0) return distanceSq(p, v1);
86
87 final double c2 = v.dot(v);
88 if (c2 <= c1) return distanceSq(p, v2);
89
90 final double b = c1 / c2;
91
92 final Vector3d pb = new Vector3d(v);
93 pb.scale(b);
94 pb.add(v1);
95 return distanceSq(p, pb);
96 }
97
98
99 static public final double distance(final Tuple3d t1, final Tuple3d t2) {
100 return Math.sqrt( Math.pow(t1.x - t2.x, 2)
101 + Math.pow(t1.y - t2.y, 2)
102 + Math.pow(t1.z - t2.z, 2) );
103 }
104
105 static public final double distanceSq(final Tuple3d t1, final Tuple3d t2) {
106 return Math.pow(t1.x - t2.x, 2)
107 + Math.pow(t1.y - t2.y, 2)
108 + Math.pow(t1.z - t2.z, 2);
109 }
110
111 static public final double distance(final double x1, final double y1, final double z1,
112 final double x2, final double y2, final double z2) {
113 return Math.sqrt(Math.pow(x1 - x2, 2)
114 + Math.pow(y1 - y2, 2)
115 + Math.pow(z1 - z2, 2));
116 }
117
118 static public final double distance(final double x1, final double y1,
119 final double x2, final double y2) {
120 return Math.sqrt(Math.pow(x1 - x2, 2)
121 + Math.pow(y1 - y2, 2));
122 }
123
124 static public final double distanceSq(final double x1, final double y1, final double z1,
125 final double x2, final double y2, final double z2) {
126 return Math.pow(x1 - x2, 2)
127 + Math.pow(y1 - y2, 2)
128 + Math.pow(z1 - z2, 2);
129 }
130
131 /** In 3D */
132 static public double distancePointToLine(final double px, final double py, final double pz, final double lx1, final double ly1, final double lz1, final double lx2, final double ly2, final double lz2 ) {
133 final double segment_length = new Vector3d(lx2 - lx1, ly2 - ly1, lz2 - lz1).length();
134 if (0 == segment_length) return 0;
135 final Vector3d cross = new Vector3d();
136 cross.cross(new Vector3d(px - lx1, py - ly1, pz - lz1),
137 new Vector3d(px - lx2, py - ly2, pz - lz2));
138 return cross.length() / segment_length;
139 }
140
141
142 /*============== 2D =============*/
143
144 /** For a point px,py return its distance to a line defined by points lx1,ly1 and lx2,ly2. Formula as in http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html at 2008-11-11 20:19 Zurich time.*/
145 static public double distancePointToLine(final double px, final double py, final double lx1, final double ly1, final double lx2, final double ly2) {
146 return Math.abs( (lx2 - lx1) * (ly1 - py) - (lx1 - px) * (ly2 - ly1) )
147 / Math.sqrt( Math.pow(lx2 - lx1, 2) + Math.pow(ly2 - ly1, 2) );
148 }
149
150
151 /** For a point px,py return its distance to a line SEGMENT defined by points lx1,ly1 and lx2,ly2. */
152 static public double distancePointToSegment(final double px, final double py, final double lx1, final double ly1, final double lx2, final double ly2) {
153 double xu = px - lx1;
154 double yu = py - ly1;
155 double xv = lx2 - lx1;
156 double yv = ly2 - ly1;
157 if (xu * xv + yu * yv < 0) return Math.sqrt( Math.pow(px - lx1, 2) + Math.pow(py - ly1, 2) );
158
159 xu = px - lx2;
160 yu = py - ly2;
161 xv = -xv;
162 yv = -yv;
163 if (xu * xv + yu * yv < 0) return Math.sqrt( Math.pow(px - lx2, 2) + Math.pow(py - ly2, 2) );
164
165 return Math.abs( (px * (ly1 - ly2) + py * (lx2 - lx1) + (lx1 * ly2 - lx2 * ly1) ) / Math.sqrt( Math.pow(lx2 - lx1, 2) + Math.pow(ly2 - ly1, 2) ) );
166 }
167
168 static public final boolean equals(final double a, final double b) {
169 return Math.abs(a - b) < Utils.FL_ERROR;
170 }
171
172 static public final Point2D.Double transform(final AffineTransform affine, final double x, final double y) {
173 final Point2D.Double pSrc = new Point2D.Double(x, y);
174 if (affine.isIdentity()) return pSrc;
175 final Point2D.Double pDst = new Point2D.Double();
176 affine.transform(pSrc, pDst);
177 return pDst;
178 }
179
180 static public final Point2D.Double inverseTransform(final AffineTransform affine, final double x, final double y) {
181 final Point2D.Double pSrc = new Point2D.Double(x, y);
182 if (affine.isIdentity()) return pSrc;
183 final Point2D.Double pDst = new Point2D.Double();
184 try {
185 affine.createInverse().transform(pSrc, pDst);
186 } catch (NoninvertibleTransformException nite) {
187 IJError.print(nite);
188 }
189 return pDst;
190 }
191
192 // from utilities.c in my CurveMorphing C module ... from C! Java is a low level language with the disadvantages of the high level languages ...
193 /** Returns the angle in radians of the given polar coordinates, correcting the Math.atan2 output.
194 * Adjusting so that 0 is 3 o'clock, PI+PI/2 is 12 o'clock, PI is 9 o'clock, and PI/2 is 6 o'clock (why atan2 doesn't output angles this way? I remember I had the same problem for Pipe.java in the A_3D_editing plugin) */
195 static public final double getAngle(final double x, final double y) {
196 // calculate angle
197 double a = Math.atan2(x, y);
198 // fix too large angles (beats me why are they ever generated)
199 if (a > 2 * Math.PI) {
200 a = a - 2 * Math.PI;
201 }
202 // fix atan2 output scheme to match my mental scheme
203 if (a >= 0.0 && a <= Math.PI/2) {
204 a = Math.PI/2 - a;
205 } else if (a < 0 && a >= -Math.PI) {
206 a = Math.PI/2 -a;
207 } else if (a > Math.PI/2 && a <= Math.PI) {
208 a = Math.PI + Math.PI + Math.PI/2 - a;
209 }
210 return a;
211 }
212
213
214 /*============== Geometry =============*/
215
216 static public final boolean isEmpty(final Area area) {
217 final Rectangle b = area.getBounds();
218 return 0 == b.width || 0 == b.height;
219 }
220
221 /** Test whether the areas intersect each other. */
222 static public final boolean intersects(final Area a1, final Area a2) {
223 final Area b = new Area(a1);
224 b.intersect(a2);
225 final java.awt.Rectangle r = b.getBounds();
226 return 0 != r.width && 0 != r.height;
227 }
228
229 static public final Area getArea(final Roi roi) {
230 if (null == roi) return null;
231 if (roi instanceof ShapeRoi) return getArea((ShapeRoi)roi);
232 return getArea(new ShapeRoi(roi));
233 }
234 static public final Area getArea(final ShapeRoi sroi) {
235 if (null == sroi) return null;
236 AffineTransform at = new AffineTransform();
237 Rectangle bounds = sroi.getBounds();
238 at.translate(bounds.x, bounds.y);
239 Area area = new Area(sroi.getShape());
240 area.transform(at);
241 return area;
242 }
243
244 /** Returns the approximated area of the given Area object; Loader can be null; if not, it's used to secure memory space. */
245 static public final double measureArea(Area area, final Loader loader) {
246 double sum = 0;
247 try {
248 Rectangle bounds = area.getBounds();
249 double scale = 1;
250 if (bounds.width > 2048 || bounds.height > 2048) { // TODO value 2048 should be reconsidered as a project property
251 scale = 2048.0 / bounds.width;
252 }
253 if (0 == scale) {
254 Utils.log("Can't measure: area too large, out of scale range for approximation.");
255 return sum;
256 }
257 AffineTransform at = new AffineTransform();
258 at.translate(-bounds.x, -bounds.y);
259 at.scale(scale, scale);
260 area = area.createTransformedArea(at);
261 bounds = area.getBounds();
262 if (0 == bounds.width || 0 == bounds.height) {
263 Utils.log("Can't measure: area too large, approximates zero.");
264 return sum;
265 }
266 if (null != loader) loader.releaseToFit(bounds.width * bounds.height * 3);
267 BufferedImage bi = new BufferedImage(bounds.width, bounds.height, BufferedImage.TYPE_BYTE_INDEXED);
268 Graphics2D g = bi.createGraphics();
269 g.setColor(Color.white);
270 g.fill(area);
271 final byte[] pixels = ((DataBufferByte)bi.getRaster().getDataBuffer()).getData(); // buffer.getData();
272 for (int i=pixels.length-1; i>-1; i--) {
273 //if (255 == (pixels[i]&0xff)) sum++;
274 if (0 != pixels[i]) sum++;
275 }
276 bi.flush();
277 g.dispose();
278 if (1 != scale) sum = sum / (scale * scale);
279 } catch (Throwable e) {
280 IJError.print(e);
281 }
282 return sum;
283 }
284
285 /** Compute the area of the triangle defined by 3 points in 3D space, returning half of the length of the vector resulting from the cross product of vectors p1p2 and p1p3. */
286 static public final double measureArea(final Point3f p1, final Point3f p2, final Point3f p3) {
287 // Distance from p1 to line p2-p3, times length of line p2-p3, divided by 2:
288 return 0.5 * M.distancePointToLine(p1.x, p1.y, p1.z,
289 p2.x, p2.y, p2.z,
290 p3.x, p3.y, p3.z)
291 * p2.distance(p3);
292 }
293
294 /** Returns true if the roi is of closed shape type like an OvalRoi, ShapeRoi, a Roi rectangle, etc. */
295 static public final boolean isAreaROI(final Roi roi) {
296 switch (roi.getType()) {
297 case Roi.POLYLINE:
298 case Roi.FREELINE:
299 case Roi.LINE:
300 case Roi.POINT:
301 return false;
302 }
303 return true;
304 }
305
306 @SuppressWarnings("null")
307 static public final Collection<Polygon> getPolygons(Area area) {
308 final ArrayList<Polygon> pols = new ArrayList<Polygon>();
309 Polygon pol = null;
310
311 final float[] coords = new float[6];
312 for (PathIterator pit = area.getPathIterator(null); !pit.isDone(); ) {
313 int seg_type = pit.currentSegment(coords);
314 switch (seg_type) {
315 case PathIterator.SEG_MOVETO:
316 pol = new Polygon();
317 //$FALL-THROUGH$
318 case PathIterator.SEG_LINETO:
319 pol.addPoint((int)coords[0], (int)coords[1]);
320 break;
321 case PathIterator.SEG_CLOSE:
322 pols.add(pol);
323 break;
324 default:
325 Utils.log2("WARNING: unhandled seg type.");
326 break;
327 }
328 pit.next();
329 if (pit.isDone()) {
330 break;
331 }
332 }
333
334 return pols;
335 }
336 static public final Collection<Polygon> getPolygonsByRounding(final Area area) {
337 final ArrayList<Polygon> pols = new ArrayList<Polygon>();
338 Polygon pol = new Polygon();
339
340 final float[] coords = new float[6];
341 for (PathIterator pit = area.getPathIterator(null); !pit.isDone(); ) {
342 int seg_type = pit.currentSegment(coords);
343 switch (seg_type) {
344 case PathIterator.SEG_MOVETO:
345 case PathIterator.SEG_LINETO:
346 pol.addPoint(Math.round(coords[0]), Math.round(coords[1]));
347 break;
348 case PathIterator.SEG_CLOSE:
349 pols.add(pol);
350 pol = new Polygon();
351 break;
352 default:
353 Utils.log2("WARNING: unhandled seg type.");
354 break;
355 }
356 pit.next();
357 if (pit.isDone()) {
358 break;
359 }
360 }
361 return pols;
362 }
363
364 /** Return a new Area resulting from applying @param ict to @param a;
365 * assumes areas consists of paths with moveTo, lineTo and close operations. */
366 static public final Area transform(final mpicbg.models.CoordinateTransform ict, final Area a) {
367 final GeneralPath path = new GeneralPath();
368 final float[] coords = new float[6];
369 final float[] fp = new float[2];
370
371 for (final PathIterator pit = a.getPathIterator(null); !pit.isDone(); ) {
372 final int seg_type = pit.currentSegment(coords);
373 fp[0] = coords[0];
374 fp[1] = coords[1];
375 ict.applyInPlace(fp);
376 switch (seg_type) {
377 case PathIterator.SEG_MOVETO:
378 path.moveTo(fp[0], fp[1]);
379 break;
380 case PathIterator.SEG_LINETO:
381 case PathIterator.SEG_CLOSE:
382 path.lineTo(fp[0], fp[1]);
383 break;
384 default:
385 Utils.log2("WARNING: unhandled seg type.");
386 break;
387 }
388 pit.next();
389 if (pit.isDone()) {
390 break;
391 }
392 }
393 return new Area(path);
394 }
395
396 /** Apply in place the @param ict to the Area @param a, but only for the part that intersects the roi. */
397 static final public void apply(final mpicbg.models.CoordinateTransform ict, final Area roi, final Area a) {
398 final Area intersection = new Area(a);
399 intersection.intersect(roi);
400 a.subtract(intersection);
401 a.add(M.transform(ict, intersection));
402 }
403
404 static final public void apply(final VectorDataTransform vdt, final Area a) {
405 apply(vdt, a, false);
406 }
407
408 /** Parts of @param a not intersected by any of @param vdt rois will be left untouched if @param remove_outside is false. */
409 static final public void apply(final VectorDataTransform vdt, final Area a, final boolean remove_outside) {
410 final Area b = new Area();
411 for (final VectorDataTransform.ROITransform rt : vdt.transforms) {
412 // Cut the intersecting part from a:
413 final Area intersection = new Area(a);
414 intersection.intersect(rt.roi);
415 a.subtract(intersection);
416 // .. and add it to b, transformed:
417 b.add(M.transform(rt.ct, intersection));
418 }
419
420 if (!M.isEmpty(a)) {
421 if (remove_outside) {
422 // Clear areas not affected any ROITransform
423 Utils.log("WARNING: parts of an area in layer " + vdt.layer + "\n did not intersect any transformation target\n and where removed.");
424 a.reset();
425 } else Utils.log("WARNING: parts of an area in layer " + vdt.layer + "\n remain untransformed.");
426 }
427
428 // Add b (the transformed parts) to what remains of a
429 a.add(b);
430 }
431
432 /** Detect if a point is not in the area, but lays inside one of its path, which is returned as a Polygon. Otherwise returns null. The given x,y must be already in the Area's coordinate system. */
433 static public final Polygon findPath(final Area area, final int x, final int y) {
434 Polygon pol = new Polygon();
435 for (PathIterator pit = area.getPathIterator(null); !pit.isDone(); ) {
436 float[] coords = new float[6];
437 int seg_type = pit.currentSegment(coords);
438 switch (seg_type) {
439 case PathIterator.SEG_MOVETO:
440 case PathIterator.SEG_LINETO:
441 pol.addPoint((int)coords[0], (int)coords[1]);
442 break;
443 case PathIterator.SEG_CLOSE:
444 if (pol.contains(x, y)) return pol;
445 // else check next
446 pol = new Polygon();
447 break;
448 default:
449 Utils.log2("WARNING: unhandled seg type.");
450 break;
451 }
452 pit.next();
453 if (pit.isDone()) {
454 break;
455 }
456 }
457 return null;
458 }
459
460 /** Converts all points in @param area to ints by casting. */
461 static public final Area areaInInts(final Area area) {
462 final Area a = new Area();
463 for (final Polygon pol : M.getPolygons(area)) {
464 a.exclusiveOr(new Area(pol));
465 }
466 return a;
467 }
468
469 /** Converts all points in @param area to ints by rounding. */
470 static public final Area areaInIntsByRounding(final Area area) {
471 final Area a = new Area();
472 for (final Polygon pol : M.getPolygonsByRounding(area)) {
473 a.add(new Area(pol));
474 }
475 return a;
476 }
477
478 /* ================================================= */
479
480 public static void quicksort(float[] data, Object[] sortAlso) throws IllegalArgumentException {
481 if (data.length != sortAlso.length) {
482 throw new IllegalArgumentException("data and sortAlso arrays don't have the same length.");
483 }
484 quicksort(data, sortAlso, 0, data.length-1);
485 }
486
487 /** Adapted from Stephan Preibisch's mpi.fruitfly.math.General homonimous method. */
488 public static void quicksort(final float[] data, final Object[] sortAlso,
489 final int left, final int right) {
490 if (data.length < 2) return;
491 int i = left, j = right;
492 float x = data[(left + right) / 2];
493 do {
494 while (data[i] < x) i++;
495 while (x < data[j]) j--;
496 if (i <= j) {
497 float temp = data[i];
498 data[i] = data[j];
499 data[j] = temp;
500
501 Object temp2 = sortAlso[i];
502 sortAlso[i] = sortAlso[j];
503 sortAlso[j] = temp2;
504
505 i++;
506 j--;
507 }
508 } while (i <= j);
509 if (left < j) quicksort(data, sortAlso, left, j);
510 if (i < right) quicksort(data, sortAlso, i, right);
511 }
512
513 static private Polygon arrowhead = null;
514
515 /** Create an arrowhead at the end of the line segment defined by x1,y1 and x2,y2. */
516 static public final Shape createArrowhead(final double x1, final double y1, final double x2, final double y2) {
517 return createArrowhead(x1, y1, x2, y2, 1.0);
518 }
519 static public final Shape createArrowhead(final double x1, final double y1, final double x2, final double y2, final double magnification) {
520 if (null == arrowhead) arrowhead = new Polygon(new int[]{-14, -13, 0, -13, -14, -9}, new int[]{-5, -5, 0, 5, 5, 0}, 6);
521 final AffineTransform aff = new AffineTransform();
522 aff.translate(x2, y2);
523 aff.rotate(M.getAngle(x2 - x1, y2 - y1));
524 if (magnification < 1.0) aff.scale(magnification, magnification);
525 return aff.createTransformedShape(arrowhead);
526 }
527
528 static final private float phi = (1 + (float)Math.sqrt(5)) / 2;
529 static final private float[][] icosahedron = { { phi, 1, 0 },
530 { -phi, 1, 0 },
531 { phi, -1, 0 },
532 { -phi, -1, 0 },
533 { 1, 0, phi },
534 { 1, 0, -phi },
535 {-1, 0, phi },
536 {-1, 0, -phi },
537 {0, phi, 1 },
538 {0, -phi, 1},
539 {0, phi, -1 },
540 {0, -phi, -1} };
541 static final private int[][] icosfaces = { { 0, 8, 4 },
542 { 0, 5, 10 },
543 { 2, 4, 9 },
544 { 2, 11, 5 },
545 { 1, 6, 8 },
546 { 1, 10, 7 },
547 { 3, 9, 6 },
548 { 3, 7, 11 },
549 { 0, 10, 8 },
550 { 1, 8, 10 },
551 { 2, 9, 11 },
552 { 3, 11, 9 },
553 { 4, 2, 0 },
554 { 5, 0, 2 },
555 { 6, 1, 3 },
556 { 7, 3, 1 },
557 { 8, 6, 4 },
558 { 9, 4, 6 },
559 { 10, 5, 7 },
560 { 11, 7, 5 } };
561
562 /** Returns a "3D Viewer"-ready list mesh, centered at 0,0,0 and with radius as the radius of the enclosing sphere. */
563 static public final List<Point3f> createIcosahedron(int subdivisions, float radius) {
564 List<Point3f> ps = new ArrayList<Point3f>();
565 for (int i=0; i<icosfaces.length; i++) {
566 for (int k=0; k<3; k++) {
567 ps.add(new Point3f(icosahedron[icosfaces[i][k]]));
568 }
569 }
570 while (subdivisions-- > 0) {
571 final List<Point3f> sub = new ArrayList<Point3f>();
572 // Take three consecutive points, which define a face, and create 4 faces out of them.
573 for (int i=0; i<ps.size(); i+=3) {
574 Point3f p0 = ps.get(i);
575 Point3f p1 = ps.get(i+1);
576 Point3f p2 = ps.get(i+2);
577
578 Point3f p01 = new Point3f((p0.x + p1.x)/2, (p0.y + p1.y)/2, (p0.z + p1.z)/2);
579 Point3f p02 = new Point3f((p0.x + p2.x)/2, (p0.y + p2.y)/2, (p0.z + p2.z)/2);
580 Point3f p12 = new Point3f((p1.x + p2.x)/2, (p1.y + p2.y)/2, (p1.z + p2.z)/2);
581 // lower left:
582 sub.add(p0);
583 sub.add(p01);
584 sub.add(p02);
585 // upper:
586 sub.add(new Point3f(p01)); // as copies
587 sub.add(p1);
588 sub.add(p12);
589 // lower right:
590 sub.add(new Point3f(p12));
591 sub.add(p2);
592 sub.add(new Point3f(p02));
593 // center:
594 sub.add(new Point3f(p01));
595 sub.add(new Point3f(p12));
596 sub.add(new Point3f(p02));
597 }
598 ps = sub;
599 }
600
601 // Project all vertices to the surface of a sphere of radius 1
602 final Vector3f v = new Vector3f();
603 for (final Point3f p : ps) {
604 v.set(p);
605 v.normalize();
606 v.scale(radius);
607 p.set(v);
608 }
609
610 return ps;
611 }
612
613 /** Reuses the @param fp to apply in place. */
614 static public final void apply(final mpicbg.models.CoordinateTransform ict, final double[][] p, final int i, final float[] fp) {
615 fp[0] = (float)p[0][i];
616 fp[1] = (float)p[1][i];
617 ict.applyInPlace(fp);
618 p[0][i] = fp[0];
619 p[1][i] = fp[1];
620 }
621
622 /** The @param ict is expected to transform the data as if this data was expressed in world coordinates,
623 * so this method returns a transformation list that prepends the transform from local to world, then the @param ict, then from world to local. */
624 static public final mpicbg.models.CoordinateTransform wrap(final AffineTransform to_world, final mpicbg.models.CoordinateTransform ict, final AffineTransform to_local) throws Exception {
625 final mpicbg.models.CoordinateTransformList<mpicbg.models.CoordinateTransform> chain
626 = new mpicbg.models.CoordinateTransformList<mpicbg.models.CoordinateTransform>(); // bravo!
627 // 1 - Prepend to world
628 final mpicbg.models.AffineModel2D toworld = new mpicbg.models.AffineModel2D();
629 toworld.set(to_world);
630 chain.add(toworld);
631 // 2 - Perform the transform in world coordinates
632 chain.add(ict);
633 // 3 - back to local
634 final mpicbg.models.AffineModel2D tolocal = new mpicbg.models.AffineModel2D();
635 tolocal.set(to_local);
636 chain.add(tolocal);
637 return chain;
638 }
639
640 /** Returns the shortest possible String representation of a float number, according to the desired decimal @param precision. */
641 static public final String shortest(final float f, final float precision) {
642 return Math.abs(f - (int)f) < precision ?
643 Integer.toString((int)f)
644 : Float.toString(f);
645 }
646
647 /** Appends to @param sb the shortest possible String representation of a float number, according to the desired decimal @param precision. */
648 static public final void appendShortest(final float f, final float precision, final StringBuilder sb) {
649 if (Math.abs(f - (int)f) < precision) {
650 sb.append((int)f);
651 } else {
652 sb.append(f);
653 }
654 }
655
656 /** @returns the point of intersection of the two segments a and b, or null if they don't intersect. */
657 public static final float[] computeSegmentsIntersection(
658 float ax0, float ay0, float ax1, float ay1,
659 float bx0, float by0, float bx1, float by1) {
660 float d = (ax0 - ax1)*(by0 - by1) - (ay0 - ay1) * (bx0 - bx1);
661 if (0 == d) return null;
662 float xi = ((bx0 - bx1)*(ax0*ay1 - ay0*ax1) - (ax0 - ax1)*(bx0*by1 - by0*bx1)) / d;
663 float yi = ((by0 - by1)*(ax0*ay1 - ay0*ax1) - (ay0 - ay1)*(bx0*by1 - by0*bx1)) / d;
664 if (xi < Math.min(ax0, ax1) || xi > Math.max(ax0, ax1)) return null;
665 if (xi < Math.min(bx0, bx1) || xi > Math.max(bx0, bx1)) return null;
666 if (yi < Math.min(ay0, ay1) || yi > Math.max(ay0, ay1)) return null;
667 if (yi < Math.min(by0, by1) || yi > Math.max(by0, by1)) return null;
668 return new float[]{xi, yi};
669 }
670
671 /** Returns an array of two Area objects, or of 1 if the @param proi doesn't intersect @param a. */
672 public static final Area[] splitArea(final Area a, final PolygonRoi proi, final Rectangle world) {
673 if (!a.getBounds().intersects(proi.getBounds())) return new Area[]{a};
674
675 int[] x = proi.getXCoordinates(),
676 y = proi.getYCoordinates();
677 final Rectangle rb = proi.getBounds();
678 final int len = proi.getNCoordinates();
679 int x0 = x[0] + rb.x;
680 int y0 = y[0] + rb.y;
681 int xN = x[len -1] + rb.x;
682 int yN = y[len -1] + rb.y;
683
684 // corners, clock-wise:
685 int[][] corners = new int[][]{new int[]{world.x, world.y}, // top left
686 new int[]{world.x + world.width, world.y}, // top right
687 new int[]{world.x + world.width, world.y + world.height}, // bottom right
688 new int[]{world.x, world.y + world.height}}; // bottom left
689 // Which corner is closest to x0,y0, and which to xN,yN
690 double min_dist0 = Double.MAX_VALUE;
691 int min_i0 = 0;
692 int i = 0;
693 double min_distN = Double.MAX_VALUE;
694 int min_iN = 0;
695 for (int[] corner : corners) {
696 double d = distance(x0, y0, corner[0], corner[1]);
697 if (d < min_dist0) {
698 min_dist0 = d;
699 min_i0 = i;
700 }
701 d = distance(xN, yN, corner[0], corner[1]);
702 if (d < min_distN) {
703 min_distN = d;
704 min_iN = i;
705 }
706 i++;
707 }
708 // Create new Area 'b' with which to intersect Area 'a':
709 int[] xb, yb;
710 int l,
711 i1,
712 i2 = -1;
713 // 0 1 2 3: when there difference is 2, there is a point in between
714 if (2 != Math.abs(min_iN - min_i0)) {
715 l = len +4;
716 xb = new int[l];
717 yb = new int[l];
718 // -4 and -1 will be the drifted corners
719 i1 = -4;
720 xb[l-4] = corners[min_iN][0];
721 yb[l-4] = corners[min_iN][1];
722 xb[l-3] = corners[min_iN][0];
723 yb[l-3] = corners[min_iN][1];
724 xb[l-2] = corners[min_i0][0];
725 yb[l-2] = corners[min_i0][1];
726 xb[l-1] = corners[min_i0][0];
727 yb[l-1] = corners[min_i0][1];
728 } else {
729 l = len +5;
730 xb = new int[l];
731 yb = new int[l];
732 // -5 and -1 will be the drifted corners
733 i1 = -5;
734 xb[l-5] = corners[min_iN][0];
735 yb[l-5] = corners[min_iN][1];
736 xb[l-4] = corners[min_iN][0];
737 yb[l-4] = corners[min_iN][1];
738 xb[l-3] = corners[(min_iN + min_i0) / 2][0];
739 yb[l-3] = corners[(min_iN + min_i0) / 2][1];
740 xb[l-2] = corners[min_i0][0];
741 yb[l-2] = corners[min_i0][1];
742 xb[l-1] = corners[min_i0][0];
743 yb[l-1] = corners[min_i0][1];
744 }
745 // Enter polyline points, corrected for proi bounds
746 for (int k=0; k<len; k++) {
747 xb[k] = x[k] + rb.x;
748 yb[k] = y[k] + rb.y;
749 }
750 // Drift corners to closest point on the sides
751 // Last:
752 int dx = Math.abs(xb[l+i1] - (x[len-1] + rb.x)),
753 dy = Math.abs(yb[l+i1] - (y[len-1] + rb.y));
754 if (dy < dx) xb[l+i1] = x[len-1] + rb.x;
755 else yb[l+i1] = y[len-1] + rb.y;
756 // First:
757 dx = Math.abs(xb[l+i2] - (x[0] + rb.x));
758 dy = Math.abs(yb[l+i2] - (y[0] + rb.y));
759 if (dy < dx) xb[l+i2] = x[0] + rb.x;
760 else yb[l+i2] = y[0] + rb.y;
761
762 Area b = new Area(new Polygon(xb, yb, xb.length));
763 Area c = new Area(world);
764 c.subtract(b);
765
766 return new Area[]{b, c};
767 }
768
769 public static final VectorString2D asVectorString2D(final Polygon pol, final double z) throws Exception {
770 double[] x = new double[pol.npoints];
771 double[] y = new double[pol.npoints];
772 for (int i=0; i<x.length; i++) {
773 x[i] = pol.xpoints[i];
774 y[i] = pol.ypoints[i];
775 }
776 return new VectorString2D(x, y, z, true);
777 }
778
779 public static final double volumeOfTruncatedCone(final double r1, final double r2, final double height) {
780 return Math.PI
781 * ( r1 * r1
782 + r1 * r2
783 + r2 * r2)
784 * height
785 / 3;
786 }
787
788 public static final double lateralAreaOfTruncatedCone(final double r1, final double r2, final double height) {
789 return Math.PI * (r1 + r2) * Math.sqrt(Math.pow(r2 - r1, 2) + height * height);
790 }
791
792 /** The lateral area plus the two circles. */
793 public static final double totalAreaOfTruncatedCone(final double r1, final double r2, final double height) {
794 return lateralAreaOfTruncatedCone(r1, r2, height)
795 + Math.PI * (r1 * r1 + r2 * r2);
796 }
797
798 public static final void convolveGaussianSigma1(final float[] in, final float[] out, final CircularSequence seq) {
799 // Weights for sigma = 1 and kernel 5x1, normalized so they add up to 1.
800 final float w2 = 0.05448869f,
801 w1 = 0.24420134f,
802 w0 = 0.40261994f;
803 // Handle boundary case
804 if (out.length < 5) {
805 for (int i=0; i<out.length; ++i) {
806 out[i] = in[seq.setPosition(i -2)] * w2
807 + in[seq.setPosition(i -1)] * w1
808 + in[seq.setPosition(i )] * w0
809 + in[seq.setPosition(i +1)] * w1
810 + in[seq.setPosition(i +2)] * w2;
811 }
812 return;
813 }
814 // Optimized general case: use the CircularSequence only at the ends
815 for (int i=0; i<2; ++i) {
816 out[i] = in[seq.setPosition(i -2)] * w2
817 + in[seq.setPosition(i -1)] * w1
818 + in[seq.setPosition(i )] * w0
819 + in[seq.setPosition(i +1)] * w1
820 + in[seq.setPosition(i +2)] * w2;
821 }
822 final int cut = out.length -2;
823 for (int i=2; i<cut; ++i) {
824 out[i] = in[i -2] * w2
825 + in[i -1] * w1
826 + in[i ] * w0
827 + in[i +1] * w1
828 + in[i +2] * w2;
829 }
830 for (int i=cut; i<out.length; ++i) {
831 out[i] = in[seq.setPosition(i -2)] * w2
832 + in[seq.setPosition(i -1)] * w1
833 + in[seq.setPosition(i )] * w0
834 + in[seq.setPosition(i +1)] * w1
835 + in[seq.setPosition(i +2)] * w2;
836 }
837 }
838
839 /** Copied from ImageJ's ij.gui.PolygonRoi.getInterpolatedPolygon, by Wayne Rasband and collaborators.
840 * The reason I copied this method is that creating a new PolygonRoi just to get an interpolated polygon
841 * processes the float[] arrays of the coordinates, subtracting the minimum x,y. Not only it is an extra
842 * operation but it is also in place, altering data arrays. Fortunately FloatPolygon does not touch the arrays. */
843 final static public FloatPolygon createInterpolatedPolygon(
844 final FloatPolygon p,
845 final double interval,
846 final boolean isLine) {
847 double length = p.getLength(isLine);
848 int npoints2 = (int)((length*1.2)/interval);
849 float[] xpoints2 = new float[npoints2];
850 float[] ypoints2 = new float[npoints2];
851 xpoints2[0] = p.xpoints[0];
852 ypoints2[0] = p.ypoints[0];
853 int n=1, n2;
854 double inc = 0.01;
855 double distance=0.0, distance2=0.0, dx=0.0, dy=0.0, xinc, yinc;
856 double x, y, lastx, lasty, x1, y1, x2=p.xpoints[0], y2=p.ypoints[0];
857 int npoints = p.npoints;
858 if (!isLine) npoints++;
859 for (int i=1; i<npoints; i++) {
860 x1=x2; y1=y2;
861 x=x1; y=y1;
862 if (i<p.npoints) {
863 x2=p.xpoints[i];
864 y2=p.ypoints[i];
865 } else {
866 x2=p.xpoints[0];
867 y2=p.ypoints[0];
868 }
869 dx = x2-x1;
870 dy = y2-y1;
871 distance = Math.sqrt(dx*dx+dy*dy);
872 xinc = dx*inc/distance;
873 yinc = dy*inc/distance;
874 lastx=xpoints2[n-1]; lasty=ypoints2[n-1];
875 //n2 = (int)(dx/xinc);
876 n2 = (int)(distance/inc);
877 if (npoints==2) n2++;
878 do {
879 dx = x-lastx;
880 dy = y-lasty;
881 distance2 = Math.sqrt(dx*dx+dy*dy);
882 //IJ.log(i+" "+IJ.d2s(xinc,5)+" "+IJ.d2s(yinc,5)+" "+IJ.d2s(distance,2)+" "+IJ.d2s(distance2,2)+" "+IJ.d2s(x,2)+" "+IJ.d2s(y,2)+" "+IJ.d2s(lastx,2)+" "+IJ.d2s(lasty,2)+" "+n+" "+n2);
883 if (distance2>=interval-inc/2.0 && n<xpoints2.length-1) {
884 xpoints2[n] = (float)x;
885 ypoints2[n] = (float)y;
886 //IJ.log("--- "+IJ.d2s(x,2)+" "+IJ.d2s(y,2)+" "+n);
887 n++;
888 lastx=x; lasty=y;
889 }
890 x += xinc;
891 y += yinc;
892 } while (--n2>0);
893 }
894 return new FloatPolygon(xpoints2, ypoints2, n);
895 }
896 }
register/edit/analyze/3D model TB of images