net.imglib2.type.numeric.complex.ComplexFloatType Java Examples

protected < T extends RealType< T > > ContentBased< T > getContentBased( final RandomAccessibleInterval< T > img, final ViewDescription desc, final ImgLoader imgLoader )
{
	final double[] sigma1 = ProcessFusion.defaultContentBasedSigma1.clone();
	final double[] sigma2 = ProcessFusion.defaultContentBasedSigma2.clone();

	final double minRes = getMinRes( desc, imgLoader );
	final VoxelDimensions voxelSize = ViewSetupUtils.getVoxelSizeOrLoad( desc.getViewSetup(), desc.getTimePoint(), imgLoader );

	if ( ProcessFusion.defaultAdjustContentBasedSigmaForAnisotropy )
	{
		for ( int d = 0; d < 2; ++d )
		{
			sigma1[ d ] /= voxelSize.dimension( d ) / minRes;
			sigma2[ d ] /= voxelSize.dimension( d ) / minRes;
		}
	}

	return new ContentBased<T>( img, bb.getImgFactory( new ComplexFloatType() ), sigma1, sigma2);
}
 

protected void assertComplexImagesEqual(final Img<ComplexFloatType> img1,
	final Img<ComplexFloatType> img2, final float delta)
{
	final Cursor<ComplexFloatType> c1 = img1.cursor();
	final Cursor<ComplexFloatType> c2 = img2.cursor();

	while (c1.hasNext()) {
		c1.fwd();
		c2.fwd();

		// assert that the inverse = the input within the error delta
		assertEquals(c1.get().getRealFloat(), c2.get().getRealFloat(), delta);
		// assert that the inverse = the input within the error delta
		assertEquals(c1.get().getImaginaryFloat(), c2.get().getImaginaryFloat(),
			delta);
	}

}
 

/**
 * test that a forward transform followed by an inverse transform gives us
 * back the original image
 */
@Test
public void testFFT3DOp() {
	final int min = expensiveTestsEnabled ? 115 : 9;
	final int max = expensiveTestsEnabled ? 120 : 11;
	for (int i = min; i < max; i++) {

		final long[] dimensions = new long[] { i, i, i };

		// create an input with a small sphere at the center
		final Img<FloatType> in = generateFloatArrayTestImg(false, dimensions);
		placeSphereInCenter(in);

		final Img<FloatType> inverse = generateFloatArrayTestImg(false,
			dimensions);

		@SuppressWarnings("unchecked")
		final Img<ComplexFloatType> out = (Img<ComplexFloatType>) ops.run(
			FFTMethodsOpF.class, in);
		ops.run(IFFTMethodsOpC.class, inverse, out);

		assertImagesEqual(in, inverse, .00005f);
	}

}
 

@Test
public void testCreateNativeType() {

	// default
	Object type = ops.run(DefaultCreateNativeType.class);
	assertEquals(type.getClass(), DoubleType.class);

	// FloatType
	type = ops.run(CreateNativeTypeFromClass.class, FloatType.class);
	assertEquals(type.getClass(), FloatType.class);

	// ComplexFloatType
	type = ops.run(CreateNativeTypeFromClass.class, ComplexFloatType.class);
	assertEquals(type.getClass(), ComplexFloatType.class);

	// DoubleType
	type = ops.run(CreateNativeTypeFromClass.class, DoubleType.class);
	assertEquals(type.getClass(), DoubleType.class);

	// ComplexDoubleType
	type = ops.run(CreateNativeTypeFromClass.class, ComplexDoubleType.class);
	assertEquals(type.getClass(), ComplexDoubleType.class);

}
 

@Override
@SuppressWarnings({ "unchecked", "rawtypes" })
public void initialize() {
	super.initialize();

	// if no type was passed in the default is ComplexFloatType
	if (fftType == null) {
		fftType = (C)ops().create().nativeType(ComplexFloatType.class);
	}

	padOp = (BinaryFunctionOp) Functions.binary(ops(), PadInputFFTMethods.class,
		RandomAccessibleInterval.class, RandomAccessibleInterval.class,
		Dimensions.class, fast);

	createOp = (UnaryFunctionOp) Functions.unary(ops(),
		CreateOutputFFTMethods.class, RandomAccessibleInterval.class,
		Dimensions.class, fftType, fast);

	fftMethodsOp = (UnaryComputerOp) Computers.nullary(ops(),
		FFTMethodsOpC.class, RandomAccessibleInterval.class,
		RandomAccessibleInterval.class);

}
 

