diff --git a/pom.xml b/pom.xml index 9be7fb4f9e..7123055458 100644 --- a/pom.xml +++ b/pom.xml @@ -273,6 +273,11 @@ scifio test + + net.imagej + imagej-ui-swing + test + org.scijava diff --git a/src/main/java/net/imagej/ops/segment/SegmentNamespace.java b/src/main/java/net/imagej/ops/segment/SegmentNamespace.java index 4d9f4ffcea..c9e536bf16 100644 --- a/src/main/java/net/imagej/ops/segment/SegmentNamespace.java +++ b/src/main/java/net/imagej/ops/segment/SegmentNamespace.java @@ -6,13 +6,13 @@ * %% * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions are met: - * + * * 1. Redistributions of source code must retain the above copyright notice, * this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright notice, * this list of conditions and the following disclaimer in the documentation * and/or other materials provided with the distribution. - * + * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE @@ -35,16 +35,23 @@ import net.imagej.ops.Namespace; import net.imagej.ops.OpMethod; import net.imagej.ops.Ops; +import net.imagej.ops.segment.hough.HoughCircle; +import net.imglib2.IterableInterval; +import net.imglib2.RandomAccessible; import net.imglib2.RandomAccessibleInterval; import net.imglib2.RealPoint; +import net.imglib2.img.Img; import net.imglib2.roi.geom.real.WritablePolyline; +import net.imglib2.type.BooleanType; +import net.imglib2.type.NativeType; import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.real.DoubleType; import org.scijava.plugin.Plugin; /** * The segment namespace contains segmentation operations. - * + * * @author Gabe Selzer */ @Plugin(type = Namespace.class) @@ -61,17 +68,20 @@ public > List detectRidges( final int ridgeLengthMin) { @SuppressWarnings("unchecked") - final List result = (List) ops().run( - Ops.Segment.DetectRidges.class, input, width, lowerThreshold, - higherThreshold, ridgeLengthMin); + final List result = + (List) ops().run( + Ops.Segment.DetectRidges.class, input, width, lowerThreshold, + higherThreshold, ridgeLengthMin); return result; } // -- detectJunctions -- - - @OpMethod(op = net.imagej.ops.segment.detectJunctions.DefaultDetectJunctions.class) - public List detectJunctions(final List lines) + + @OpMethod( + op = net.imagej.ops.segment.detectJunctions.DefaultDetectJunctions.class) + public List detectJunctions( + final List lines) { @SuppressWarnings("unchecked") final List result = (List) ops().run( @@ -79,10 +89,11 @@ public List detectJunctions(final List li return result; } - - @OpMethod(op = net.imagej.ops.segment.detectJunctions.DefaultDetectJunctions.class) - public List detectJunctions(final List lines, - final double threshold) + + @OpMethod( + op = net.imagej.ops.segment.detectJunctions.DefaultDetectJunctions.class) + public List detectJunctions( + final List lines, final double threshold) { @SuppressWarnings("unchecked") final List result = (List) ops().run( @@ -91,6 +102,137 @@ public List detectJunctions(final List li return result; } + // -- detectHoughCircleDoG -- + + @OpMethod(op = net.imagej.ops.segment.hough.HoughCircleDetectorDogOp.class) + public & NativeType> List + detectHoughCircleDoG(final RandomAccessibleInterval in, + final double circleThickness, final double minRadius, + final double stepRadius, final double sigma) + { + @SuppressWarnings("unchecked") + final List result = (List) ops().run( + net.imagej.ops.Ops.Segment.DetectHoughCircleDoG.class, in, + circleThickness, minRadius, stepRadius, sigma); + return result; + } + + @OpMethod(op = net.imagej.ops.segment.hough.HoughCircleDetectorDogOp.class) + public & NativeType> List + detectHoughCircleDoG(final RandomAccessibleInterval in, + final double circleThickness, final double minRadius, + final double stepRadius, final double sigma, final double sensitivity) + { + @SuppressWarnings("unchecked") + final List result = (List) ops().run( + net.imagej.ops.Ops.Segment.DetectHoughCircleDoG.class, in, + circleThickness, minRadius, stepRadius, sigma, sensitivity); + return result; + } + + // -- detectHoughCircleLE -- + + @OpMethod( + op = net.imagej.ops.segment.hough.HoughCircleDetectorLocalExtremaOp.class) + public & NativeType> List + detectHoughCircleLE(final RandomAccessibleInterval in, + final double minRadius, final double stepRadius) + { + @SuppressWarnings("unchecked") + final List result = (List) ops().run( + net.imagej.ops.Ops.Segment.DetectHoughCircleLE.class, in, minRadius, + stepRadius); + return result; + } + + @OpMethod( + op = net.imagej.ops.segment.hough.HoughCircleDetectorLocalExtremaOp.class) + public & NativeType> List + detectHoughCircleLE(final RandomAccessibleInterval in, + final double minRadius, final double stepRadius, final double sensitivity) + { + @SuppressWarnings("unchecked") + final List result = (List) ops().run( + net.imagej.ops.Ops.Segment.DetectHoughCircleLE.class, in, minRadius, + stepRadius, sensitivity); + return result; + } + + // -- transformHoughCircle -- + + @OpMethod(op = net.imagej.ops.segment.hough.HoughTransformOpNoWeights.class) + public > Img transformHoughCircle( + final IterableInterval in, final long minRadius, final long maxRadius) + { + @SuppressWarnings("unchecked") + final Img result = (Img) ops().run( + net.imagej.ops.Ops.Segment.TransformHoughCircle.class, in, minRadius, + maxRadius); + return result; + } + + @OpMethod(op = net.imagej.ops.segment.hough.HoughTransformOpNoWeights.class) + public > Img transformHoughCircle( + final Img out, final IterableInterval in, + final long minRadius, final long maxRadius) + { + @SuppressWarnings("unchecked") + final Img result = (Img) ops().run( + net.imagej.ops.Ops.Segment.TransformHoughCircle.class, out, in, minRadius, + maxRadius); + return result; + } + + @OpMethod(op = net.imagej.ops.segment.hough.HoughTransformOpNoWeights.class) + public > Img transformHoughCircle( + final Img out, final IterableInterval in, + final long minRadius, final long maxRadius, final long stepRadius) + { + @SuppressWarnings("unchecked") + final Img result = (Img) ops().run( + net.imagej.ops.Ops.Segment.TransformHoughCircle.class, out, in, minRadius, + maxRadius, stepRadius); + return result; + } + + @OpMethod(op = net.imagej.ops.segment.hough.HoughTransformOpWeights.class) + public , R extends RealType> Img + transformHoughCircle(final IterableInterval in, final long minRadius, + final long maxRadius, final RandomAccessible weights) + { + @SuppressWarnings("unchecked") + final Img result = (Img) ops().run( + net.imagej.ops.Ops.Segment.TransformHoughCircle.class, in, minRadius, + maxRadius, weights); + return result; + } + + @OpMethod(op = net.imagej.ops.segment.hough.HoughTransformOpWeights.class) + public , R extends RealType> Img + transformHoughCircle(final Img out, + final IterableInterval in, final long minRadius, final long maxRadius, + final RandomAccessible weights) + { + @SuppressWarnings("unchecked") + final Img result = (Img) ops().run( + net.imagej.ops.Ops.Segment.TransformHoughCircle.class, out, in, minRadius, + maxRadius, weights); + return result; + } + + @OpMethod(op = net.imagej.ops.segment.hough.HoughTransformOpWeights.class) + public , R extends RealType> Img + transformHoughCircle(final Img out, + final IterableInterval in, final long minRadius, final long maxRadius, + final long stepRadius, final RandomAccessible weights) + { + @SuppressWarnings("unchecked") + final Img result = (Img) ops().run( + net.imagej.ops.Ops.Segment.TransformHoughCircle.class, out, in, minRadius, + maxRadius, stepRadius, weights); + return result; + } + // -- Namespace methods -- @Override diff --git a/src/main/java/net/imagej/ops/segment/hough/HoughCircle.java b/src/main/java/net/imagej/ops/segment/hough/HoughCircle.java new file mode 100644 index 0000000000..216aa91964 --- /dev/null +++ b/src/main/java/net/imagej/ops/segment/hough/HoughCircle.java @@ -0,0 +1,81 @@ +package net.imagej.ops.segment.hough; + +import net.imglib2.RealLocalizable; +import net.imglib2.roi.geom.real.AbstractWritableSphere; + +/** + * Representation of a circle segmented via Hough circle transformation. + *

