Skip to content

Commit 6a60e5a

Browse files
gselzerwiedenm
gselzer
authored and
wiedenm
committed
WIP: Port segment namespace
TODO: ensure passing tests
1 parent b9c81cb commit 6a60e5a

6 files changed

Lines changed: 1407 additions & 0 deletions

File tree

‎src/main/java/net/imagej/ops/segment/detectJunctions/DefaultDetectJunctions.java‎

Lines changed: 333 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,333 @@
1+
/*
2+
* #%L
3+
* ImageJ software for multidimensional image processing and analysis.
4+
* %%
5+
* Copyright (C) 2014 - 2018 ImageJ developers.
6+
* %%
7+
* Redistribution and use in source and binary forms, with or without
8+
* modification, are permitted provided that the following conditions are met:
9+
*
10+
* 1. Redistributions of source code must retain the above copyright notice,
11+
* this list of conditions and the following disclaimer.
12+
* 2. Redistributions in binary form must reproduce the above copyright notice,
13+
* this list of conditions and the following disclaimer in the documentation
14+
* and/or other materials provided with the distribution.
15+
*
16+
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
17+
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18+
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
19+
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
20+
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
21+
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
22+
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
23+
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24+
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
25+
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
26+
* POSSIBILITY OF SUCH DAMAGE.
27+
* #L%
28+
*/
29+
30+
package net.imagej.ops.segment.detectJunctions;
31+
32+
import java.util.ArrayList;
33+
import java.util.List;
34+
import java.util.function.BiFunction;
35+
import java.util.function.Function;
36+
37+
import net.imagej.ops.segment.detectRidges.DefaultDetectRidges;
38+
import net.imglib2.Interval;
39+
import net.imglib2.RealInterval;
40+
import net.imglib2.RealLocalizable;
41+
import net.imglib2.RealPoint;
42+
import net.imglib2.RealPositionable;
43+
import net.imglib2.roi.geom.real.WritablePolyline;
44+
import net.imglib2.roi.util.RealLocalizableRealPositionable;
45+
import net.imglib2.util.Intervals;
46+
47+
import org.scijava.ops.OpDependency;
48+
import org.scijava.ops.core.Op;
49+
import org.scijava.param.Parameter;
50+
import org.scijava.plugin.Plugin;
51+
import org.scijava.struct.ItemIO;
52+
53+
/**
54+
* Finds the junctions between a {@link ArrayList} of {@link WritablePolyline},
55+
* intended to be used optionally after running {@link DefaultDetectRidges} but
56+
* applicable to all groups of polylines.
57+
* <p>
58+
* TODO refactor the op to determine junction points between n-d
59+
* {@link WritablePolyline}
60+
* </p>
61+
*
62+
* @author Gabe Selzer
63+
*/
64+
@Plugin(type = Op.class, name = "segment.detectJunctions")
65+
@Parameter(key = "lines")
66+
@Parameter(key = "threshold", description = "Maximum distance between polylines to be considered a junction")
67+
@Parameter(key = "junctions", type = ItemIO.OUTPUT)
68+
public class DefaultDetectJunctions implements BiFunction<List<? extends WritablePolyline>, Double, List<RealPoint>> {
69+
70+
// @Parameter(required = false)
71+
private double threshold = 2;
72+
73+
private boolean areClose(RealPoint p1, RealPoint p2) {
74+
return getDistance(p1, p2) <= threshold;
75+
}
76+
77+
private boolean areClose(RealPoint p1, List<RealPoint> points) {
78+
for (RealPoint p : points) {
79+
if (areClose(p1, p) == true)
80+
return true;
81+
}
82+
return false;
83+
}
84+
85+
private static Interval slightlyEnlarge(RealInterval realInterval, long border) {
86+
return Intervals.expand(Intervals.smallestContainingInterval(realInterval), border);
87+
}
88+
89+
private double getDistance(double[] point1, RealLocalizable point2) {
90+
return Math.sqrt(Math.pow(point2.getDoublePosition(0) - point1[0], 2)
91+
+ Math.pow(point2.getDoublePosition(1) - point1[1], 2));
92+
}
93+
94+
private double getDistance(RealLocalizable point1, RealLocalizable point2) {
95+
return Math.sqrt(Math.pow(point2.getDoublePosition(0) - point1.getDoublePosition(0), 2)
96+
+ Math.pow(point2.getDoublePosition(1) - point1.getDoublePosition(1), 2));
97+
}
98+
99+
private RealPoint makeRealPoint(RealLocalizableRealPositionable input) {
100+
return new RealPoint(input);
101+
}
102+
103+
@Override
104+
public List<RealPoint> apply(final List<? extends WritablePolyline> input, final Double threshold) {
105+
106+
// check arguments for validity
107+
if (input.size() < 1)
108+
return new ArrayList<RealPoint>();
109+
if (input.get(0).vertex(0).numDimensions() != 2)
110+
throw new IllegalArgumentException("Only 2-dimensional WritablePolylines are supported!");
111+
112+
if (threshold != null)
113+
this.threshold = threshold;
114+
115+
// output that allows for both split polyline inputs and a
116+
// realPointCollection for our junctions.
117+
List<RealPoint> output = new ArrayList<>();
118+
119+
for (int first = 0; first < input.size() - 1; first++) {
120+
WritablePolyline firstLine = input.get(first);
121+
for (int second = first + 1; second < input.size(); second++) {
122+
WritablePolyline secondLine = input.get(second);
123+
// interval containing both plines
124+
Interval intersect = Intervals.intersect(slightlyEnlarge(firstLine, 2), slightlyEnlarge(secondLine, 2));
125+
// if the two do not intersect, then don't bother checking them against
126+
// each other.
127+
if (Intervals.isEmpty(intersect))
128+
continue;
129+
130+
// create an arraylist to contain all of the junctions for these two
131+
// lines (so that we can filter the junctions before putting them in
132+
// output).
133+
ArrayList<RealPoint> currentPairJunctions = new ArrayList<>();
134+
135+
for (int p = 0; p < firstLine.numVertices() - 1; p++) {
136+
for (int q = 0; q < secondLine.numVertices() - 1; q++) {
137+
RealLocalizableRealPositionable p1 = firstLine.vertex(p);
138+
RealLocalizableRealPositionable p2 = firstLine.vertex(p + 1);
139+
RealLocalizableRealPositionable q1 = secondLine.vertex(q);
140+
RealLocalizableRealPositionable q2 = secondLine.vertex(q + 1);
141+
142+
// special cases if both lines are vertical
143+
boolean pVertical = Math.round(p1.getDoublePosition(0)) == Math.round(p2.getDoublePosition(0));
144+
boolean qVertical = Math.round(q1.getDoublePosition(0)) == Math.round(q2.getDoublePosition(0));
145+
146+
// intersection point between the lines created by line segments p
147+
// and q.
148+
double[] intersectionPoint = new double[2];
149+
150+
// if both p and q are vertical, then p and q cannot intersect,
151+
// since they are parallel and cannot be the same.
152+
if (pVertical && qVertical) {
153+
parallelRoutine(p1, p2, q1, q2, currentPairJunctions, true);
154+
continue;
155+
} else if (pVertical) {
156+
double mq = (q2.getDoublePosition(1) - q1.getDoublePosition(1))
157+
/ (q2.getDoublePosition(0) - q1.getDoublePosition(0));
158+
double bq = (q1.getDoublePosition(1) - mq * q1.getDoublePosition(0));
159+
double x = p1.getDoublePosition(0);
160+
double y = mq * x + bq;
161+
intersectionPoint[0] = x;
162+
intersectionPoint[1] = y;
163+
} else if (qVertical) {
164+
double mp = (p2.getDoublePosition(1) - p1.getDoublePosition(1))
165+
/ (p2.getDoublePosition(0) - p1.getDoublePosition(0));
166+
double bp = (p1.getDoublePosition(1) - mp * p1.getDoublePosition(0));
167+
double x = q1.getDoublePosition(0);
168+
double y = mp * x + bp;
169+
intersectionPoint[0] = x;
170+
intersectionPoint[1] = y;
171+
} else {
172+
173+
double mp = (p2.getDoublePosition(1) - p1.getDoublePosition(1))
174+
/ (p2.getDoublePosition(0) - p1.getDoublePosition(0));
175+
double mq = (q2.getDoublePosition(1) - q1.getDoublePosition(1))
176+
/ (q2.getDoublePosition(0) - q1.getDoublePosition(0));
177+
178+
if (mp == mq) {
179+
parallelRoutine(p1, p2, q1, q2, currentPairJunctions, false);
180+
continue;
181+
}
182+
183+
double bp = (p2.getDoublePosition(1) - mp * p2.getDoublePosition(0));
184+
double bq = (q2.getDoublePosition(1) - mq * q2.getDoublePosition(0));
185+
186+
// point of intersection of lines created by line segments p and
187+
// q.
188+
double x = (bq - bp) / (mp - mq);
189+
double y = mp * x + bp;
190+
intersectionPoint[0] = x;
191+
intersectionPoint[1] = y;
192+
}
193+
194+
// find the distance from the intersection point to both line
195+
// segments, and the length of the line segments.
196+
double distp1 = getDistance(intersectionPoint, p1);
197+
double distp2 = getDistance(intersectionPoint, p2);
198+
double distq1 = getDistance(intersectionPoint, q1);
199+
double distq2 = getDistance(intersectionPoint, q2);
200+
201+
// max distance from line segment to intersection point
202+
double maxDist = Math.max(Math.min(distp1, distp2), Math.min(distq1, distq2));
203+
204+
// if the maximum distance is close enough to the two lines, then
205+
// these lines are close enough to form a junction
206+
if (maxDist <= threshold)
207+
currentPairJunctions.add(new RealPoint(intersectionPoint));
208+
}
209+
}
210+
// filter out the current pair's junctions by removing duplicates and
211+
// then averaging all remaining nearby junctions
212+
filterJunctions(currentPairJunctions);
213+
214+
// add the filtered junctions to the output list.
215+
for (RealPoint point : currentPairJunctions)
216+
output.add(point);
217+
}
218+
}
219+
220+
// filter the junctions -- for each set of junctions that seem vaguely
221+
// similar, pick out the best one-
222+
filterJunctions(output);
223+
224+
return output;
225+
}
226+
227+
private void filterJunctions(List<RealPoint> list) {
228+
// filter out all vaguely similar junction points.
229+
for (int i = 0; i < list.size() - 1; i++) {
230+
ArrayList<RealPoint> similars = new ArrayList<>();
231+
similars.add(list.get(i));
232+
list.remove(i);
233+
for (int j = 0; j < list.size(); j++) {
234+
if (areClose(list.get(j), similars)) {
235+
similars.add(list.get(j));
236+
list.remove(j);
237+
j--;
238+
}
239+
}
240+
if (list.size() > 0)
241+
list.add(i, averagePoints(similars));
242+
else
243+
list.add(averagePoints(similars));
244+
}
245+
}
246+
247+
private RealPoint averagePoints(ArrayList<RealPoint> list) {
248+
double[] pos = { 0, 0 };
249+
for (RealPoint p : list) {
250+
pos[0] += p.getDoublePosition(0);
251+
pos[1] += p.getDoublePosition(1);
252+
}
253+
pos[0] /= list.size();
254+
pos[1] /= list.size();
255+
return new RealPoint(pos);
256+
}
257+
258+
private <L extends RealLocalizable & RealPositionable> void parallelRoutine(RealLocalizableRealPositionable p1,
259+
RealLocalizableRealPositionable p2, RealLocalizableRealPositionable q1, RealLocalizableRealPositionable q2,
260+
List<RealPoint> junctions, boolean areVertical) {
261+
262+
// find out whether or not they are on the same line
263+
boolean sameLine = false;
264+
if (areVertical && Math.round(p1.getDoublePosition(0)) == Math.round(q1.getDoublePosition(0)))
265+
sameLine = true;
266+
else {
267+
double m = (q2.getDoublePosition(1) - q1.getDoublePosition(1))
268+
/ (q2.getDoublePosition(0) - q1.getDoublePosition(0));
269+
double bp = (p2.getDoublePosition(1) - m * p2.getDoublePosition(0));
270+
double bq = (q2.getDoublePosition(1) - m * q2.getDoublePosition(0));
271+
272+
if (bp == bq)
273+
sameLine = true;
274+
}
275+
276+
// if the two line segments do not belong to the same line, then if the
277+
// minimum distance between the two points is greater than the threshold,
278+
// there is no junction
279+
if (!sameLine && Math.min(Math.min(getDistance(p1, q1), getDistance(p2, q1)),
280+
Math.min(getDistance(p1, q2), getDistance(p2, q2))) > threshold)
281+
return;
282+
283+
int foundJunctions = 0;
284+
double lengthp = getDistance(p1, p2);
285+
double lengthq = getDistance(q1, q2);
286+
// if p and q are segments on the same line, then p1, p2, q1, and q2 can all
287+
// be junctions. There can be at most 2 junctions between these two line
288+
// segments.
289+
// check p1 to be a junction
290+
if ((getDistance(p1, q1) < lengthq && getDistance(p1, q2) < lengthq && sameLine)
291+
|| Math.min(getDistance(p1, q1), getDistance(p1, q2)) < threshold) {
292+
junctions.add(makeRealPoint(p1));
293+
foundJunctions++;
294+
}
295+
// check p2 to be a junction
296+
if ((getDistance(p2, q1) < lengthq && getDistance(p2, q2) < lengthq && sameLine)
297+
|| Math.min(getDistance(p2, q1), getDistance(p2, q2)) < threshold) {
298+
junctions.add(makeRealPoint(p2));
299+
foundJunctions++;
300+
}
301+
302+
// check q1 to be a junction
303+
if (((getDistance(q1, p1) < lengthp && getDistance(q1, p2) < lengthp && sameLine)
304+
|| (Math.min(getDistance(q1, p1), getDistance(q1, p2)) < threshold)) && foundJunctions < 2) {
305+
junctions.add(makeRealPoint(q1));
306+
foundJunctions++;
307+
}
308+
309+
// check q2 to be a junction
310+
if (((getDistance(q2, p1) < lengthp && getDistance(q2, p2) < lengthp && sameLine)
311+
|| (Math.min(getDistance(q2, p1), getDistance(q2, p2)) < threshold)) && foundJunctions < 2) {
312+
junctions.add(makeRealPoint(q2));
313+
foundJunctions++;
314+
}
315+
}
316+
317+
}
318+
319+
@Plugin(type = Op.class, name = "segment.detectJunctions")
320+
@Parameter(key = "lines")
321+
@Parameter(key = "threshold", description = "Maximum distance between polylines to be considered a junction")
322+
@Parameter(key = "junctions", type = ItemIO.OUTPUT)
323+
class SimpleDetectJunctions implements Function<List<? extends WritablePolyline>, List<RealPoint>> {
324+
325+
@OpDependency(name = "segment.detectJunctions")
326+
private BiFunction<List<? extends WritablePolyline>, Double, List<RealPoint>> junctionDetector;
327+
328+
@Override
329+
public List<RealPoint> apply(List<? extends WritablePolyline> t) {
330+
return junctionDetector.apply(t, null);
331+
}
332+
333+
}

0 commit comments

Comments
 (0)