/**
 * Convenience wrapper to create an {@link Img} of type {@link ComplexFloatType}
 * with a Gabor kernel.
 */
@OpMethod(
	op = net.imagej.ops.create.kernelGabor.DefaultCreateKernelGabor.class)
public RandomAccessibleInterval<ComplexFloatType>
	kernelGaborComplexFloat(final double[] sigmas, final double... period)
{
	@SuppressWarnings("unchecked")
	final RandomAccessibleInterval<ComplexFloatType> result =
		(RandomAccessibleInterval<ComplexFloatType>) ops().run(
			net.imagej.ops.create.kernelGabor.DefaultCreateKernelGabor.class,
			sigmas, period, new ComplexFloatType());
	return result;
}
 

protected Img< FloatType > approximateEntropy(
		final RandomAccessibleInterval< FloatType > input,
		final ImgFactory< ComplexFloatType > imgFactory,
		final double[] sigma1,
		final double[] sigma2 )

{
	// the result
	ImgFactory<FloatType> f;
	try { f = imgFactory.imgFactory( new FloatType() ); } catch (IncompatibleTypeException e) { f = new ArrayImgFactory< FloatType >(); }
	
	final Img< FloatType > conv = f.create( input, new FloatType() );
	
	// compute I*sigma1
	FFTConvolution< FloatType > fftConv = new FFTConvolution<FloatType>( input, createGaussianKernel( sigma1 ), conv, imgFactory );
	fftConv.convolve();
	
	// compute ( I - I*sigma1 )^2
	final Cursor< FloatType > c = conv.cursor();
	final RandomAccess< FloatType > r = input.randomAccess();
	
	while ( c.hasNext() )
	{
		c.fwd();
		r.setPosition( c );
		
		final float diff = c.get().get() - r.get().get();
		c.get().set( diff * diff );
	}
	
	// compute ( ( I - I*sigma1 )^2 ) * sigma2
	fftConv = new FFTConvolution<FloatType>( conv, createGaussianKernel( sigma2 ), imgFactory );
	fftConv.convolve();

	// normalize to [0...1]
	FusionHelper.normalizeImage( conv );

	return conv;
}
 

public ContentBased(
		final RandomAccessibleInterval< T > input,
		final ImgFactory< ComplexFloatType > imgFactory,
		final double[] sigma1,
		final double[] sigma2 )
{
	this.n = input.numDimensions();
	
	this.contentBasedImg = approximateEntropy(
			new ConvertedRandomAccessibleInterval< T, FloatType >( input, new RealFloatConverter< T >(),  new FloatType() ),
			imgFactory,
			sigma1,
			sigma2 );		
}
 

protected RandomAccessibleInterval<C> getFFTInput() {
	if (fftType == null) {
		fftType = (ComplexType<C>) ops().create().nativeType(
			ComplexFloatType.class);
	}
	return fftInput;
}
 

@Override
@SuppressWarnings({ "unchecked", "rawtypes" })
public void initialize() {
	super.initialize();

	/**
	 * Op used to pad the input
	 */
	padOp = (BinaryFunctionOp) Functions.binary(ops(), PadInputFFTMethods.class,
		RandomAccessibleInterval.class, RandomAccessibleInterval.class,
		Dimensions.class, true, obfInput);

	/**
	 * Op used to pad the kernel
	 */
	padKernelOp = (BinaryFunctionOp) Functions.binary(ops(),
		PadShiftKernelFFTMethods.class, RandomAccessibleInterval.class,
		RandomAccessibleInterval.class, Dimensions.class, true);

	if (fftType == null) {
		fftType = (ComplexType<C>) ops().create().nativeType(
			ComplexFloatType.class);
	}

	/**
	 * Op used to create the complex FFTs
	 */
	createOp = (UnaryFunctionOp) Functions.unary(ops(),
		CreateOutputFFTMethods.class, RandomAccessibleInterval.class,
		Dimensions.class, fftType, true);
}
 

/**
 * Convenience wrapper to create an {@link Img} of type {@link ComplexFloatType}
 * with an isotropic Gabor kernel.
 */