+ * Instances are {@link RealLocalizable} and have a radius (in pixel units) and + * a sensitivity. The sensitivity factor is used to report a quality of the + * Hough circle detection, the smallest the better the circle. The sensitivity + * can be used to sort a list of circles. + * + * @author Jean-Yves Tinevez. + * @author Gabe Selzer. + * + */ +public class HoughCircle extends AbstractWritableSphere implements + Comparable +{ + + private final double sensitivity; + + public HoughCircle( final RealLocalizable pos, final double radius, final double sensitivity ) + { + super( getCenter(pos) , radius); + this.radius = radius; + this.sensitivity = sensitivity; + } + + private static double[] getCenter(RealLocalizable l) { + double[] center = new double[l.numDimensions()]; + l.localize(center); + return center; + } + + @Override + public String toString() + { + final StringBuilder sb = new StringBuilder(); + char c = '('; + for ( int i = 0; i < numDimensions(); i++ ) + { + sb.append( c ); + sb.append( String.format( "%.1f", center[ i ] ) ); + c = ','; + } + sb.append( ")" ); + return String.format( "%s\tR=%.1f\tSensitivity=%.1f", sb.toString(), radius, sensitivity ); + } + + public double getRadius() + { + return radius; + } + + public double getSensitivity() + { + return sensitivity; + } + + @Override + public int compareTo( final HoughCircle o ) + { + return sensitivity < o.sensitivity ? -1 : sensitivity > o.sensitivity ? +1 : 0; + } + + @Override + public boolean test( final RealLocalizable point ) + { + final double dx = center[0] - point.getDoublePosition( 0 ); + final double dy = center[1] - point.getDoublePosition( 1 ); + final double dr2 = dx * dx + dy * dy; + + if ( dr2 < radius * radius ) + return false; + + return true; + } + +} diff --git a/src/main/java/net/imagej/ops/segment/hough/HoughCircleDetectorDogOp.java b/src/main/java/net/imagej/ops/segment/hough/HoughCircleDetectorDogOp.java new file mode 100644 index 0000000000..0b22808102 --- /dev/null +++ b/src/main/java/net/imagej/ops/segment/hough/HoughCircleDetectorDogOp.java @@ -0,0 +1,113 @@ +package net.imagej.ops.segment.hough; + +import java.util.ArrayList; +import java.util.Collections; +import java.util.List; +import java.util.concurrent.ExecutorService; + +import org.scijava.plugin.Parameter; +import org.scijava.plugin.Plugin; +import org.scijava.thread.ThreadService; + +import net.imagej.ops.Ops; +import net.imagej.ops.special.function.AbstractUnaryFunctionOp; +import net.imglib2.Point; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.RealPoint; +import net.imglib2.algorithm.dog.DogDetection; +import net.imglib2.algorithm.localextrema.RefinedPeak; +import net.imglib2.type.NativeType; +import net.imglib2.type.numeric.RealType; +import net.imglib2.util.Util; +import net.imglib2.view.Views; + +/** + * A Hough circle detector based on running sub-pixel DoG detection in the vote + * image. + *

+ * This detector offers typically a more robust and accurate detection than its + * local maxima counterpart, at the cost of longer processing time. + * + * @author Jean-Yves Tinevez + * + * @param + * the type of the source vote image. + */ +@Plugin( type = Ops.Segment.DetectHoughCircleDoG.class ) +public class HoughCircleDetectorDogOp & NativeType> + extends + AbstractUnaryFunctionOp, List> + implements Ops.Segment.DetectHoughCircleDoG +{ + + private static final double K = 1.6; + + @Parameter + private ThreadService threadService; + + @Parameter( required = true, min = "1" ) + private double circleThickness; + + @Parameter( required = true, min = "1" ) + private double minRadius; + + @Parameter( required = true, min = "1" ) + private double stepRadius; + + @Parameter + private double sigma; + + /** + * Hough circle detection is a detection algorithm. Many of them returns a + * quality measure for what they detect (spot, circle, line,...) that + * reports how likely it is that it is not a spurious detection. The greater + * the quality, the more likely that it is a real one. + *

+ * The sensitivity used here which varies as the inverse of this quality. + * The sensitivity of a circle appears as the minimal value the sensitivity + * settings must be set to incorporate it in the results. + */ + @Parameter( required = false, min = "0.1" ) + private double sensitivity = 20.; + + @Override + public List< HoughCircle > calculate( final RandomAccessibleInterval< T > input ) + { + final int numDimensions = input.numDimensions(); + final ExecutorService es = threadService.getExecutorService(); + + final double threshold = 2. * Math.PI * minRadius * circleThickness / sensitivity; + final double[] calibration = Util.getArrayFromValue( 1., numDimensions ); + final DogDetection< T > dog = new DogDetection<>( + Views.extendZero( input ), + input, + calibration, + sigma / K, + sigma, + DogDetection.ExtremaType.MINIMA, + threshold, + false ); + dog.setExecutorService( es ); + final ArrayList< RefinedPeak< Point > > refined = dog.getSubpixelPeaks(); + + /* + * Create circles. + */ + + final ArrayList< HoughCircle > circles = new ArrayList<>( refined.size() ); + for ( final RefinedPeak< Point > peak : refined ) + { + // Minima are negative. + final RealPoint center = new RealPoint( numDimensions - 1 ); + for ( int d = 0; d < numDimensions - 1; d++ ) + center.setPosition( peak.getDoublePosition( d ), d ); + + final double radius = minRadius + ( peak.getDoublePosition( numDimensions - 1 ) ) * stepRadius; + final double ls = -2. * Math.PI * minRadius * circleThickness / peak.getValue(); + circles.add( new HoughCircle( center, radius, ls ) ); + } + + Collections.sort( circles ); + return circles; + } +} diff --git a/src/main/java/net/imagej/ops/segment/hough/HoughCircleDetectorLocalExtremaOp.java b/src/main/java/net/imagej/ops/segment/hough/HoughCircleDetectorLocalExtremaOp.java new file mode 100644 index 0000000000..6389b72821 --- /dev/null +++ b/src/main/java/net/imagej/ops/segment/hough/HoughCircleDetectorLocalExtremaOp.java @@ -0,0 +1,233 @@ + +package net.imagej.ops.segment.hough; + +import java.util.ArrayList; +import java.util.Collections; +import java.util.List; + +import net.imagej.ops.Ops; +import net.imagej.ops.special.function.AbstractUnaryFunctionOp; +import net.imglib2.Localizable; +import net.imglib2.Point; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.RealLocalizable; +import net.imglib2.RealPoint; +import net.imglib2.Sampler; +import net.imglib2.algorithm.localextrema.LocalExtrema; +import net.imglib2.algorithm.localextrema.LocalExtrema.LocalNeighborhoodCheck; +import net.imglib2.algorithm.localextrema.RefinedPeak; +import net.imglib2.algorithm.localextrema.SubpixelLocalization; +import net.imglib2.algorithm.neighborhood.Neighborhood; +import net.imglib2.type.NativeType; +import net.imglib2.type.numeric.RealType; +import net.imglib2.util.Util; + +import org.scijava.Cancelable; +import org.scijava.plugin.Parameter; +import org.scijava.plugin.Plugin; +import org.scijava.thread.ThreadService; + +/** + * A Hough circle detector based on finding local extrema in the vote image, + * followed by non-maxima suppression. + *

+ * This detector offers typically a faster processing than its DoG-based + * counterpart, at the cost of accuracy and robustness. + *

+ * It involves non-maxima suppression. The local extrema-based detector directly + * looks for local maxima in the vote image. The vote image might be 'noisy': + * several local maxima might correspond to the same true circle but be + * separated by a couple of pixels. This will results in having two or more + * circles detected very close to each other with a very similar radius. + * Non-maxima suppression is an ad-hoc workaround for these artifacts. You order + * all the circles by some quality factor (we use the sensitivity) and if two + * circles are found too close (in my case, one inside one another) you just + * keep the one with the best quality (lower sensitivity). This is not always + * desirable if you know that some smaller circles are supposed to be inside + * larger one. In that case you have to use the other detector. + * + * @author Jean-Yves Tinevez + * @param the type of the source vote image. + */ +@Plugin(type = Ops.Segment.DetectHoughCircleLE.class) +public class HoughCircleDetectorLocalExtremaOp & NativeType> + extends + AbstractUnaryFunctionOp, List> + implements Cancelable, Ops.Segment.DetectHoughCircleLE +{ + + @Parameter + private ThreadService threadService; + + @Parameter(min = "1") + private double minRadius; + + @Parameter(min = "1") + private double stepRadius; + + /** + * Hough circle detection is a detection algorithm. Many of them returns a + * quality measure for what they detect (spot, circle, line,...) that reports + * how likely it is that it is not a spurious detection. The greater the + * quality, the more likely that it is a real one. + *

+ * The sensitivity used here which varies as the inverse of this quality. The + * sensitivity of a circle appears as the minimal value the sensitivity + * settings must be set to incorporate it in the results. + */ + @Parameter(required = false, min = "0.1") + private double sensitivity = 20.; + + @Override + public List calculate(final RandomAccessibleInterval input) { + final int numDimensions = input.numDimensions(); + + /* + * Find local extrema. + */ + + final double threshold = 2. * Math.PI * minRadius / sensitivity; + final T t = Util.getTypeFromInterval(input).createVariable(); + t.setReal(threshold); + final LocalNeighborhoodCheck maximumCheck = + new LocalNeighborhoodCheck() + { + + @Override + public > Circle check(final C center, + final Neighborhood neighborhood) + { + final T c = center.get(); + if (t.compareTo(c) > 0) return null; + + for (final T t1 : neighborhood) + if (t1.compareTo(c) > 0) return null; + + final double val = c.getRealDouble(); + final double radius = minRadius + (center.getDoublePosition( + numDimensions - 1)) * stepRadius; + final double ls = 2. * Math.PI * radius / val; + return new Circle(center, radius, ls); + } + }; + final ArrayList peaks = LocalExtrema.findLocalExtrema(input, + maximumCheck, threadService.getExecutorService()); + + if (isCanceled()) return Collections.emptyList(); + + /* + * Non-maxima suppression. + * + * Rule: when one circle has a center inside one another, we discard the + * one with the highest sensitivity. + */ + + // Sort by ascending sensitivity. + Collections.sort(peaks); + + final List retained = new ArrayList<>(); + NEXT_CIRCLE: + for (final Circle tested : peaks) { + for (final Circle kept : retained) { + if (kept.contains(tested)) continue NEXT_CIRCLE; + } + + // Was not found in any circle, so we keep it. + retained.add(tested); + } + + if (isCanceled()) return Collections.emptyList(); + + /* + * Refine local extrema. + */ + + final SubpixelLocalization spl = new SubpixelLocalization<>( + numDimensions); + spl.setAllowMaximaTolerance(true); + spl.setMaxNumMoves(10); + final ArrayList> refined = spl.process(retained, input, + input); + + if (isCanceled()) return Collections.emptyList(); + + /* + * Create circles. + */ + + final ArrayList circles = new ArrayList<>(refined.size()); + for (final RefinedPeak peak : refined) { + final double radius = minRadius + (peak.getDoublePosition(numDimensions - + 1)) * stepRadius; + final double ls = 2. * Math.PI * radius / peak.getValue(); + if (ls < 0 || ls > sensitivity) continue; + final RealPoint center = new RealPoint(numDimensions - 1); + for (int d = 0; d < numDimensions - 1; d++) + center.setPosition(peak.getDoublePosition(d), d); + + circles.add(new HoughCircle(center, radius, ls)); + } + + Collections.sort(circles); + return circles; + } + + // -- Cancelable methods -- + + /** Reason for cancelation, or null if not canceled. */ + private String cancelReason; + + @Override + public boolean isCanceled() { + return cancelReason != null; + } + + /** Cancels the command execution, with the given reason for doing so. */ + @Override + public void cancel(final String reason) { + cancelReason = reason == null ? "" : reason; + } + + @Override + public String getCancelReason() { + return cancelReason; + } + + /* + * INNER CLASSES + */ + + /** + * A circle class defined on integer position with a sensitivity value. + * + * @author Jean-Yves Tinevez + */ + private final class Circle extends Point implements Comparable { + + private final double radius; + + private final double sensitivity1; + + public Circle(final Localizable pos, final double radius, + final double sensitivity) + { + super(pos); + this.radius = radius; + this.sensitivity1 = sensitivity; + } + + @Override + public int compareTo(final Circle o) { + return sensitivity1 < o.sensitivity1 ? -1 : sensitivity1 > o.sensitivity1 + ? +1 : 0; + } + + public boolean contains(final RealLocalizable point) { + final double dx = getDoublePosition(0) - point.getDoublePosition(0); + final double dy = getDoublePosition(1) - point.getDoublePosition(1); + final double dr2 = dx * dx + dy * dy; + return dr2 <= radius * radius; + } + + } +} diff --git a/src/main/java/net/imagej/ops/segment/hough/HoughTransformOpNoWeights.java b/src/main/java/net/imagej/ops/segment/hough/HoughTransformOpNoWeights.java new file mode 100644 index 0000000000..f28c5efceb --- /dev/null +++ b/src/main/java/net/imagej/ops/segment/hough/HoughTransformOpNoWeights.java @@ -0,0 +1,122 @@ + +package net.imagej.ops.segment.hough; + +import net.imagej.ops.Contingent; +import net.imagej.ops.Ops; +import net.imagej.ops.special.hybrid.AbstractUnaryHybridCF; +import net.imglib2.Cursor; +import net.imglib2.Dimensions; +import net.imglib2.FinalDimensions; +import net.imglib2.IterableInterval; +import net.imglib2.img.Img; +import net.imglib2.img.ImgFactory; +import net.imglib2.type.BooleanType; +import net.imglib2.type.numeric.real.DoubleType; +import net.imglib2.view.IntervalView; +import net.imglib2.view.Views; + +import org.scijava.Priority; +import org.scijava.app.StatusService; +import org.scijava.plugin.Parameter; +import org.scijava.plugin.Plugin; + +/** + * Hough transform for binary images. + * + * @author Jean-Yves Tinevez + * @param the type of source image. Must extend boolean type. + */ +@Plugin(type = Ops.Segment.TransformHoughCircle.class, priority = Priority.HIGH) +public class HoughTransformOpNoWeights> extends + AbstractUnaryHybridCF, Img> implements + Contingent, Ops.Segment.TransformHoughCircle +{ + + @Parameter + protected StatusService statusService; + + /** + * The minimum radius to search for. + */ + @Parameter(label = "Min circle radius", + description = "Minimal radius, in pixel units, for the transform.", + min = "1") + protected long minRadius = 1; + + /** + * The maximum radius to search for. + */ + @Parameter(label = "Max circle radius", + description = "Maximal radius, in pixel units, for the transform.", + min = "1") + protected long maxRadius = 50; + + /** + * The increment jumped in radius size after every pass. + */ + @Parameter(required = false, label = "Step radius", + description = "Radius step, in pixel units, for the transform.", min = "1") + protected long stepRadius = 1; + + @Override + public boolean conforms() { + return in().numDimensions() == 2; + } + + @Override + public Img createOutput(final IterableInterval input) { + final long maxR = Math.max(minRadius, maxRadius); + final long minR = Math.min(minRadius, maxRadius); + maxRadius = maxR; + minRadius = minR; + final long nRadiuses = (maxRadius - minRadius) / stepRadius + 1; + + /* + * Voting image. + */ + + final int numDimensions = input.numDimensions(); + final long[] dims = new long[numDimensions + 1]; + for (int d = 0; d < numDimensions; d++) + dims[d] = input.dimension(d); + dims[numDimensions] = nRadiuses; + final Dimensions dimensions = FinalDimensions.wrap(dims); + final ImgFactory factory = ops().create().imgFactory( + dimensions); + final Img votes = factory.create(dimensions, new DoubleType()); + return votes; + } + + @Override + public void compute(final IterableInterval input, + final Img output) + { + final long maxR = Math.max(minRadius, maxRadius); + final long minR = Math.min(minRadius, maxRadius); + maxRadius = maxR; + minRadius = minR; + final long nRadiuses = (maxRadius - minRadius) / stepRadius + 1; + + /* + * Hough transform. + */ + + final double sum = ops().stats().sum(input).getRealDouble(); + int progress = 0; + + final Cursor cursor = input.localizingCursor(); + while (cursor.hasNext()) { + cursor.fwd(); + if (!cursor.get().get()) continue; + + for (int i = 0; i < nRadiuses; i++) { + final IntervalView slice = Views.hyperSlice(output, input + .numDimensions(), i); + final long r = minRadius + i * stepRadius; + MidPointAlgorithm.inc(Views.extendZero(slice), cursor, r); + } + + statusService.showProgress(++progress, (int) sum); + } + } +} diff --git a/src/main/java/net/imagej/ops/segment/hough/HoughTransformOpWeights.java b/src/main/java/net/imagej/ops/segment/hough/HoughTransformOpWeights.java new file mode 100644 index 0000000000..c93bd62901 --- /dev/null +++ b/src/main/java/net/imagej/ops/segment/hough/HoughTransformOpWeights.java @@ -0,0 +1,71 @@ + +package net.imagej.ops.segment.hough; + +import net.imagej.ops.Ops; +import net.imglib2.Cursor; +import net.imglib2.IterableInterval; +import net.imglib2.RandomAccess; +import net.imglib2.RandomAccessible; +import net.imglib2.img.Img; +import net.imglib2.type.BooleanType; +import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.real.DoubleType; +import net.imglib2.view.IntervalView; +import net.imglib2.view.Views; + +import org.scijava.plugin.Parameter; +import org.scijava.plugin.Plugin; + +/** + * Hough transform for binary images, with weights for each pixel. + * + * @author Jean-Yves Tinevez + * @param the type of source image. Must extend boolean type. + */ +@Plugin(type = Ops.Segment.TransformHoughCircle.class) +public class HoughTransformOpWeights, R extends RealType> + extends HoughTransformOpNoWeights +{ + + @Parameter(label = "Weights", + description = "Weight image for the vote image.") + private RandomAccessible weights; + + @Override + public void compute(final IterableInterval input, + final Img output) + { + final long maxR = Math.max(minRadius, maxRadius); + final long minR = Math.min(minRadius, maxRadius); + maxRadius = maxR; + minRadius = minR; + final long nRadiuses = (maxRadius - minRadius) / stepRadius + 1; + + /* + * Hough transform. + */ + + final DoubleType weight = new DoubleType(Double.NaN); + final RandomAccess raWeight = weights.randomAccess(input); + final double sum = ops().stats().sum(input).getRealDouble(); + int progress = 0; + + final Cursor cursor = input.localizingCursor(); + while (cursor.hasNext()) { + cursor.fwd(); + if (!cursor.get().get()) continue; + + raWeight.setPosition(cursor); + weight.set(raWeight.get().getRealDouble()); + + for (int i = 0; i < nRadiuses; i++) { + final IntervalView slice = Views.hyperSlice(output, input + .numDimensions(), i); + final long r = minRadius + i * stepRadius; + MidPointAlgorithm.add(Views.extendZero(slice), cursor, r, weight); + } + + statusService.showProgress(++progress, (int) sum); + } + } +} diff --git a/src/main/java/net/imagej/ops/segment/hough/MidPointAlgorithm.java b/src/main/java/net/imagej/ops/segment/hough/MidPointAlgorithm.java new file mode 100644 index 0000000000..0bd6232811 --- /dev/null +++ b/src/main/java/net/imagej/ops/segment/hough/MidPointAlgorithm.java @@ -0,0 +1,363 @@ +package net.imagej.ops.segment.hough; + +import net.imglib2.Cursor; +import net.imglib2.Localizable; +import net.imglib2.RandomAccess; +import net.imglib2.RandomAccessible; +import net.imglib2.Sampler; +import net.imglib2.type.numeric.RealType; + +/** + * Write circles in an image using the mid-point algorithm. + * + * @author Jean-Yves Tinevez + * @see Midpoint + * circle algorithm on Wikipedia. + * + */ +public class MidPointAlgorithm +{ + + public static final < T > Cursor< T > cursor( final RandomAccessible< T > rai, final Localizable center, final long radius ) + { + return new MidPointCursor< T >( rai, center, radius, 0, 1 ); + } + + /** + * Writes a circle in the target {@link RandomAccessible}. The circle is + * written by incrementing the pixel values by 1 along the circle. + * + * @param rai + * the random accessible. It is the caller responsibility to + * ensure it can be accessed everywhere the circle will be + * iterated. + * @param center + * the circle center. Must be at least of dimension 2. Dimensions + * 0 and 1 are used to specify the circle center. + * @param radius + * the circle radius. The circle is written in a plane in + * dimensions 0 and 1. + * @param + * the type of the target image. + */ + public static < T extends RealType< T > > void inc( final RandomAccessible< T > rai, final Localizable center, final long radius ) + { + final Cursor< T > cursor = cursor( rai, center, radius ); + while ( cursor.hasNext() ) + { + cursor.fwd(); + cursor.get().inc(); + } + } + + /** + * Writes a circle in the target {@link RandomAccessible}. The circle is + * written by setting the pixel values with the specified value. + * + * @param rai + * the target random accessible. It is the caller responsibility + * to ensure it can be accessed everywhere the circle will be + * iterated. + * @param center + * the circle center. Must be at least of dimension 2. Dimensions + * 0 and 1 are used to specify the circle center. + * @param radius + * the circle radius. The circle is written in a plane in + * dimensions 0 and 1. + * @param value + * the value to write along the circle. + * @param + * the type of the target image. + */ + public static < T extends RealType< T > > void set( final RandomAccessible< T > rai, final Localizable center, final long radius, final T value ) + { + final Cursor< T > cursor = cursor( rai, center, radius ); + while ( cursor.hasNext() ) + { + cursor.fwd(); + cursor.get().set( value ); + } + } + + /** + * Writes a circle in the target {@link RandomAccessible}. The circle is + * written by adding the specified value to the pixel values already + * in the image. + * + * @param rai + * the random accessible. It is the caller responsibility to + * ensure it can be accessed everywhere the circle will be + * iterated. + * @param center + * the circle center. Must be at least of dimension 2. Dimensions + * 0 and 1 are used to specify the circle center. + * @param radius + * the circle radius. The circle is written in a plane in + * dimensions 0 and 1. + * @param value + * the value to add along the circle. + * @param + * the type of the target image. + */ + public static < T extends RealType< T > > void add( final RandomAccessible< T > rai, final Localizable center, final long radius, final T value ) + { + final Cursor< T > cursor = cursor( rai, center, radius ); + while ( cursor.hasNext() ) + { + cursor.fwd(); + cursor.get().add( value ); + } + } + + private MidPointAlgorithm() + {} + + private static final class MidPointCursor< T > implements Cursor< T > + { + + private final RandomAccessible< T > rai; + + private final RandomAccess< T > ra; + + private final Localizable center; + + private final long radius; + + private final int dimX; + + private final int dimY; + + private long x; + + private long y; + + private long dx; + + private long dy; + + private long f; + + private Octant octant; + + private boolean hasNext; + + + private static enum Octant + { + INIT, EAST, NORTH, WEST, SOUTH, O1, O2, O3, O4, O5, O6, O7, O8; + } + + public MidPointCursor( final RandomAccessible< T > rai, final Localizable center, final long radius, final int dimX, final int dimY ) + { + this.rai = rai; + this.center = center; + this.radius = radius; + this.dimX = dimX; + this.dimY = dimY; + this.ra = rai.randomAccess(); + reset(); + } + + @Override + public void reset() + { + x = 0; + y = radius; + f = 1 - radius; + dx = 1; + dy = -2 * radius; + octant = Octant.INIT; + ra.setPosition( center ); + hasNext = true; + } + + @Override + public void fwd() + { + switch ( octant ) + { + default: + case INIT: + ra.setPosition( center.getLongPosition( dimY ) + radius, dimY ); + octant = Octant.NORTH; + break; + + case NORTH: + ra.setPosition( center.getLongPosition( dimY ) - radius, dimY ); + octant = Octant.SOUTH; + break; + + case SOUTH: + ra.setPosition( center.getLongPosition( dimX ) - radius, dimX ); + ra.setPosition( center.getLongPosition( dimY ), dimY ); + octant = Octant.WEST; + break; + + case WEST: + ra.setPosition( center.getLongPosition( dimX ) + radius, dimX ); + octant = Octant.EAST; + break; + + case EAST: + x = x + 1; + dx = dx + 2; + f = f + dx; + ra.setPosition( center.getLongPosition( dimX ) + x, dimX ); + ra.setPosition( center.getLongPosition( dimY ) + y, dimY ); + octant = Octant.O1; + break; + + case O1: + ra.setPosition( center.getLongPosition( dimX ) - x, dimX ); + octant = Octant.O2; + break; + + case O2: + ra.setPosition( center.getLongPosition( dimY ) - y, dimY ); + octant = Octant.O3; + break; + + case O3: + ra.setPosition( center.getLongPosition( dimX ) + x, dimX ); + octant = Octant.O4; + // Stop here if x==y, lest the 45ยบ will be iterated twice. + if ( x >= y ) + hasNext = false; + break; + + case O4: + ra.setPosition( center.getLongPosition( dimX ) + y, dimX ); + ra.setPosition( center.getLongPosition( dimY ) - x, dimY ); + octant = Octant.O5; + break; + + case O5: + ra.setPosition( center.getLongPosition( dimX ) - y, dimX ); + octant = Octant.O6; + break; + + case O6: + ra.setPosition( center.getLongPosition( dimY ) + x, dimY ); + octant = Octant.O7; + break; + + case O7: + ra.setPosition( center.getLongPosition( dimX ) + y, dimX ); + octant = Octant.O8; + // Stop here if dx would cross y. + if ( x >= y - 1 ) + hasNext = false; + break; + + case O8: + if ( f > 0 ) + { + y = y - 1; + dy = dy + 2; + f = f + dy; + } + x = x + 1; + dx = dx + 2; + f = f + dx; + ra.setPosition( center.getLongPosition( dimX ) + x, dimX ); + ra.setPosition( center.getLongPosition( dimY ) + y, dimY ); + octant = Octant.O1; + break; + } + } + + @Override + public boolean hasNext() + { + return hasNext; + } + + @Override + public void localize( final float[] position ) + { + ra.localize( position ); + } + + @Override + public void localize( final double[] position ) + { + ra.localize( position ); + } + + @Override + public float getFloatPosition( final int d ) + { + return ra.getFloatPosition( d ); + } + + @Override + public double getDoublePosition( final int d ) + { + return ra.getDoublePosition( d ); + } + + @Override + public int numDimensions() + { + return ra.numDimensions(); + } + + @Override + public T get() + { + return ra.get(); + } + + @Override + public Sampler< T > copy() + { + return ra.copy(); + } + + @Override + public void jumpFwd( final long steps ) + { + for ( int i = 0; i < steps; i++ ) + fwd(); + } + + @Override + public T next() + { + fwd(); + return get(); + } + + @Override + public void localize( final int[] position ) + { + ra.localize( position ); + } + + @Override + public void localize( final long[] position ) + { + ra.localize( position ); + } + + @Override + public int getIntPosition( final int d ) + { + return ra.getIntPosition( d ); + } + + @Override + public long getLongPosition( final int d ) + { + return ra.getLongPosition( d ); + } + + @Override + public Cursor< T > copyCursor() + { + return new MidPointCursor<>( rai, center, radius, dimX, dimY ); + } + + } +} diff --git a/src/main/templates/net/imagej/ops/Ops.list b/src/main/templates/net/imagej/ops/Ops.list index e38e8c72f6..f3f18aacf0 100644 --- a/src/main/templates/net/imagej/ops/Ops.list +++ b/src/main/templates/net/imagej/ops/Ops.list @@ -341,6 +341,9 @@ namespaces = ``` [name: "segment", iface: "Segment", ops: [ [name: "detectRidges", iface: "DetectRidges"], [name: "detectJunctions", iface: "DetectJunctions"], + [name: "detectHoughCircleDoG", iface: "DetectHoughCircleDoG"], + [name: "detectHoughCircleLE", iface: "DetectHoughCircleLE"], + [name: "transformHoughCircle", iface: "TransformHoughCircle"], ]], [name: "stats", iface: "Stats", ops: [ [name: "geometricMean", iface: "GeometricMean"], diff --git a/src/test/java/net/imagej/ops/segment/hough/HoughDemo.java b/src/test/java/net/imagej/ops/segment/hough/HoughDemo.java new file mode 100644 index 0000000000..d8977dcba7 --- /dev/null +++ b/src/test/java/net/imagej/ops/segment/hough/HoughDemo.java @@ -0,0 +1,75 @@ + +package net.imagej.ops.segment.hough; + +import java.util.List; + +import net.imagej.ops.OpService; +import net.imagej.ops.Ops; +import net.imagej.ops.special.function.Functions; +import net.imagej.ops.special.hybrid.Hybrids; +import net.imglib2.img.Img; +import net.imglib2.type.logic.BitType; +import net.imglib2.type.numeric.integer.UnsignedByteType; +import net.imglib2.type.numeric.real.DoubleType; + +import org.scijava.Context; +import org.scijava.app.AppService; +import org.scijava.io.IOService; +import org.scijava.ui.UIService; + +public class HoughDemo { + + public static void main(final String... args) throws Exception { + final Context context = new Context(); + final UIService ui = context.service(UIService.class); + final IOService io = context.service(IOService.class); + final OpService ops = context.service(OpService.class); + final AppService app = context.service(AppService.class); + + ui.showUI(); + + final String path = app.getApp().getBaseDirectory() + + // TEMP: The need for "/../.." is a bug, which is + // fixed on SJC master now, but for now we need it. + "/../../src/test/resources/net/imagej/ops/segment/hough/circles.png"; + System.out.println(path); + final Img input = (Img) io.open(path); + final Img img = (Img) ops.run(Ops.Convert.Bit.class, + input); + + /* + * Hough transform. + */ + final long minCircleRadius = 40; + final long maxCircleRadius = 250; + final long sigma = 1; + final long stepRadius = 3; + @SuppressWarnings({ "unchecked", "rawtypes" }) + final HoughTransformOpNoWeights hough = + (HoughTransformOpNoWeights) Hybrids.unaryCF(ops, + HoughTransformOpNoWeights.class, Img.class, img, minCircleRadius, + maxCircleRadius, stepRadius); + + final Img voteImg = hough.calculate(img); + + ui.show("Two circles", img); + ui.show("Vote", voteImg); + + /* + * Hough detection. + */ + + final double sensitivity = 1.; + @SuppressWarnings({ "rawtypes", "unchecked" }) + final HoughCircleDetectorDogOp houghDetector = + (HoughCircleDetectorDogOp) Functions.unary(ops, + HoughCircleDetectorDogOp.class, List.class, voteImg, minCircleRadius, + stepRadius, sigma, sensitivity); + final List circles = houghDetector.calculate(voteImg); + System.out.println("Found " + circles.size() + " circles:"); + for (final HoughCircle circle : circles) { + System.out.println(circle); + } + } + +} diff --git a/src/test/resources/net/imagej/ops/segment/hough/circles.png b/src/test/resources/net/imagej/ops/segment/hough/circles.png new file mode 100644 index 0000000000..fbac4fe3b0 Binary files /dev/null and b/src/test/resources/net/imagej/ops/segment/hough/circles.png differ