diff --git a/src/main/java/fr/pasteur/ida/zellige/utils/LocalMinimumDetection.java b/src/main/java/fr/pasteur/ida/zellige/utils/LocalMinimumDetection.java new file mode 100644 index 0000000000000000000000000000000000000000..f54f264769c5bfa141305e0a7854ad12376f6664 --- /dev/null +++ b/src/main/java/fr/pasteur/ida/zellige/utils/LocalMinimumDetection.java @@ -0,0 +1,173 @@ +package fr.pasteur.ida.zellige.utils; + + +import fr.pasteur.ida.zellige.surfaceConstruction.element.Pixels; +import net.imglib2.RandomAccess; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.img.Img; +import net.imglib2.img.ImgFactory; +import net.imglib2.img.display.imagej.ImageJFunctions; +import net.imglib2.type.NativeType; +import net.imglib2.type.numeric.RealType; +import net.imglib2.util.ImgUtil; +import net.imglib2.view.IntervalView; +import net.imglib2.view.Views; + + +/** + * @author ctrebeau + * This class is used to select local maximums of a stack image (3 dimensions) by using the static method + * findMaximums(). + * First a gaussian convolution (anisotropic ) is performed then for each orthogonal section (XZ section) + * the local maximums are found by the determination of the partial derivative in the section first dimension + * AKA dimension Y. + */ + +public class LocalMinimumDetection +{ + + /** + * Finds all the local maximums of an input stack image, + * according to given local and global thresholds by first using a gaussian blur. + * The local maximums are stored in a 2D {@link Pixels} array. + * + * @param zStack - the z-stack to analysed. + * @param <T> - the type of the z-stack. + */ + + + public static < T extends RealType< T > & NativeType< T > > Img< T > + findMinimums( RandomAccessibleInterval< T > zStack, ImgFactory< T > factory ) + { + Img< T > minimums = factory.create( zStack.dimension( 0 ), zStack.dimension( 1 ), zStack.dimension( 2 ) ); +// ImgUtil.copy( zStack, minimums ); + // For each individual orthogonal XZ sections.*/ + for ( int x = 0; x <= zStack.dimension( 0 ) - 1; x++ )// iteration on dimension X + { + try + { + // The YZ section. + IntervalView< T > intervalView = Views.hyperSlice( zStack, 0, x ); +// if (x == 15) +// { +// ImageJFunctions.show(intervalView); +// } + // Each maximum found is stored in double Pixel array . + findListOfMinimum( intervalView, factory, x, minimums ); + } + catch ( Exception e ) + { + e.printStackTrace(); + } + } + return minimums; + } + + + /** + * Finds all the local maximums of an 2D image, + * according to given local and global thresholds by first using a gaussian blur and then a partial derivative + * The local maximums are stored in a 2D {@link Pixels} array. + * + * @param YZProjection - An orthogonal section (2 dimensions). + * @param x - the index of the orthogonal section AKA the stack index in dimension X. + * @param <T> - the YZProjection type. + */ + private static < T extends RealType< T > & NativeType< T > > void + findListOfMinimum( RandomAccessibleInterval< T > YZProjection, ImgFactory< T > factory, int x, + RandomAccessibleInterval< T > maximums ) + { + // Output for partial derivative. + // Addition of one supplementary line, in case the last slice contains some local maximums + RandomAccessibleInterval< T > partialDerivative = factory.create( YZProjection.dimension( 0 ) + 1, + YZProjection.dimension( 1 ) + 1 ); + // Partial derivative computation */ + Utils.derivative( YZProjection, partialDerivative ); +// if (x == 15) +// { +// ImageJFunctions.show(partialDerivative); +// } + getListOfLocalMinimum( maximums, partialDerivative,YZProjection, x ); + } + + + /** + * @param max the output image containing the local maximums as a {@link RandomAccessibleInterval} + * @param partialDerivative the x-positioned partial derivative section of the input image as a {@link RandomAccessibleInterval} + * @param x the index in dimension X + * @param <T> the type of both {@link RandomAccessibleInterval} + */ + private static < T extends RealType< T > & NativeType< T > > void + getListOfLocalMinimum + ( RandomAccessibleInterval< T > max, + RandomAccessibleInterval< T > partialDerivative, + int x + ) + { + RandomAccess< T > derivativeAccess = partialDerivative.randomAccess(); + RandomAccess< T > maxAccess = max.randomAccess(); + + for ( int y = 0; y <= max.dimension( 1 ) - 1; y++ ) + { + derivativeAccess.setPosition( y, 0 ); + for ( int z = 0; z <= max.dimension( 2 ) - 1; z++ ) + { + if ( isALocalMinimum( derivativeAccess, z ) ) + { + Utils.setPosition( maxAccess, x, y, z ); + maxAccess.get().setOne(); + } + } + } + } + + private static < T extends RealType< T > & NativeType< T > > void + getListOfLocalMinimum + ( RandomAccessibleInterval< T > minimums, + RandomAccessibleInterval< T > partialDerivative, + RandomAccessibleInterval <T> original, + int x + ) + { + RandomAccess< T > derivativeAccess = partialDerivative.randomAccess(); + RandomAccess< T > minAccess = minimums.randomAccess(); + RandomAccess< T > oriAccess = original.randomAccess(); + + for ( int y = 0; y <= minimums.dimension( 1 ) - 1; y++ ) + { + derivativeAccess.setPosition( y, 0 ); + for ( int z = 0; z <= minimums.dimension( 2 ) - 1; z++ ) + { + if ( isALocalMinimum( derivativeAccess, z ) ) + { + Utils.setPosition( oriAccess, y, z ); + Utils.setPosition( minAccess, x, y, z ); + minAccess.get().set( oriAccess.get() ); +// maxAccess.get().setOne(); + } + } + } + } + + /** + * Method to find if a pixel of a given stack is a local maximum. + * <p> + * (X or Y for respectively XZ sections and YZ sections). + * + * @param derivative - the first {@link RandomAccess} for the partial derivative image + * set at index v for dimension 0. + * @param z - the index of the Z dimension + * @param <T> - the type of the image + * @return - true if the position (v, z) is a local maximum false otherwise. + */ + + private static < T extends RealType< T > & NativeType< T > > boolean isALocalMinimum + ( RandomAccess< T > derivative, int z ) + { + derivative.setPosition( z, 1 ); + final double tZero = derivative.get().getRealDouble(); + derivative.setPosition( ( z + 1 ), 1 ); + final double tOne = derivative.get().getRealDouble(); + return ( tZero <0 && tOne >= 0 ); + } +}