@OpMethod(
	op = net.imagej.ops.create.kernelGabor.CreateKernelGaborIsotropic.class)
public RandomAccessibleInterval<ComplexFloatType>
	kernelGaborComplexFloat(final Double sigma, final double... period)
{
	@SuppressWarnings("unchecked")
	final RandomAccessibleInterval<ComplexFloatType> result =
		(RandomAccessibleInterval<ComplexFloatType>) ops().run(
			net.imagej.ops.create.kernelGabor.CreateKernelGaborIsotropic.class,
			sigma, period, new ComplexFloatType());
	return result;
}
 

@OpMethod(op = net.imagej.ops.convert.ConvertTypes.ComplexToCfloat32.class)
public <C extends ComplexType<C>> ComplexFloatType cfloat32(
	final ComplexFloatType out, final C in)
{
	final ComplexFloatType result = (ComplexFloatType) ops().run(
		Ops.Convert.Cfloat32.class, out, in);
	return result;
}
 

@OpMethod(op = net.imagej.ops.convert.ConvertImages.Cfloat32.class)
public <C extends ComplexType<C>> Img<ComplexFloatType> cfloat32(
	final Img<ComplexFloatType> out, final IterableInterval<C> in)
{
	@SuppressWarnings("unchecked")
	final Img<ComplexFloatType> result = (Img<ComplexFloatType>) ops().run(
		Ops.Convert.Cfloat32.class, out, in);
	return result;
}
 

@OpMethod(op = net.imagej.ops.convert.ConvertImages.Cfloat32.class)
public <C extends ComplexType<C>> Img<ComplexFloatType> cfloat32(
	final IterableInterval<C> in)
{
	@SuppressWarnings("unchecked")
	final Img<ComplexFloatType> result = (Img<ComplexFloatType>) ops().run(
		Ops.Convert.Cfloat32.class, in);
	return result;
}
 

public static void main(String[] args) {
	
	new ImageJ();
	
	Img<FloatType> img1 = ImgLib2Util.openAs32Bit(new File("src/main/resources/img1singleplane.tif"));
	Img<FloatType> img2 = ImgLib2Util.openAs32Bit(new File("src/main/resources/img2singleplane.tif"));
	
	ExecutorService service = Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors());
	
	RandomAccessibleInterval<FloatType> pcm = calculatePCM(img1, img2, new ArrayImgFactory<FloatType>(), new FloatType(),
			new ArrayImgFactory<ComplexFloatType>(), new ComplexFloatType(), service );
	
	PhaseCorrelationPeak2 shiftPeak = getShift(pcm, img1, img2);
	
	RandomAccessibleInterval<FloatType> res = PhaseCorrelation2Util.dummyFuse(img1, img2, shiftPeak,service);
			
	ImageJFunctions.show(res);
}
 

@OpMethod(op = net.imagej.ops.convert.ConvertTypes.ComplexToCfloat32.class)
public <C extends ComplexType<C>> ComplexFloatType cfloat32(final C in) {
	final ComplexFloatType result = (ComplexFloatType) ops().run(
		Ops.Convert.Cfloat32.class, in);
	return result;
}
 

@Override
public void compute(final C input, final ComplexFloatType output) {
	output.set(input.getRealFloat(), input.getImaginaryFloat());
}
 

/**
 * test the fast FFT
 */
