Skip to content

Commit 74f2b15

Browse files
gselzerctrueden
gselzer
authored and
ctrueden
committed
WIP: port topology namespace
TODO: ensure tests pass
1 parent 4f20ae9 commit 74f2b15

11 files changed

Lines changed: 2005 additions & 0 deletions

File tree

‎src/main/java/net/imagej/ops/topology/BoxCount.java‎

Lines changed: 299 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,299 @@
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.topology;
31+
32+
import java.util.ArrayList;
33+
import java.util.Arrays;
34+
import java.util.List;
35+
import java.util.Spliterator;
36+
import java.util.stream.IntStream;
37+
import java.util.stream.LongStream;
38+
import java.util.stream.Stream;
39+
import java.util.stream.StreamSupport;
40+
41+
import net.imglib2.RandomAccessibleInterval;
42+
import net.imglib2.roi.labeling.BoundingBox;
43+
import net.imglib2.type.BooleanType;
44+
import net.imglib2.type.numeric.integer.LongType;
45+
import net.imglib2.type.numeric.real.DoubleType;
46+
import net.imglib2.util.ValuePair;
47+
import net.imglib2.view.IntervalView;
48+
import net.imglib2.view.Views;
49+
50+
import org.scijava.ops.core.Op;
51+
import org.scijava.ops.core.function.Function5;
52+
import org.scijava.param.Parameter;
53+
import org.scijava.plugin.Plugin;
54+
import org.scijava.struct.ItemIO;
55+
56+
/**
57+
* An N-dimensional box counting that can be used to estimate the fractal
58+
* dimension of an interval
59+
* <p>
60+
* The algorithm repeatedly lays a fixed grid on the interval, and counts the
61+
* number of sections that contain foreground. After each step the grid is made
62+
* finer by a factor of {@link #scaling}. If the objects in the interval are
63+
* fractal, the proportion of foreground sections should increase as the grid
64+
* gets finer.
65+
* </p>
66+
* <p>
67+
* Produces a set of points (log(foreground count), -log(section size)) for
68+
* curve fitting. The slope of the function gives the fractal dimension of the
69+
* interval.
70+
* </p>
71+
*
72+
* @author Richard Domander (Royal Veterinary College, London)
73+
* @author Per Christian Henden
74+
* @author Jens Bache-Wiig
75+
*/
76+
public class BoxCount {
77+
78+
// /** Starting size of the grid sections in pixels */
79+
// private Long maxSize = 48L;
80+
//
81+
// /** Minimum size of the grid sections in pixels */
82+
// private Long minSize = 6L;
83+
//
84+
// /** Grid downscaling factor */
85+
// private Double scaling = 1.2;
86+
//
87+
// /**
88+
// * Number of times the grid is moved in each dimension to find the best fit
89+
// * <p>
90+
// * The best fitting grid covers the objects in the interval with the least
91+
// * amount of sections.
92+
// * </p>
93+
// * <p>
94+
// * NB Additional moves multiply algorithm's time complexity by n^d!
95+
// * </p>
96+
// */
97+
// private Long gridMoves = 0L;
98+
99+
/**
100+
* Counts the number of foreground sections in the interval repeatedly with
101+
* different size sections
102+
*
103+
* @param input
104+
* an n-dimensional binary interval
105+
* @return A list of (log(foreground count), -log(section size))
106+
* {@link ValuePair} objects for curve fitting
107+
*/
108+
public static <B extends BooleanType<B>> List<ValuePair<DoubleType, DoubleType>> apply(
109+
final RandomAccessibleInterval<B> input, final Long maxSize, final Long minSize, final Double scaling,
110+
final Long gridMoves) {
111+
final List<ValuePair<DoubleType, DoubleType>> points = new ArrayList<>();
112+
final int dimensions = input.numDimensions();
113+
final long[] sizes = new long[dimensions];
114+
final long numTranslations = 1 + gridMoves;
115+
input.dimensions(sizes);
116+
for (long sectionSize = maxSize; sectionSize >= minSize; sectionSize /= scaling) {
117+
final long translationAmount = Math.max(1, sectionSize / numTranslations);
118+
final Stream<long[]> translations = translationStream(numTranslations, translationAmount, dimensions - 1,
119+
new long[dimensions]);
120+
final LongStream foregroundCounts = countTranslatedGrids(input, translations, sizes, sectionSize);
121+
final long foreground = foregroundCounts.min().orElse(0);
122+
final double logSize = -Math.log(sectionSize);
123+
final double logCount = Math.log(foreground);
124+
final ValuePair<DoubleType, DoubleType> point = new ValuePair<>(new DoubleType(logSize),
125+
new DoubleType(logCount));
126+
points.add(point);
127+
}
128+
return points;
129+
}
130+
131+
/**
132+
* Count foreground sections in all grids created from the translations
133+
*
134+
* @param input
135+
* N-dimensional binary interval
136+
* @param translations
137+
* Stream of translation coordinates in n-dimensions
138+
* @param sizes
139+
* Sizes of the interval's dimensions in pixels
140+
* @param sectionSize
141+
* Size of a section in the grids
142+
* @return Foreground sections counted in each grid
143+
*/
144+
private static <B extends BooleanType<B>> LongStream countTranslatedGrids(final RandomAccessibleInterval<B> input,
145+
final Stream<long[]> translations, final long[] sizes, final long sectionSize) {
146+
final int lastDimension = sizes.length - 1;
147+
final LongType foreground = new LongType();
148+
final long[] sectionPosition = new long[sizes.length];
149+
return translations.mapToLong(gridOffset -> {
150+
foreground.setZero();
151+
Arrays.fill(sectionPosition, 0);
152+
countGrid(input, lastDimension, sizes, gridOffset, sectionPosition, sectionSize, foreground);
153+
return foreground.get();
154+
});
155+
}
156+
157+
/**
158+
* Creates a {@link net.imglib2.View} of the given grid section in the interval
159+
* <p>
160+
* Fits the view inside the bounds of the interval.
161+
* </p>
162+
*
163+
* @param interval
164+
* An n-dimensional interval with binary elements
165+
* @param sizes
166+
* Sizes of the interval's dimensions
167+
* @param coordinates
168+
* Starting coordinates of the section
169+
* @param sectionSize
170+
* Size of the section (n * n * ... n)
171+
* @return A view of the interval spanning n pixels in each dimension from the
172+
* coordinates. Null if view couldn't be set inside the interval
173+
*/
174+
private static <B extends BooleanType<B>> IntervalView<B> sectionView(final RandomAccessibleInterval<B> interval,
175+
final long[] sizes, final long[] coordinates, final long sectionSize) {
176+
final int n = sizes.length;
177+
final long[] startPosition = IntStream.range(0, n).mapToLong(i -> Math.max(0, coordinates[i])).toArray();
178+
final long[] endPosition = IntStream.range(0, n)
179+
.mapToLong(i -> Math.min((sizes[i] - 1), (coordinates[i] + sectionSize - 1))).toArray();
180+
final boolean badBox = IntStream.range(0, n).anyMatch(
181+
d -> (startPosition[d] >= sizes[d]) || (endPosition[d] < 0) || (endPosition[d] < startPosition[d]));
182+
if (badBox) {
183+
return null;
184+
}
185+
final BoundingBox box = new BoundingBox(n);
186+
box.update(startPosition);
187+
box.update(endPosition);
188+
return Views.offsetInterval(interval, box);
189+
}
190+
191+
/** Checks if the view has any foreground elements */
192+
private static <B extends BooleanType<B>> boolean hasForeground(IntervalView<B> view) {
193+
final Spliterator<B> spliterator = view.spliterator();
194+
return StreamSupport.stream(spliterator, false).anyMatch(BooleanType::get);
195+
}
196+
197+
/**
198+
* Recursively counts the number of foreground sections in the grid over the
199+
* given interval
200+
*
201+
* @param interval
202+
* An n-dimensional interval with binary elements
203+
* @param dimension
204+
* Current dimension processed, start from the last
205+
* @param sizes
206+
* Sizes of the interval's dimensions in pixels
207+
* @param translation
208+
* Translation of grid start in each dimension
209+
* @param sectionPosition
210+
* The accumulated position of the current grid section (start from
211+
* [0, 0, ... 0])
212+
* @param sectionSize
213+
* Size of a grid section (n * n * ... n)
214+
* @param foreground
215+
* Number of foreground sections found so far (start from 0)
216+
*/
217+
private static <B extends BooleanType<B>> void countGrid(final RandomAccessibleInterval<B> interval,
218+
final int dimension, final long[] sizes, final long[] translation, final long[] sectionPosition,
219+
final long sectionSize, final LongType foreground) {
220+
for (int p = 0; p < sizes[dimension]; p += sectionSize) {
221+
sectionPosition[dimension] = translation[dimension] + p;
222+
if (dimension == 0) {
223+
final IntervalView<B> box = sectionView(interval, sizes, sectionPosition, sectionSize);
224+
if (box != null && hasForeground(box)) {
225+
foreground.inc();
226+
}
227+
} else {
228+
countGrid(interval, dimension - 1, sizes, translation, sectionPosition, sectionSize, foreground);
229+
}
230+
}
231+
}
232+
233+
/**
234+
* Creates a {@link Stream} of t * 2^n translations in n-dimensions
235+
* <p>
236+
* The elements in the stream are arrays of coordinates [0, 0, .. 0], [-i, 0, ..
237+
* 0], [0, -i, 0, .. 0], .. [-i, -i, .. -i], [-2i, 0, .. 0] .. [-ti, -ti, ..
238+
* -ti], where each array has n elements, i = number of pixels translated, and t
239+
* = number of translations. If the translations were positive, a part of the
240+
* interval would not get inspected, because it always starts from [0, 0, ...
241+
* 0].
242+
* </p>
243+
* <p>
244+
* The order of arrays in the stream is not guaranteed.
245+
* </p>
246+
*
247+
* @param numTranslations
248+
* Number of translations (1 produces Stream.of(long[]{0, 0, .. 0}))
249+
* @param amount
250+
* Number of pixels shifted in translations
251+
* @param dimension
252+
* Current translation dimension (start from last)
253+
* @param translation
254+
* The accumulated position of the current translation (start from
255+
* {0, 0, .. 0})
256+
* @return A stream of coordinates of the translations
257+
*/
258+
private static Stream<long[]> translationStream(final long numTranslations, final long amount, final int dimension,
259+
final long[] translation) {
260+
final Stream.Builder<long[]> builder = Stream.builder();
261+
generateTranslations(numTranslations, amount, dimension, translation, builder);
262+
return builder.build();
263+
}
264+
265+
/**
266+
* Adds translations to the given {@link Stream.Builder}
267+
*
268+
* @see #translationStream(long, long, int, long[])
269+
*/
270+
private static void generateTranslations(final long numTranslations, final long amount, final int dimension,
271+
final long[] translation, final Stream.Builder<long[]> builder) {
272+
for (int t = 0; t < numTranslations; t++) {
273+
translation[dimension] = -t * amount;
274+
if (dimension == 0) {
275+
builder.add(translation.clone());
276+
} else {
277+
generateTranslations(numTranslations, amount, dimension - 1, translation, builder);
278+
}
279+
}
280+
}
281+
}
282+
283+
@Plugin(type = Op.class, name = "topology.boxCount")
284+
@Parameter(key = "input")
285+
@Parameter(key = "maxSize")
286+
@Parameter(key = "minSize")
287+
@Parameter(key = "scaling")
288+
@Parameter(key = "gridMoves")
289+
@Parameter(key = "output", type = ItemIO.OUTPUT)
290+
class DefaultBoxCount<B extends BooleanType<B>> implements
291+
Function5<RandomAccessibleInterval<B>, Long, Long, Double, Long, List<ValuePair<DoubleType, DoubleType>>> {
292+
293+
@Override
294+
public List<ValuePair<DoubleType, DoubleType>> apply(RandomAccessibleInterval<B> input, Long maxSize, Long minSize,
295+
Double scaling, Long gridMoves) {
296+
return BoxCount.apply(input, maxSize, minSize, scaling, gridMoves);
297+
}
298+
299+
}

0 commit comments

Comments
 (0)