From 934f45dfcef895e4e2a4b4a80e3cdf973a60ca98 Mon Sep 17 00:00:00 2001 From: Jean-Yves Tinevez Date: Mon, 3 Apr 2017 17:21:36 +0200 Subject: [PATCH] Add Hough circle transform and detector ops The transform is: - 2D version only. - weigthed or not. The detectors retrieve the Hough circles from a vote image generated by a HoughCircleTransformOp. The circles are returned as a collection of HoughCircles. Two versions: - one with a DoG detector and sub-pixel accuracy. - one with plain local extrema. Signed-off-by: Curtis Rueden --- pom.xml | 5 + .../imagej/ops/segment/SegmentNamespace.java | 168 +++++++- .../imagej/ops/segment/hough/HoughCircle.java | 81 ++++ .../hough/HoughCircleDetectorDogOp.java | 113 ++++++ .../HoughCircleDetectorLocalExtremaOp.java | 233 +++++++++++ .../hough/HoughTransformOpNoWeights.java | 122 ++++++ .../hough/HoughTransformOpWeights.java | 71 ++++ .../ops/segment/hough/MidPointAlgorithm.java | 363 ++++++++++++++++++ src/main/templates/net/imagej/ops/Ops.list | 3 + .../imagej/ops/segment/hough/HoughDemo.java | 75 ++++ .../net/imagej/ops/segment/hough/circles.png | Bin 0 -> 16578 bytes 11 files changed, 1221 insertions(+), 13 deletions(-) create mode 100644 src/main/java/net/imagej/ops/segment/hough/HoughCircle.java create mode 100644 src/main/java/net/imagej/ops/segment/hough/HoughCircleDetectorDogOp.java create mode 100644 src/main/java/net/imagej/ops/segment/hough/HoughCircleDetectorLocalExtremaOp.java create mode 100644 src/main/java/net/imagej/ops/segment/hough/HoughTransformOpNoWeights.java create mode 100644 src/main/java/net/imagej/ops/segment/hough/HoughTransformOpWeights.java create mode 100644 src/main/java/net/imagej/ops/segment/hough/MidPointAlgorithm.java create mode 100644 src/test/java/net/imagej/ops/segment/hough/HoughDemo.java create mode 100644 src/test/resources/net/imagej/ops/segment/hough/circles.png 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 0000000000000000000000000000000000000000..fbac4fe3b03d48c6408408ad9f352e6be4e61f1a GIT binary patch literal 16578 zcmajHWmr_*7dU!m7+~l_2nfQQK~fQ#Pn`G^@49h>bF10CeB~BdGwwz3z~GJd*rYoW21IWZ7C`h^xv zKx4g!cL`{9bQlXr0@Y``)>My9siV`8Q?y$D543kK^WZ!T{k2bo(AWfC4E^3b8cQ}y z#BgBISQHL|QPbC^nF;VDA$SlLMC1t{>BZK%M;t5Jgxk8QD7vd_87Dxd7+-4okT)u| zpT?@*9QDTG3eCF3Mn0(C{Hs&yjPE)Ti%y4pVlMJEG)?31=&vkjq#7QrE=qdlU}m&} zPzTi)pPY(aeRM3Vv-cyCInA*m9F@#pUMGt@XI4j7Y@iui-F679B+%1+Z2Rd}b+-mO zJ^#>Te>dS+*7Vo&27@i3Msughlnx}|wc6o`hEeAjUtA)RyY*a#6ZM7amLfQJD179j_AJ=S02xVULt!u9-*nIq zD=@f*CC`B7qRU5Y00AwvBlSkx$7Y=D$%-YLQc-vV){jJ#7_cb(>n6AMXm0h}@Q;ya z&#NB((JQ)+4zxg@xG#i!Op3iVovMRQc>q$Qe7?kI%`H$}PR8u)Q zMF9Q3-wO7(;|nxU>E*Dtnw-xuR?DYvt82kzz<5i-&G~=cT(FAGq}2)GYdUBo{y(e) zjPH&T0R%|@P4(}&n*;)_{+lC?`W$`{Ay9d&;AoR>?It^l2*C=9i{k984#HthXyUI~ z&i-t6BGh2^{7Z(sUD*x_LKpT5eaw0C+)W+k;&^Au~017+!=o5((ItQ>GJODkfn`_m3MjIt1Le=ej*!@R2*I>-hYo+dxCy83fiN9N>JeAXiA>MKN%SRpJ8XU#}lr%Sk( z(X&6Xi|w%`iY_qF+kG>(==vk9<-zP~)@*P z*;P(LpjsmAT>&4K*R^1__pmNR5E% zzj!+yJ9L`ky5&Kq14#@IzRW^1_TLG(^zbPe^YHV19_ah2gfxNk4L{6QH&s+A|OMY78=p#Zqu1mr8|nNk#=B zEZ-+SX9I*gR2pYaPkQ{kNEdHU)z-w4rg)QtK|dV)B{a2tLmyj7#V{l)jM5>UC5tnvtP7kTYYHKKVIh|#VjDsGil>gnzGa1~sUe3- z=>rOPK_s>Zv=2;lyBB~d1uR!(n=3(B$Prev?6I<8X!~9Qv`m$5i#Dwhm~2yPnzJtS ztVXsC2}ukOQsF=)(}xN8$b++5?X~JcJ_zX}`{Dc7#J-*b-^Wm1HT^JR6Kqk)Qn8yn z;8flY1JV5+IwEv+3R#=BH=8unmL29@x3R0NLv&_a>9LZ^LVAVCcDD{bbZ-wRlJIPg zC^;_E8cfF0sy0FGD|3>gIFg#-{!e zu(Rb(eQxgDWW&hM^MWh4jN1h?Gi$1H3(VuC#@C%gjbD{yl$&Ihid2Ha45uOv0`67$ z*4|bLMo4>CM9~f4GTSDm3*JTC*4^c``&sHS`b_^U?q>OhfetW*5n=Oq_0N7-DI9c=16jJ;pgwBovYkbZACyzHAD^s_5R3e6U zM6E`@o}S@Ks6E*pj(Nh9w6#B#hq6z^?zF`|7MR@dJ_`4WxuZrmf5cC8tsEUPzb2tQ z_p6g{`=?5I#!KT}BHF*WXRy?1i{d+w!Y<`BGV+ZNUyBeZjqaA~DW>T#lztJxvKI3I zJTb0wX^Pb$;EoS@aYqIynjq=WrY*?p*j`Wwi_UMkwzHAx@@o%z@OJ+WfnR4RXQ*wwA%DxXpC`$dSu9QCWl(c z#$&5m`oN8?o_z^lx^XtBiPo}Z>B;u6*eb7W%D>z#=Iv{--tyPDa3ljN%9G~_jbB4l ziJx$H?jAu_o*_X&CF90(T}eRjA)5pfe$%GAK{hc+%Dih3Tlc}}kd#-)5*!*aho|?x zmtg-VFx6Zba9obR8e+-<5mL*Rxo)wB4DKD6u-QNjzAa4nTM|yY0yq&IQY|^cYH~AI zBp*}tWRwrqOpz#0r#;a~I(Ii0tnu`dhCi#ei z_ZV?6X;#*%Qfbz(#Y-2E@F`vjM46IAmu1$LCuVK#_;&`$p5}Xjgc=y}%VW!NeWN-f z*Hf{PdaOviPfSvYxyp%HUl-HuKAQSi8JYCG)0`Sq35qfLJh_6<@gO9z)LUFI(SW9(eLA`5;rT#p|(+voQ{PI)|?cUXxQ3zF2yZL_9 z%}QE@4f8tx1^TIzlJ|EM9zX*phav8bo6j-Aef?c=@Ki0F5(3S`DpZtZ&NIOHFJ}*_ zvE`qWN{Saob?9N)mT$EEW7KbBkvmkDe0=!kpOviqsLtLmAH!qefM)&>&9+~xBn~Bf z`Oyce78blS#K+)^x{!p-0dkX=c1V?XZ4CF-MCt_;>G&(--f5 zd_t#QV>5BSdx3k<-HV-}sYLAqZ9puSo>MTs?Fp=&8(mTa!1|YyDzL~~uvuGob zuP1xM_tRL2beW64)Li&I;r-S14+9NIC;|F>d+YlN_{^g-pnRjEpiT#OG>dxC_U_V1 zDADwnNC411^s0RE?a>sSR4^mOicZC|dd^plWOaFW= zl)n_yZK1IF*QIF5+bQ;y1t6$R=UI5hX1$24dRpxyo(F^PknVl%nt}t4|C>JY5IX(BcKsmY{J;A}DjnB|c~$JW13b8WKWmxPlBbx!?`lH|^(^c6?MNKI zzjO~vtOv$0Vk!`o_)fohhIuB{~GGhD#%JNXVzhvK!?+#W-w-Grh zcTh5-SR~;SbQ7yPFM@u`w0(0(gdU`I8nl2Hbv|lm^Z8!tFcGTXr&KZDRj2I3q+*4) z6bgMB(fDLG&9^&+m1sQz$R~~;dSw zyNDU8=8VNmARUzsD0;NX!>`ZGeah#Bs@YnN+%+V1^TLD8XAsaZbh)~})(q%tzNC?Q z4$qbd?I;v5%Xi(T8^518V^!iE64Wa#QWh2KlrDxSB0?Onk1w=1S(i;k{PTXUuErS0 z30ZaX(^G!IxA&FOmvVnR0ISVtCIX(WGd_3H4&r+-B7X3S+{9nKh^s|=XuG9duPiI`=}5K^|XP@Rx# zq3UzDo3JH|{amYT+K&7Y8=K0JMcw@CLXHWsn%SfKqLjaNsbKzr857oS`QE9q!({}1 z8ex_WboIMio02a3@u+SEhuZpDCY_w`*Wl7BYp({9xI<1Db#{#0l7H5tWSsey-}~A> zv)bs5$vl$;H6-L-)8ihKZ0GW1Ks`k@Z!mN5F|OuFvbW9i&C-J_-;%f?NqUe$hxwaw z=${Tbas3hAe>X0@Jx|yzN;?TXJ@?Lps7cj)b}xdijk~YIdEGd3YDTl-co)woRaCZ+ zJ0rX^pQwLp!{D_FHxAJ4LQcN9^JI_~G&cs!y~W3o&>Rq714PkN28Gc3P?5!l_FaiQ z{a&lg?*m+gY3bHQ0Jh`{wQg_uN@7tj&&^)#g0}wO4qpvi0YN`R9jk5%rN-NbP`^4f z8Nk;adYTFlvm`c9>)OHB(s2ef7&rHxT33^}5M>#|W6_$AK7Rps^-UyF@Vbn{4bR&> zlUQ~RNE3}0#S0+FuMfDX_n~i;54Q2?vdR`(^ya|08_J7TlJHO*-|Hya$D9c0*1Vod zzX#LM;po1suT)9qH2;(&KoB(M>5q}m* zqhvu>hFsCz-n*oW4K44a#8cIqqKe7EA?Me6=rW$z5r2Pe%ETcktvAvZ58sc}oF2Rk zBsH5pZCR*@n+XkUh9On?8ySEeOExYxV zEr5ogAoFQjK+onzY2t}4`sJp(XaB)6Okh6_?1QFzE-5j>XKdut*e*NMG7x_UXgdXw zQyc-ueh;y6H$Uf-d;qz?m&NNZ-r#mskb^;_^) zaze?T!7=BvO|+OW5zFz9f37$m?OKYhE`N@-g}H3yw|}7p-ceC-o!iDymoCqv$;JZ` zVoeuAsEOdG89B4AQ$0#Zx>kd>9Ik5T7)73-c%elC%?uBbqfJYfXz^sDg{4eNV^_Jz!>-{F$POF?sBF9(VS}!^RJ_Y+Se|7xy{dl zSx)<&C?p<@A^M?of0TZp)p0jgg&j{`_T3Tq=`3uj)mO&tfZxdGcA(d@{(B$xVboJG+(e+eULe9Q=b-)QKt7)T(@g2G$|51ZO9!St|C zexl|c7-H~$CGI{k-!`U(jwSE|?5`BaQf2K&%b#v`JJdputXY+RNKc^wr;}=Gh&u2s z9FhU-V_AXmVpPQs+%b(eFkZlofQ)Q*`Jw2$fch-ZK&y4Qb%daJP&MB=H8hpR?!EA+ z+3ES|$ZEA#5U{W6{M>@$wMc+a%=qLU>=Ny~->++L6MOO!igQ3UMttX$>|t&sK16ne z3m@NM4XBUW_4Y)!)@mMHM+1en3!%sJGT4pNKz8gKRUqrrfjjX*ZMgl{D6DJetfKKz z>U%2C;PW98+(zKhyhOi_pP^UKClk7^SJFgY0L0{2qF!cG|F|^5FvcJ_bN7q{!u|xK z{kd!Tz>UOr1>#1uyICe{;`UR$3sqw3|Dg*H{3W%bCkIf_5*tFi%qCw;Jv8cIi) z@nGMo?%zHfqZLATb=PhK@1iO75{IXcmc!Y|*WIR7*xhCI6!O_6A5n_c4%g(Jw+VFW zd4M1-{;zPC3QQ)ac?GqT z1&s!{`JgS7>6q(3w}a5Mcee$(lHL!kg3x*xy|j7$2){y?>TYT|vZ5!H@Qnr9K}tVL+-?nXbMHTLCHL}2VKn~WLD zpHv{Y`f|i$cVlaJgB$(+_M=op-Ce>5APnN$;q3^MxY1zvp(al>IAVcLgYhpl;YS+9g36ql!_-_Qh!0U$cOB}AG9X)twHNYN4PfgRD8_H_4b#94A zzi|q;J94CnAc>f9rFY8QZfN)$va~s0&RccVc4yI4$#~L{=D|7qH^aPTk)QXsHa;k* z9?$?^*;T*M;5vQHe(1^c%c>ThL=qr)ckQ0+VeZSVfTAi!l1m|`l1{)uX z$baw72s=-6k@17_a@5=OMi<9x9Wzr3N6(Ptig|I}Wz!SicWbX~EJv_IM93rSP8ACS z@zf&Ud0Jf1K=IKs53p+GBObjga(1WwL`_4BqeM&Wp7}GCkdHp^qd-R=Tf<#ZeDb!! zr3whgZ|?#$cPzd4l<9T4b~^BP##57%k2bbfstG;I6>HagzENQ#5b6R2jWHW`-wbA& zA~Kj2?js!?Pu-J&&-%1memXk{9L;DaGsq;V1vNKB3*F<+3@0Xi;oV{NeJw|Qj1D?!NGq zir&f5p-!UEQLll*59k*ae@y;9$G9kZ?PP+Jnj-NsfW0-xdAIVL_3uQ=bb7W~J}-P@ zsdi69{8qNmJGpqrlIy(F)>QpD+*Of@`2#qWG-m*B_jk`9%HI-H!3{iJxv$Jmf`IrF zSi`;)y_LC#17Blb@sf7q0fNb|o~t!9=gobQ4QFv>S>uawC^QLvR2>+`g>aELCYug4 zy1=ZrD_iyfNU&ILWjM`H?(vw;ZFcuGJtFjB&|*unYlZF-EjIe|(45|v*jm{^+zBSa zqW|4hNz+oK<4T_WZ$2op>6FW9Tg_b*r)ccjXk@q@BI`-20)HS(1V?>l87r&kZME%l|rNJk@qtBv%Eh!L@ous8JtmQe%Q;T{Q-uubB$Byf-HW9VL@y z#u(tBhNn81woe(0&b>sk*{aE9WrZ{Wk{P>WrM6gIY$}_kQ^don72sG?15my|Tc5)C zR3=6;R0I9AN;4%o?{-?oibfgNM-!5bPxYUh!FhMY2qvtgNMzP6I<|JoXNWup8-tG^ zivMUpncSl?A6g7i#YoS-XQz7xMRx3ntht=^rDN-f%4e&&8fspekn9 z&ok_36wYe{C;~`N2^1HUAdvFMeIJgfN0|twt7!1jjeD~-*Y^}71wD5P3#q9L(P69r z>AhmaJ(G&1Nv7D?d=39|N-up2k{wXQ*KR*!kXbf#Mcw}FB@^{~XM~AfD5J;kuZSslNu0#5CSXP$Q-&=Ia-4rJKTUv|#nh9zAWxZ@gunB%;XbeO z^!Abu`2-9_6O!Ny*c*L{sISR_TaWu=-wL&G=*|H*7vP`j#l&FM%3{qjf9H2c>Cbjz zh}}c62cB#V=WsSZD#QaVdZ!}nl$omD7{m-z$~-h7y`FLrEFTHJJM~FVpCTVS6nNmJ z)jae{Aoj`Lnt2u+vO}*Tp~w2*cGpWo8Pe}xR53=pVua4Wh!g&_1~W+iXpt|2K@2O^ zY4@wsJUuHyyT*@#1~X>Ms^w9HKpmgoZ_ynM!vSaW_qR$5pO=jMEUS&1zCwH|jUv77 z*^k^6-qMIGA4nbCDlL6F+we^C)m&XoFa(>s)IsTQQ)g`B&}`r$j>@H-g-`F-ZH<}* z%)www28y7#J?6!mJJvil8Z3KJY9}=bliyJ?5fKypckAe93wL=e4e9kpX8+Kuk^SXN z|2gqKRkuZ~lJCEcVlQ(Kj<^r88emDlyXhDxAlhao)Bg-%v`OUwSNRq$pn0Uan1LRa zt~4Y^A%dT!%xz06MtOh*=e_Ldey!Ib8bp(NSL=BR?iy0)Qm35-A1kHNCW8pc!?lG$ zhrZwZ5Aa|KH&hD>m~Y?Nz~y(W2aZMP$sBy! zzu-4fzf|lMJMyA*u!#iSk;9+DYSu1yeIhe2+e+jvq|7w!^+&Wtq#SBRKAns0oCH=^ zuUFPiXp~`xY)7Wp4aA${+4X0besG0duOjsIwg)6dYgN&d)38y zQ8$IMKlwZOe%eqtCq|d+tINb1|>H93M;1GU94yx^$EVXpktvA!3~WuG-4VXw+&q z#)xt4_bWzBiUOKocZv+rQ?&M%nMT^oqZH>r$^8pkC)S&CH3j}Nh}hi9^~Xj4^0=tu z&+Rhn08!0d33u*|uG4BU(Xp$UH@-Jq(4B!$ro`{v^?R&86#hJaEFtDY9R6Ewst@L4 z#xfT;B7{m#AiXM=B-2OdtozBEdUBIMBHvmSe~ z>~xIX%#r}c>jYI-CsLqcT5c5K;}4INqmCkne#^QQL5aX`{#Ww1mzyuq&eo{$W~0KO z98E;V633nW+-EG4meyI{s~Y=g5*0#F0}Ui72Z>j^dsioZGn6}rve!>f3xr`1^s;CY zlFE_yh0-tg@DTQ$BooW6Hp!PaKZp4?->`-2Kwz`5#&c`9KQ)$6l(OrRlIZPw>;8pj zz3;*+3lMb$)!f<>EKX{!bWnMb=@p`RJW+H0Jc2T1CpM9~QlL7Kf=QPzgwPp43j@Ws zJi)Ec&%%0dMknW@IGP0L+7YXF-pJt>o=72a>TdP~aPsF0&>YZ|;52gK_e;93u@=YE zASl(q&Qwx1aIfGabKU=7{yfoNW8P}8Y9^*Qcs#nXz4B;lT<(}5o{QkOKRMnq6wL4_72Y+%P;0i+;>Jw*5 z5CvoKUG^&7D!jV@f;Zo-D-3)->eyKn`^NHt@4Bl!HsO2X35mj&#w!=_t}z)j{((0#Jv#4DXuEyhbJUy)hEncN+b~f$h$L!dlkPWLqPB0fW2Js;% z$&KAQo+2^Q<2C%tgPT3}>`YflcPTOFry7JJ?F`V;N$>afI}_)?H?VWjL_z9iqzvD# zvcYM_?=05K>L5ofh85HCm}Ajcxu;N~XpN9E*|e*gd|e3L1fx5iBa}3l>Qpwj)4cU} za!dQF55^Ecm)&1L((uouz2>ccel@EQ*ZwbbmLxog!WTUVyqR4*MfPAG zx_siNbUh2X+cjxSACxbIi%zb{1~YgyxK&lomFy%XV5$_{?2qNT@l8)3{2)6hK}`N< z_jJ`^?~(HL|IO}1s%0n|@zO!Wm+W7-{(KC8C#ZY8mH;`I+n)RP_Gyen=IR8|nkOud z4TS(RLKs)Jo-dUYG|nWp$IqSR^(~Kqw=!mi*C@r4>6KvjuxXS0Qcug5J9Z}*xl>>g zilC70R@tw5j5jc8k{FP|;`9t&TcV#>X_=}w9IJzwa41UFR8zY3RVpN`bI>PWAKI9Y zLlLx#TRq>J7Q|QXIltjs;!e3Jj>5BP4&U#%{EAC~$`R-AVddU;;fZGuyeM5SgwS>G zOuruq$xBKdsJ$oY6_o7N;@32r0IvixOYroX@rYr?2`G>~Ys=))t1M9h=)6URw&!?} z%G$l#uxr$IiD6FK-%%5e3qe4)a!q~H#_SI_+tRhcJ~}lcE5|Pm-uhsJfSUJu_>JJ* z;@wc6Kk=4uK)x3%MgUvs!L+eh+Ks5_#oHV5Cx?R_qnM;2IBEm}QM^Z&bt}a*l+xU; zJ?ByF(T+Nhn&jscg1?xdvkii*w&I&T%#^95xVPgIEO^68nagyfaxB1BfM7uG7R-hk z4gydm{hQw1m#X6ARb+!TLoeCOpXF-3vE-F{3hT0YE!lK zGe6&2BGwj{sc1%e?Kb&hhwC_(d1oobIQD)sPjOZCqvwS|gS;AW=~siIO4+^xz*fzg zOV!oGb+q{KSNy-k`M0A}B8^3Ve_eR`fK3pgg2ZboYA3CJTX9&}ugu&knz`zql3V5T zBSmTyix$%5^t%wD%SjaV;yz+ZdHhytSfL`Za6zQhhKZHRA^KB9V@fsi+}N0cr;QJ89&I_m%C{E7DL?|R+g zzMh2}*xNlMPt4x*UvY)L!zF={-^+zP@hp(AhWzxC>ICQ=pI(1Q@Y&MoUo$JiYAJ^Y zfmr$(-;Woi#*bCnr2epq+2?eMar91U3M?o>G~_&XSiE6Ma{Ly-LKYoklUvftm-N1W zH!1CSV3pgiBj3>uX_cj$_kwWc2s`IT(B1Q6sYJ>g`G`&x8%{!Asug_|N+V4A92%z-}tM%3jeXoS$SuN+{|!Qy@}n7pc@rzd_-}CHbYQfULx)O3I1U*joJ_Z6cbuV5kZZ zLR7w9pzJGW<;tANH~x1Q<(Y$uq5;unt=g)pFB567t6P@OBv(l z&97YrBL7XD&gLTi@^VFw7c#Ca+t=-}WKdBw09b{JO&9d|6YT2ILGyyb4Y)Sp0*cT+ z^ipE00aUDP%s7vaeyepE{$xLGIR~Rj+qWX+4_q*SsTJ`;i)sXPr^Rva+hgieVbCVH zvZb1>O`KeinSY1?OFG~34&H1fLd|rMHolZUu{((4{Py@-X1(pVCuq{CDP3~M7kW5i z3$P_)4rpZV=pAqDX~hU4BN8QYuRiwc)VpsVwJ_;XPhXp4F*F7Mh&3@pbsjV7`SczmbkG zA6~Xe+Wlx8iUsUPKIVwgK)ey8xKLkK@}W9_f{%7T>gKXOs~i!U0LGd}xJYI^psN|; z7>h)Jh9L2Edf4^9gM6s;4$wmymP61CQFQPjk1v1yO&g)jWeW7~2z zaxK+ytpxxcQ8)Kb)i%8_DJf*Y!>X^C5|cOoYLhwx{~HFzs{zH&y}mrKeijfwSbV|* z@5V?t8?dLTR5~Y5>I{DaAG|Q(Jv*Ysxi!_MhX9%rjekK(HQz~q?%dh|jt9l6o8Q8h z7Zj|0Gj-$}4Fy)R9qNAqJg0H1s`JG_(tLc!4fvKrg7LF&4FpO1|*;Lw^nZQgWC;k}uPF;XySh#m%}CR2Q7@jG2t2`Ynv z998JJTrMW9Dmn&YLXo3BxGvAJ3#{HGY6?eEu|ywmqf-k{K;*}RiJ1IFkH*=#E%*un z2s$sq=^KG6{mUt7xtF8Ol+VCRwP3+8TIxea<74?)QdxKU%~-XC&U@6Lrw(E~Knpf# zh*j{be+r5=myhdP#(QE9e!Od9BdcAEWw{^QDIdWQ03u=cS`$I3)oPN70-nR=E5}z9 zHU7mDC!o$RN7lXyv3fp4fEg1=Y_R=|pm@4}=jYs?{Lzsn6PmoNEJt;tx+@s?U;&y{6qV9#L^yMm*+Q3X| zqPH*87Eb0Y0^MpsT;H=1R47qlj0u9XZO7}pHTE8{S~kQq`L0s_T`nT1AicEQUNPd~ z4c~PEHc^)oa0ah^u*m??eXxVSztp31hVWgx<^qxq_sj<&D85rL$VMIcFjUBv&c6I4 zdj%jyM|n2lMg!&OBcxHFmOd*$r(x?r6nBvW2sh~|aep5CuK%mcs%$X@XNpXckd6)W z$Y+ptz1s4{YAHmKy{?1?99n}h(}xmGfH=1q_!!6zu&{t5e}GvraMH7ll_+Wg=hV^D z?E|%$BS;1|J)_1;C=KR(y@8%cfFO)L}kgwkJzsa zD+Dn#8#C0#A3SbZ_W(!-qJG2{8m#D917=JRRsLk-3G+AcNH_-dP=W>zMFtw!B=7%9 zd#1!d`9!TtIqjWtpXW!=Bs(kyT*>lVN@@Nt?Z=|8wUqS$nTP%Fnra8qOaxhnnKYAQ;M%x#eS;TKH{s=TB z4NjAhf-E>CNhbXG{kvmG$BVsMboTfg^b}Es%Wx@2#969_hJB7RHWT?fZxWxXwOqnuG#Rvko4N1fN79|bj&Xo{Wf)?5Vb;31=u*Dq0P4zia*1?V#oImlwza#cDV*K4L64V3~DfQS2b2CD(RO68- zVo5BZcnXmUAPF(+8>@!C6|d{wxXpaZ0>lc1dtPqs`|%3Yj#EQ3fc>5DE}*DKNDS{R zu4k=5C-0y9fPoO>z4=gfYa8mYx!QIMWNfh6UdXHJ1dCYb78jl% zp}`?{DnOBn^23OGz0G;Xd|><>Tt-nYL|e~?>aVr=onT)YlS@>s{&8li=cHlb%U=eo z&md-}CgI-2jk=GF3X<7Wn@bY;S+I+TerkMtNUX9zG=~ruFc({ZqBe0cykTFH1eIAb zTf^-f(D5vO^G*4EMLssjdN3ZDl^ExEmJp=fD%?A)_<1GfmptzB?_Y2zf?D?~vAF|W zj#NSybx?0jm$~T(3*+TM!F?J|>gJdi9AqeD0tNAQ1uc{sdDezVU%OWvMvbR2xA2rk zRX?O7kc6)?y^b6x%{_N1`8BU3rH1B2CTiuB&?^hV+JOC|8;GJ8kmoIQVuffi`r>8M zCuTMh>W6cES@d|e$IWHe>j(;+M`A?(=nGFkP&5(K*Au4f^Tw?x?54+E_xkuy(kB*H zlB?dk7-sRCXYjr5g9Tw%e&Vh^gL>DRgUPP;ZY?>?Hq*l=PbfXe9luS*Ky0uZq0Nbd znsj4B^=3adCdVS%ui|ov(@mFBB>iWGripIe0n$O5a8h`b*^hrG;67$N7Klv1!CjIjebdgf1l3MrSATm+t zO}!TiE=-#qo&5NL!2@`G>IoJ2l-2j1Pipq*rWDy#KI(hhan#T-THwKVHI!wtN-mpb zRFP8e%iq+x$rLZ3;nfY=DcV;Sxn3f?B*(vvh0}M5?eWEDE8n~odi!N&)u0Y^yg9wD zq&6B;O}Aa}%2^t_W6g*cMvO;7o?0?!`y_=VP_x49Pm1>TJ?Y=4#!9Qoe46r_D>sfJkEf%wz)AZ2rmE7jc0$c6D9e}m^{AaKwc#RE7zS`WuVPsq*cQW!g12c0+k-vj zt5@idIH37Fz7n0W@%41DK_7b?}MDBWE?hE&ocIZ#2=O>r~*{AwA60#_#0|kEYgGusmCXB4@cv z%ckJF^U<(GU3s>ozK#WuW&kh#_$0F(xtb9QN+)Gtj%;`(0)^5=4OlQHC2RJ1l({z z6sPl8@Q}TOoik%5a8BM@|FL*!E|KlRT)|Ix-XhWJuQ;9yYNoZ(8kC*y->YAkGpHf# z;{GFK@Yp`x$|u>~^GH%ogF@annGil$b-obwin-$4J^&e*o}XZNc_`u!FG8i_?vT7% z4s7|>QnsM(%gr@?&H5r|?EXWUFz*r3{he+}SjDS4GizRjVSMcFDz_X)-3kq8 z4C-z5@?B9oRsYy4)4BfKW4)0obs?1LaZNiBVyWmHAGhBec_Z9>ZRW|Kq0siGxfSt~ zi{U^`Qt8@->9e~WhUV_1NW7>za{^{SRzMhYR?nEyw$mW>a_aOcY1wyYLZ9TkM{4Q9Kt|C z-hK`OtYS7u`s0bJa4Cp!yzG_qIY8=8cFU=zXT{GjP8Zv$g@5k4ijijze}=!6kl8sA zCacwa|B^koSM6%1A)S4fThy9WBowi|vSnpXqnNfRuGL)b(~MJ3DBRDXqWt`4ctpGy zSJP2-Yq90a)NNF@5P5S+X;n&M}Q*!JV-hUYx!f*eoGN;CXq! zwr4h=tvyCEIV*nzFf!o zcV{@}>y=&A@H4&qPigW6bOp(bT6+|9J=q$%D94~~_2*CCum$``}p;9{zAe&rB@hyq!Jon)U6!u6kci^Qs3pl)0w18wRdbJ@3_CU*r9&Y?XbDD@SSR-e(1Uo2UiZTaq51 zvfZqMSCs4Js{);8p2dX9xc5(Cwa`uKeC?^zeB4{HA{qMs~yhjh({s1Zy|otx$f>2 zPrR9u4|bnQyA}pnOG&51IvUhH`+UDB-HelmNghS0(4bn*c&jOb zUzp5V8nKx;KlHOZ42()Y8f(3B*0;9+PGOaRrq8&8YFbg$FGm~^eAnSe9Dn#L^CSoj zhtKr;w0_fucYb~4lym~d9W2hg&LLuVv#yy0q2c_8scSAMstTvo0Krg>ffYWcXc`jS zbNWcL@Y}`>I2i_aoJQO48@clTeF&jXq;ip02;sebw{M(*n~>?j@$GN4cBda0v6H_( z0MJiq2%tXSanpx-zD)U@OCtgR63%OcVbE~gfE;o%XfRC&O@tDP(dqkm0PZ10BXP$u z9yFhHz;YP;oFr(W1(tA8AVA2$vm(|08zU20{^Qv3QPLO(x3m&OIn*5wdEo}zO5`&Q z^!9&@xRI(RKUR#n(LYXW?j51A`XD`cIT_nY1xM9LowOYM1z$v3{$n@vW$Y14P$%2* zfQgLOHn;R~uqBRX#1{vzP!_D7Zg+?5Dsz-<9A*|7ldAl%7H?&N1M?|FP0#;`{He=$DJ1IL*E zZ?akg1l(V!o&X1B;7)eJC@tK{&aW(aR|{r&@-O0;fCQh8Sh6#_;7SiH-_3V)qQ`O` zH%7uktv*;7^IbGw2t+nfBMDa0!wg)7;9^sCGyyzjJ)2Nd{~zixmrk+A5g53I{v|}e z)dsiFOLJjIVK9U~*44Xp0dAl~pTI`i9Pz_!c(h=fqp^>`R+wHr;xElI%#H(_05c%B z|Hd7Ufw1e;(J5bn|9owp0`O mx^(dW3;=`wfB!l)djveoY51*O?Slj0KYCh5njh7#ll~8&4nJQ2 literal 0 HcmV?d00001