@Test
public void testFastFFT3DOp() {

	final int min = expensiveTestsEnabled ? 120 : 9;
	final int max = expensiveTestsEnabled ? 130 : 11;
	final int size = expensiveTestsEnabled ? 129 : 10;
	for (int i = min; i < max; i++) {

		// define the original dimensions
		final long[] originalDimensions = new long[] { i, size, size };

		// arrays for the fast dimensions
		long[] fastDimensions = new long[3];
		long[] fftDimensions;

		// compute the dimensions that will result in the fastest FFT time
	
		long[][] temp=ops.filter().fftSize(new FinalDimensions(originalDimensions), true);
		fastDimensions=temp[0];
	
		// create an input with a small sphere at the center
		final Img<FloatType> inOriginal = generateFloatArrayTestImg(false,
			originalDimensions);
		placeSphereInCenter(inOriginal);

		// create a similar input using the fast size
		final Img<FloatType> inFast = generateFloatArrayTestImg(false,
			fastDimensions);
		placeSphereInCenter(inFast);

		// call FFT passing false for "fast" (in order to pass the optional
		// parameter we have to pass null for the
		// output parameter).
		@SuppressWarnings("unchecked")
		final RandomAccessibleInterval<ComplexFloatType> fft1 =
			(RandomAccessibleInterval<ComplexFloatType>) ops.run(
				FFTMethodsOpF.class, inOriginal, null, false);

		// call FFT passing true for "fast" The FFT op will pad the input to the
		// fast
		// size.
		@SuppressWarnings("unchecked")
		final RandomAccessibleInterval<ComplexFloatType> fft2 =
			(RandomAccessibleInterval<ComplexFloatType>) ops.run(
				FFTMethodsOpF.class, inOriginal, null, true);

		// call fft using the img that was created with the fast size
		@SuppressWarnings("unchecked")
		final RandomAccessibleInterval<ComplexFloatType> fft3 =
			(RandomAccessibleInterval<ComplexFloatType>) ops.run(
				FFTMethodsOpF.class, inFast);

		// create an image to be used for the inverse, using the original
		// size
		final Img<FloatType> inverseOriginalSmall = generateFloatArrayTestImg(
			false, originalDimensions);

		// create an inverse image to be used for the inverse, using the
		// original
		// size
		final Img<FloatType> inverseOriginalFast = generateFloatArrayTestImg(
			false, originalDimensions);

		// create an inverse image to be used for the inverse, using the
		// fast size
		final Img<FloatType> inverseFast = generateFloatArrayTestImg(false,
			fastDimensions);

		// invert the "small" FFT
		ops.run(IFFTMethodsOpC.class, inverseOriginalSmall, fft1);

		// invert the "fast" FFT. The inverse will should be the original
		// size.
		ops.run(IFFTMethodsOpC.class, inverseOriginalFast, fft2);

		// invert the "fast" FFT that was acheived by explicitly using an
		// image
		// that had "fast" dimensions. The inverse will be the fast size
		// this
		// time.
		ops.run(IFFTMethodsOpC.class, inverseFast, fft3);

		// assert that the inverse images are equal to the original
		assertImagesEqual(inverseOriginalSmall, inOriginal, .0001f);
		assertImagesEqual(inverseOriginalFast, inOriginal, .00001f);
		assertImagesEqual(inverseFast, inFast, 0.00001f);
	}
}
 

@Override
public ComplexFloatType createOutput(final C input) {
	return new ComplexFloatType();
}
 

@Test
public void testPCNegativeShift() {
	
	// TODO: very large shifts (nearly no overlap) lead to incorrect shift determination (as expected)
	// maybe we can optimize behaviour in this situation
	Img< FloatType > img = ArrayImgs.floats( 200, 200 );
	Random rnd = new Random( seed );
	
	for( FloatType t : img )
		t.set( rnd.nextFloat());
	
	long shiftX = -2;
	long shiftY = -2;
	
	FinalInterval interval1 = new FinalInterval(new long[] {50, 50});
	FinalInterval interval2 = Intervals.translate(interval1, shiftX, 0);
	interval2 = Intervals.translate(interval2, shiftY, 1);

	int [] extension = new int[img.numDimensions()];
	Arrays.fill(extension, 10);
	
	RandomAccessibleInterval<FloatType> pcm = PhaseCorrelation2.calculatePCM(
			Views.zeroMin(Views.interval(Views.extendZero( img ), interval1)),
			Views.zeroMin(Views.interval(Views.extendZero( img ), interval2)),
			extension, new ArrayImgFactory<FloatType>(), 
			new FloatType(), new ArrayImgFactory<ComplexFloatType>(), new ComplexFloatType(), Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors()));
	
	PhaseCorrelationPeak2 shiftPeak = PhaseCorrelation2.getShift(pcm,
			Views.zeroMin(Views.interval(Views.extendZero( img ), interval1)),
			Views.zeroMin(Views.interval(Views.extendZero( img ), interval2)),
			20, 0, false, Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors()));
	
	long[] expected = new long[]{shiftX, shiftY};
	long[] found = new long[img.numDimensions()];
	
	
	
	shiftPeak.getShift().localize(found);
	System.out.println( Util.printCoordinates( found ) );
	
	assertArrayEquals(expected, found);
	
}
 

public static void main( String[] args ) throws IncompatibleTypeException
{
	new ImageJ();
	
	ImagePlus imp = new ImagePlus( "/Users/preibischs/workspace/TestLucyRichardson/src/resources/dros-1.tif" );
	
	Img< FloatType > img = ImageJFunctions.wrap( imp );

	final double[] sigma1 = new double[]{ 20, 20 };
	final double[] sigma2 = new double[]{ 30, 30 };
	
	ContentBased< FloatType > cb = new ContentBased<FloatType>( img, img.factory().imgFactory( new ComplexFloatType() ), sigma1, sigma2 );
	
	ImageJFunctions.show( cb.getContentBasedImg() );
}
 

@Test
public void testPC() {
	
	// TODO: very large shifts (nearly no overlap) lead to incorrect shift determination (as expected)
	// maybe we can optimize behaviour in this situation
	Img< FloatType > img = ArrayImgs.floats( 200, 200 );
	Random rnd = new Random( seed );
	
	for( FloatType t : img )
		t.set( rnd.nextFloat());
	
	long shiftX = 28;
	long shiftY = 0;
	
	FinalInterval interval1 = new FinalInterval(new long[] {50, 50});
	FinalInterval interval2 = Intervals.translate(interval1, shiftX, 0);
	interval2 = Intervals.translate(interval2, shiftY, 1);

	
	int [] extension = new int[img.numDimensions()];
	Arrays.fill(extension, 10);
	
	RandomAccessibleInterval<FloatType> pcm = PhaseCorrelation2.calculatePCM(Views.interval(img, interval1), Views.zeroMin(Views.interval(img, interval2)), extension, new ArrayImgFactory<FloatType>(), 
			new FloatType(), new ArrayImgFactory<ComplexFloatType>(), new ComplexFloatType(), Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors()));
	
	PhaseCorrelationPeak2 shiftPeak = PhaseCorrelation2.getShift(pcm, Views.interval(img, interval1), Views.zeroMin(Views.interval(img, interval2)), 20, 0, false, Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors()));
	
	long[] expected = new long[]{shiftX, shiftY};
	long[] found = new long[img.numDimensions()];
	
	
	
	shiftPeak.getShift().localize(found);
	
	assertArrayEquals(expected, found);
	
}
 

@Test
public void testPCRealShift() {
	
	// TODO: very large shifts (nearly no overlap) lead to incorrect shift determination (as expected)
	// maybe we can optimize behaviour in this situation
	Img< FloatType > img = ArrayImgs.floats( 200, 200 );
	Random rnd = new Random( seed );
	
	for( FloatType t : img )
		t.set( rnd.nextFloat());
	
	double shiftX = -20.9;
	double shiftY = 1.9;
	
	// to test < 0.5 px off
	final double eps = 0.5;
	
	FinalInterval interval2 = new FinalInterval(new long[] {50, 50});
	
	

	AffineRandomAccessible< FloatType, AffineGet > imgTr = RealViews.affine( Views.interpolate( Views.extendZero( img ), new NLinearInterpolatorFactory<>() ), new Translation2D( shiftX, shiftY ));
	IntervalView< FloatType > img2 = Views.interval( Views.raster( imgTr ), interval2);
	
	int [] extension = new int[img.numDimensions()];
	Arrays.fill(extension, 10);
	
	RandomAccessibleInterval<FloatType> pcm = PhaseCorrelation2.calculatePCM(Views.zeroMin(img2), Views.zeroMin(Views.interval(img, interval2)), extension, new ArrayImgFactory<FloatType>(), 
			new FloatType(), new ArrayImgFactory<ComplexFloatType>(), new ComplexFloatType(), Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors()));
	
	PhaseCorrelationPeak2 shiftPeak = PhaseCorrelation2.getShift(pcm, Views.zeroMin(img2), Views.zeroMin(Views.interval(img, interval2)), 20, 0, true, Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors()));
	
	
	double[] expected = new double[]{shiftX, shiftY};
	double[] found = new double[img.numDimensions()];
	
	
	
	
	shiftPeak.getSubpixelShift().localize(found);
	
	System.out.println( Util.printCoordinates( found ) );
	
	
	for (int d = 0; d < expected.length; d++)
		assertTrue( Math.abs( expected[d] - found[d] ) < eps );
	
}