net.imglib2.FinalDimensions Java Examples

@Test
@SuppressWarnings({ "unchecked", "rawtypes" })
public void copyRAIDifferentSizeTest() {

	// create a copy op
	final UnaryHybridCF<IntervalView<UnsignedByteType>, RandomAccessibleInterval<UnsignedByteType>> copy =
		(UnaryHybridCF) Hybrids.unaryCF(ops, CopyRAI.class,
			RandomAccessibleInterval.class, IntervalView.class);

	assertNotNull(copy);

	final Img<UnsignedByteType> out = ops.create().img(new FinalDimensions(
		size2), new UnsignedByteType());

	// copy view to output and assert that is equal to the mean of the view
	copy.compute(view, out);
	assertEquals(ops.stats().mean(out).getRealDouble(), 100.0, delta);

	// also try with a planar image
	final Img<UnsignedByteType> outFromPlanar = ops.create().img(
		new FinalDimensions(size2), new UnsignedByteType());

	copy.compute(viewPlanar, outFromPlanar);
	assertEquals(ops.stats().mean(outFromPlanar).getRealDouble(), 100.0, delta);

}
 

@SuppressWarnings("unchecked")
@Test
public void testImageFactory() {

	final Dimensions dim = new FinalDimensions( 10, 10, 10 );

	assertEquals("Labeling Factory: ", ArrayImgFactory.class,
		((Img<?>) ((ImgLabeling<String, ?>) ops.run(
			DefaultCreateImgLabeling.class, dim, null,
			new ArrayImgFactory<IntType>())).getIndexImg()).factory().getClass());

	assertEquals("Labeling Factory: ", CellImgFactory.class,
		((Img<?>) ((ImgLabeling<String, ?>) ops.run(
			DefaultCreateImgLabeling.class, dim, null,
			new CellImgFactory<IntType>())).getIndexImg()).factory().getClass());

}
 

public static void main(String[] args) {
	
	double o1 = 6;
	double o2 = Double.NEGATIVE_INFINITY;
	int np1 = 30;
	int np2 = 20;

	System.out.println( Double.isInfinite( o2 ));

	int ccCompare = Double.compare(o1, o2);
	if (ccCompare != 0)
		System.out.println( ccCompare );
	else 
		System.out.println( (int)(np1 - np2) );
	
	System.exit( 0 );
	PhaseCorrelationPeak2 peaks = new PhaseCorrelationPeak2(new Point(new int[] {10, 10}), 1.0);
	Dimensions pcmDims = new FinalDimensions(new int[] {50, 50});
	Dimensions imgDims = new FinalDimensions(new int[] {30, 30});
	PhaseCorrelation2Util.expandPeakToPossibleShifts(peaks, pcmDims, imgDims, imgDims);
	
}
 

public void initialize(RandomAccessibleInterval<T> yk_iterated) {
	if (yk_prediction == null) {

		long[] temp = new long[yk_iterated.numDimensions()];
		yk_iterated.dimensions(temp);

		FinalDimensions dims = new FinalDimensions(temp);

		yk_prediction = create.calculate(dims);
		xkm1_previous = create.calculate(dims);
		yk_prediction = create.calculate(dims);
		gk = create.calculate(dims);
		hk_vector = create.calculate(dims);

	}

}
 

/**
 * Creates the List of {@link ViewSetup} for the {@link SpimData} object.
 * The {@link ViewSetup} are defined independent of the {@link TimePoint},
 * each {@link TimePoint} should have the same {@link ViewSetup}s. The {@link MissingViews}
 * class defines if some of them are missing for some of the {@link TimePoint}s
 *
 * @return
 */
protected ArrayList< ViewSetup > createViewSetups( final DHMMetaData meta )
{
	final ArrayList< Channel > channels = new ArrayList< Channel >();
	channels.add( new Channel( meta.getAmpChannelId(), meta.getAmplitudeDir() ) );
	channels.add( new Channel( meta.getPhaseChannelId(), meta.getPhaseDir() ) );

	final ArrayList< Illumination > illuminations = new ArrayList< Illumination >();
	illuminations.add( new Illumination( 0, String.valueOf( 0 ) ) );

	final ArrayList< Angle > angles = new ArrayList< Angle >();
	angles.add( new Angle( 0, String.valueOf( 0 ) ) );

	final ArrayList< ViewSetup > viewSetups = new ArrayList< ViewSetup >();
	for ( final Channel c : channels )
		for ( final Illumination i : illuminations )
			for ( final Angle a : angles )
			{
				final VoxelDimensions voxelSize = new FinalVoxelDimensions( meta.calUnit, meta.calX, meta.calY, meta.calZ );
				final Dimensions dim = new FinalDimensions( new long[]{ meta.getWidth(), meta.getHeight(), meta.getDepth() } );
				viewSetups.add( new ViewSetup( viewSetups.size(), null, dim, voxelSize, c, a, i ) );
			}

	return viewSetups;
}
 

@Test
public void testImageFactory() {
	final Dimensions dim = new FinalDimensions(10, 10, 10);

	@SuppressWarnings("unchecked")
	final Img<DoubleType> arrayImg = (Img<DoubleType>) ops.run(
		CreateImgFromDimsAndType.class, dim, new DoubleType(),
		new ArrayImgFactory<DoubleType>());
	final Class<?> arrayFactoryClass = arrayImg.factory().getClass();
	assertEquals("Image Factory: ", ArrayImgFactory.class, arrayFactoryClass);

	@SuppressWarnings("unchecked")
	final Img<DoubleType> cellImg = (Img<DoubleType>) ops.run(
		CreateImgFromDimsAndType.class, dim, new DoubleType(),
		new CellImgFactory<DoubleType>());
	final Class<?> cellFactoryClass = cellImg.factory().getClass();
	assertEquals("Image Factory: ", CellImgFactory.class, cellFactoryClass);
}
 

@Test
@SuppressWarnings("unchecked")
public void testPadShiftKernel() {
	long[] dims = new long[] { 1024, 1024 };
	Img<ComplexDoubleType> test = ops.create().img(new FinalDimensions(dims),
		new ComplexDoubleType());

	RandomAccessibleInterval<ComplexDoubleType> shift =
		(RandomAccessibleInterval<ComplexDoubleType>) ops.run(
			PadShiftKernel.class, test, new FinalDimensions(dims));

	RandomAccessibleInterval<ComplexDoubleType> shift2 =
		(RandomAccessibleInterval<ComplexDoubleType>) ops.run(
			PadShiftKernelFFTMethods.class, test, new FinalDimensions(dims));

	// assert there was no additional padding done by PadShiftKernel
	assertEquals(1024, shift.dimension(0));
	// assert that PadShiftKernelFFTMethods padded to the FFTMethods fast size
	assertEquals(1120, shift2.dimension(0));

}
 

@Override
public RandomAccessibleInterval<C> calculate(
	final RandomAccessibleInterval<T> input)
{
	// calculate the padded size
	long[] paddedSize = new long[in().numDimensions()];

	for (int d = 0; d < in().numDimensions(); d++) {
		paddedSize[d] = in().dimension(d);

		if (borderSize != null) {
			paddedSize[d] += borderSize[d];
		}
	}

	Dimensions paddedDimensions = new FinalDimensions(paddedSize);

	// create the complex output
	RandomAccessibleInterval<C> output = createOp.calculate(paddedDimensions);

	// pad the input
	RandomAccessibleInterval<T> paddedInput = padOp.calculate(input,
		paddedDimensions);

	// compute and return fft
	fftMethodsOp.compute(paddedInput, output);

	return output;

}
 

/**
 * Test the coordinate op version of the equation using 4 dimensions
 */
@Test
public void testEquation4DOp() {

	final long[] size4D = new long[] { 5, 5, 5, 5 };

	final Dimensions dimensions4D = new FinalDimensions(size4D);

	final Img<ShortType> image = ops.create().img(dimensions4D,
		new ShortType());

	// implement c[0]+10*c[1]+100*c[3]+1000*c[4]
	final UnaryFunctionOp<long[], Double> op =
		new AbstractUnaryFunctionOp<long[], Double>()
		{

			@Override
			public Double calculate(final long[] coords) {
				final double result = coords[0] + 10 * coords[1] + 100 * coords[2] +
					1000 * coords[3];

				return result;
			}
		};

	ops.run(DefaultCoordinatesEquation.class, image, op);

	final RandomAccess<ShortType> ra = image.randomAccess();

	ra.setPosition(new long[] { 1, 2, 2, 3 });

	assertEquals(1 + 10 * 2 + 100 * 2 + 1000 * 3, ra.get().getRealFloat(),
		0.000001);

}
 

/**
 * Validate generation of a 3D diffraction kernel against a comparable result
 * from <a href="http://bigwww.epfl.ch/algorithms/psfgenerator/">
 * PSFGenerator</a>.
 * <p>
 * It is worth noting that results are only comparable between the two when
 * using a particle position relative to the coverslip of 0. This is because
 * imagej-ops automatically crops and centers kernels produced by the Fast
 * Gibson-Lanni implementation in {@link DefaultCreateKernelGibsonLanni} while
 * the Gibson &amp; Lanni kernels produced by PSFGenerator are not. See this
 * github issue
 * <a href="https://github.com/imagej/imagej-ops/issues/550">thread</a> for
 * more details.
 * </p>
 * <p>
 * It is also worth noting that when using a particle position of 0, the model
 * degenerates to a standard Born &amp; Wolf PSF model [1].
 * </p>
 * <h3>References:</h3>
 * <ol>
 * <li>Jizhou Li, Feng Xue, and Thierry Blu, "Fast and accurate
 * three-dimensional point spread function computation for fluorescence
 * microscopy," J. Opt. Soc. Am. A 34, 1029-1034 (2017)</li>
 * </ol>
 */
@Test
public void testKernelDiffraction3D() {
	final Dimensions dims = new FinalDimensions(16, 16, 8);
	final double NA = 1.4; // numerical aperture
	final double lambda = 610E-09; // wavelength
	final double ns = 1.33; // specimen refractive index
	final double ni = 1.5; // immersion refractive index, experimental
	final double resLateral = 100E-9; // lateral pixel size
	final double resAxial = 250E-9; // axial pixel size
	final DoubleType type = new DoubleType(); // pixel type of created kernel
	
	// NB: It is important that this remain 0 for comparison to PSFGenerator.
	final double pZ = 0D; // position of particle

	final Img<DoubleType> kernel = 
		ops.create().kernelDiffraction(dims, NA, lambda, ns, ni, resLateral,
			resAxial, pZ, type);
	
	// The following image used for comparison was generated via PSFGenerator
	// with these arguments (where any not mentioned were left at default):
	// Optical Model: "Gibson & Lanni 3D Optical Model"
	// Particle position Z: 0
	// Wavelength: 610 (nm)
	// Pixelsize XY: 100 (nm)
	// Z-step: 250 (nm)
	// Size XYZ: 16, 16, 8
	Img<DoubleType> expected = openDoubleImg("kern3d1.tif");
	
	assertArrayEquals(asArray(expected), asArray(kernel), 1e-4);
}
 

@Test
public void testImageType() {
	final Dimensions dim = new FinalDimensions(10, 10, 10);

	assertEquals("Image Type: ", BitType.class, ((Img<?>) ops.run(
		CreateImgFromDimsAndType.class, dim, new BitType())).firstElement()
			.getClass());

	assertEquals("Image Type: ", ByteType.class, ((Img<?>) ops.run(
		CreateImgFromDimsAndType.class, dim, new ByteType())).firstElement()
			.getClass());

	assertEquals("Image Type: ", UnsignedByteType.class, ((Img<?>) ops.create()
		.img(dim, new UnsignedByteType())).firstElement().getClass());

	assertEquals("Image Type: ", IntType.class, ((Img<?>) ops.run(
		CreateImgFromDimsAndType.class, dim, new IntType())).firstElement()
			.getClass());

	assertEquals("Image Type: ", FloatType.class, ((Img<?>) ops.run(
		CreateImgFromDimsAndType.class, dim, new FloatType())).firstElement()
			.getClass());

	assertEquals("Image Type: ", DoubleType.class, ((Img<?>) ops.run(
		CreateImgFromDimsAndType.class, dim, new DoubleType())).firstElement()
			.getClass());
}
 

@Test
public void testImgFromImg() {
	// create img
	final Img<ByteType> img =
		ops.create().img(new FinalDimensions(1), new ByteType());
	@SuppressWarnings("unchecked")
	final Img<ByteType> newImg = (Img<ByteType>) ops.run(CreateImgFromImg.class,
		img);

	// should both be ByteType. New Img shouldn't be DoubleType (default)
	assertEquals(img.firstElement().getClass(), newImg.firstElement()
		.getClass());
}
 

@Before
public void createImages() {
	final FinalDimensions dims = FinalDimensions.wrap(new long[] {10, 10});
	in = ops.create().img(dims, new ShortType());
	addNoise(in);
	out = ops.create().img(dims, new ByteType());
}
 

@Override
@SuppressWarnings("unchecked")
public O calculate(final I input, final Dimensions paddedDimensions) {

	final long[] paddedSize = new long[paddedDimensions.numDimensions()];
	paddedDimensions.dimensions(paddedSize);

	O inputInterval = (O) FFTMethods.paddingIntervalCentered(input,
		FinalDimensions.wrap(paddedSize));

	return inputInterval;
}
 

/**
 * Calculates complex FFT size for real to complex FFT
 * 
 * @param fast if true calculate size for fast FFT
 * @param inputDimensions original real dimensions
 * @return complex FFT dimensions
 */
public static Dimensions getFFTDimensionsRealToComplex(final boolean fast,
	final Dimensions inputDimensions)
{
	final long[] paddedSize = new long[inputDimensions.numDimensions()];
	final long[] fftSize = new long[inputDimensions.numDimensions()];

	dimensionsRealToComplex(fast, inputDimensions, paddedSize, fftSize);

	return new FinalDimensions(fftSize);

}
 

/**
 * Calculates padding size size for real to complex FFT
 * 
 * @param fast if true calculate size for fast FFT
 * @param inputDimensions original real dimensions
 * @return padded real dimensions
 */
public static Dimensions getPaddedInputDimensionsRealToComplex(
	final boolean fast, final Dimensions inputDimensions)
{
	final long[] paddedSize = new long[inputDimensions.numDimensions()];
	final long[] fftSize = new long[inputDimensions.numDimensions()];

	dimensionsRealToComplex(fast, inputDimensions, paddedSize, fftSize);

	return new FinalDimensions(paddedSize);

}
 

/**
 * Updates the cached imageMetaData
 */
protected void updateMetaDataCache( final ViewId viewId,
		final int w, final int h, final int d,
		final double calX, final double calY, final double calZ )
{
	imageMetaDataCache.put( viewId, new ValuePair< Dimensions, VoxelDimensions >(
			new FinalDimensions( new long[] { w, h, d } ),
			new FinalVoxelDimensions( "", calX, calY, calZ ) ) );

	// links the viewSetupId to the last added viewId, overwrites earlier entries
	viewIdLookUp.put( viewId.getViewSetupId(), viewId );
}
 

public static void main( String[] args )
	{
		final Point p = new Point( 90, 90 ); // identical to (-10,-10), so subpixel localization can move on periodic condition outofbounds
		PhaseCorrelationPeak2 pcp = new PhaseCorrelationPeak2( p, 5 );

		Dimensions pcmDims = new FinalDimensions( 100, 100 );
		Dimensions p1 = new FinalDimensions( 80, 81 );
		Dimensions p2 = new FinalDimensions( 91, 90 );

		final List<PhaseCorrelationPeak2> peaks = expandPeakToPossibleShifts( pcp, pcmDims, p1, p2 );

		for ( final PhaseCorrelationPeak2 pc : peaks )
			System.out.println( Util.printCoordinates( pc.getShift() ) );

		Img< FloatType > a = ImgLib2Util.openAs32Bit( new File( "73.tif.zip" ) );
		Img< FloatType > b = ImgLib2Util.openAs32Bit( new File( "74.tif.zip" ) );

//		BenchmarkHelper.benchmarkAndPrint( 10, true, new Runnable()
//		{
//			@Override
//			public void run()
//			{
//				System.out.println( getCorrelation ( a, b ) );
//			}
//		} );

		BenchmarkHelper.benchmarkAndPrint( 10, true, new Runnable()
		{
			@Override
			public void run()
			{
				PairwiseStitching.getShift( a, b, new Translation3D(), new Translation3D(),
						new PairwiseStitchingParameters(), Executors.newFixedThreadPool( Runtime.getRuntime().availableProcessors() ) );
			}
		} );

//		System.out.println( getCorrelation ( a, b ) );
//		System.out.println( getCorrelation ( a, c ) );
//		System.out.println( getCorrelation ( b, c ) );
	}
 

/**
 * calculate the size of an extended image big enough to hold dim1 and dim2
 * with each dimension also enlarged by extension pixels on each side (but at most by the original image size)
 * @param dim1 first Dimensions
 * @param dim2 second Dimensions
 * @param extension: number of pixels to add at each side in each dimension
 * @return extended dimensions
 */
public static FinalDimensions getExtendedSize(Dimensions dim1, Dimensions dim2, int [] extension) {
	long[] extDims = new long[dim1.numDimensions()];
	for (int i = 0; i <dim1.numDimensions(); i++){
		extDims[i] = dim1.dimension(i) > dim2.dimension(i) ? dim1.dimension(i) : dim2.dimension(i);
		long extBothSides = extDims[i] < extension[i] ? extDims[i] * 2 : extension[i] * 2;
		extDims[i] += extBothSides;
	}
	return new FinalDimensions(extDims);		
}
 

/**
 * Duplicates all Angles and Illuminations that are processed.
 * The size of the List> ViewSetup < is therefore equal to
 * the number of channels * number of processed angles * number of processed illuminations
 * 
 * @param spimData
 * @param viewIdsToProcess
 * @param bb
 * @return
 */
public static Map< ViewSetup, ViewSetup > assembleNewViewSetupsSequential(
		final SpimData2 spimData,
		final List< ViewId > viewIdsToProcess,
		final BoundingBoxGUI bb )
{
	final HashMap< ViewSetup, ViewSetup > map = new HashMap< ViewSetup, ViewSetup >();

	// make many new view setups with new angles&illuminations
	int maxViewSetupIndex = -1;
	int maxAngleIndex = -1;
	int maxIllumIndex = -1;
	
	for ( final ViewSetup v : spimData.getSequenceDescription().getViewSetups().values() )
		maxViewSetupIndex = Math.max( maxViewSetupIndex, v.getId() );

	for ( final Angle a : spimData.getSequenceDescription().getAllAngles().values() )
		maxAngleIndex = Math.max( maxAngleIndex, a.getId() );

	for ( final Illumination i : spimData.getSequenceDescription().getAllIlluminations().values() )
		maxIllumIndex = Math.max( maxIllumIndex, i.getId() );

	final String unit = spimData.getSequenceDescription().getViewSetupsOrdered().get( 0 ).getVoxelSize().unit();
	
	// get the minimal resolution of all calibrations relative to the downsampling
	final double minResolution = Apply_Transformation.assembleAllMetaData(
			spimData.getSequenceDescription(),
			spimData.getSequenceDescription().getViewDescriptions().values() ) * bb.getDownSampling();

	// every combination of old angle and old illumination of each channel gets a new viewsetup
	final List< Angle > oldAngles = SpimData2.getAllAnglesSorted( spimData, viewIdsToProcess );
	final List< Illumination > oldIllums = SpimData2.getAllIlluminationsSorted( spimData, viewIdsToProcess );

	final HashMap< Angle, Angle > mapOldToNewAngles = new HashMap< Angle, Angle >();
	final HashMap< Illumination, Illumination > mapOldToNewIlluminations = new HashMap< Illumination, Illumination >();

	//final List< Pair< Angle, Angle > > mapOldToNewAngles = new ArrayList< Pai r< Angle, Angle > >();
	//final List< Pair< Illumination, Illumination > > mapOldToNewIlluminations = new ArrayList< Pair< Illumination, Illumination > >();

	for ( int i = 0; i < oldAngles.size(); ++i )
		mapOldToNewAngles.put(
			oldAngles.get( i ),
			new Angle(
				maxAngleIndex + i + 1,
				"Transf_" + oldAngles.get( i ).getName() + "_" + maxAngleIndex + i + 1,
				oldAngles.get( i ).getRotationAngleDegrees(),
				oldAngles.get( i ).getRotationAxis() ) );

	for ( int i = 0; i < oldIllums.size(); ++i )
		mapOldToNewIlluminations.put(
			oldIllums.get( i ),
			new Illumination(
				maxIllumIndex + i + 1,
				"Transf_" + oldIllums.get( i ).getName() + "_" + maxIllumIndex + i + 1 ) );

	// now every viewsetup of every viewdescription that is fused (maybe for several timepoints) gets a new viewsetup
	for ( final ViewId viewId : viewIdsToProcess )
	{
		final ViewDescription oldVD = spimData.getSequenceDescription().getViewDescription( viewId );

		if ( oldVD.isPresent() )
		{
			final ViewSetup oldSetup = oldVD.getViewSetup();

			// no new viewsetup defined for this old viewsetup
			if ( !map.containsKey( oldSetup ) )
			{
				// get the new angle & illumination object
				final Channel channel = oldSetup.getChannel();
				final Angle oldAngle = oldSetup.getAngle();
				final Illumination oldIllumination = oldSetup.getIllumination();

				// make the new viewsetup
				final Angle newAngle = mapOldToNewAngles.get( oldAngle );
				final Illumination newIllumination = mapOldToNewIlluminations.get( oldIllumination );

				final ViewSetup newSetup = new ViewSetup( 
						++maxViewSetupIndex,
						null,
						new FinalDimensions( bb.getDimensions() ), 
						new FinalVoxelDimensions ( unit, new double[]{ minResolution, minResolution, minResolution } ),
						channel,
						newAngle,
						newIllumination );

				map.put( oldSetup, newSetup );
			}
		}
	}

	return map;
}
 

/**
 * Creates the List of {@link ViewSetup} for the {@link SpimData} object.
 * The {@link ViewSetup} are defined independent of the {@link TimePoint},
 * each {@link TimePoint} should have the same {@link ViewSetup}s. The {@link MissingViews}
 * class defines if some of them are missing for some of the {@link TimePoint}s
 *
 * @return
 */
protected ArrayList<ViewSetup> createViewSetups(final SlideBook6MetaData meta) {
    // TODO: query rotation angle of each SlideBook channel
    final double[] yaxis = new double[]{0, 1, 0};
    final ArrayList<Angle> angles = new ArrayList<Angle>();
    final Angle angleA = new Angle(0, "Path_A");
    angleA.setRotation(yaxis, defaultAngles[0]);
    angles.add(angleA);

    final Angle angleB = new Angle(1, "Path_B");
    angleB.setRotation(yaxis, defaultAngles[1]);
    angles.add(angleB);

    // define multiple illuminations, one for every capture in the slide file
    final ArrayList<ViewSetup> viewSetups = new ArrayList<ViewSetup>();

    int firstCapture = defaultCapture;
    int numCaptures = 1;
    if (defaultCapture == -1) {
        firstCapture = 0;
        numCaptures = meta.numCaptures();
    }

    // create one illumination for each capture if defaultCapture == -1, otherwise just one capture
    for (int capture = firstCapture; capture < firstCapture + numCaptures; capture++) {
        final int channels = meta.numChannels(capture);

        // multi-angle captures must have pairs of channels
        if (channels > 1 && channels % 2 == 0) {
            String imageName = meta.imageName(capture);
            imageName = imageName.replaceAll("[^a-zA-Z0-9_-]", "_"); // convert illegal characters to _
            imageName = imageName.toLowerCase(); // convert to lowercase
            // up to 8 illuminations (SlideBook channels) per SlideBook capture
            final Illumination i = new Illumination(capture * 8, imageName);

            for (int ch = 0; ch < channels / 2; ch++) {
                // use name of first channel, SlideBook channels are diSPIM angles and each SlideBook image should have two angles per channel
                String channelName = meta.channels(capture)[ch * 2];
                channelName = channelName.replaceAll("[^a-zA-Z0-9_-]", "_"); // convert illegal characters to _
                channelName = channelName.toLowerCase(); // convert to lowercase

                final Channel channel = new Channel(ch, channelName);

                for (final Angle a : angles) {
                    float voxelSizeUm = defaultCalibrations[0];
                    float zSpacing = defaultCalibrations[2];

                    final VoxelDimensions voxelSize = new FinalVoxelDimensions("um", voxelSizeUm, voxelSizeUm, zSpacing);
                    final Dimensions dim = new FinalDimensions(new long[]{meta.imageSize(capture)[0], meta.imageSize(capture)[1], meta.imageSize(capture)[2]});

                    viewSetups.add(new ViewSetup(viewSetups.size(), a.getName(), dim, voxelSize, channel, a, i));
                }
            }
        }
    }

    return viewSetups;
}
 

/**
 * 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);
	}
}
 

/**
 * Creates one new Angle and one new Illumination for the fused dataset.
 * The size of the List> ViewSetup < is therefore equal to the number of channels
 * 
 * @param spimData
 * @param viewIdsToProcess
 * @param bb
 * @param newAngleName
 * @param newIlluminationName
 * @return
 */
public static Map< ViewSetup, ViewSetup > assembleNewViewSetupsFusion(
		final SpimData2 spimData,
		final List< ViewId > viewIdsToProcess,
		final BoundingBoxGUI bb,
		final String newAngleName,
		final String newIlluminationName )
{
	final HashMap< ViewSetup, ViewSetup > map = new HashMap< ViewSetup, ViewSetup >();

	int maxViewSetupIndex = -1;
	int maxAngleIndex = -1;
	int maxIllumIndex = -1;

	// make a new viewsetup
	for ( final ViewSetup v : spimData.getSequenceDescription().getViewSetups().values() )
		maxViewSetupIndex = Math.max( maxViewSetupIndex, v.getId() );

	// that has a new angle
	for ( final Angle a : spimData.getSequenceDescription().getAllAngles().values() )
		maxAngleIndex = Math.max( maxAngleIndex, a.getId() );

	// and a new illumination
	for ( final Illumination i : spimData.getSequenceDescription().getAllIlluminations().values() )
		maxIllumIndex = Math.max( maxIllumIndex, i.getId() );

	final Angle newAngle = new Angle( maxAngleIndex + 1, newAngleName + "_" + ( maxAngleIndex + 1 ) );
	final Illumination newIllum = new Illumination( maxIllumIndex + 1, newIlluminationName + "_" + ( maxIllumIndex + 1 ) );
	
	final String unit = spimData.getSequenceDescription().getViewSetupsOrdered().get( 0 ).getVoxelSize().unit();
	
	// get the minimal resolution of all calibrations relative to the downsampling
	final double minResolution = Apply_Transformation.assembleAllMetaData(
			spimData.getSequenceDescription(),
			spimData.getSequenceDescription().getViewDescriptions().values() ) * bb.getDownSampling();

	// one new new viewsetup for every channel that is fused (where fused means combining one or more angles&illuminations into one new image)
	for ( final Channel channel : SpimData2.getAllChannelsSorted( spimData, viewIdsToProcess ) )
	{
		final ViewSetup newSetup = new ViewSetup(
				++maxViewSetupIndex,
				null,
				new FinalDimensions( bb.getDimensions() ),
				new FinalVoxelDimensions ( unit, new double[]{ minResolution, minResolution, minResolution } ),
				channel,
				newAngle,
				newIllum );
		
		// all viewsetups of the current channel that are fused point to the same new viewsetup
		for ( final ViewId viewId : viewIdsToProcess )
		{
			final ViewDescription oldVD = spimData.getSequenceDescription().getViewDescription( viewId );
			final ViewSetup oldSetup = oldVD.getViewSetup();

			// is this old setup from the current channel? then point to it
			if ( oldVD.isPresent() && oldSetup.getChannel().getId() == channel.getId() )
				map.put( oldSetup, newSetup );
		}
	}

	return map;
}
 

public static SpimData twoAngles()
{
	final ArrayList< ViewSetup > setups = new ArrayList< ViewSetup >();
	final ArrayList< ViewRegistration > registrations = new ArrayList< ViewRegistration >();

	final Channel c0 = new Channel( 0, "test" );
	final Angle a0 = new Angle( 0 );
	final Angle a1 = new Angle( 1 );
	final Illumination i0 = new Illumination( 0 );

	final Dimensions d0 = new FinalDimensions( 512l, 512l, 86l );
	final VoxelDimensions vd0 = new FinalVoxelDimensions( "px", 0.4566360, 0.4566360, 2.0000000 );

	setups.add( new ViewSetup( 0, "setup 0", d0, vd0, c0, a0, i0 ) );
	setups.add( new ViewSetup( 1, "setup 1", d0, vd0, c0, a1, i0 ) );

	final ArrayList< TimePoint > t = new ArrayList< TimePoint >();
	t.add( new TimePoint( 0 ) );
	final TimePoints timepoints = new TimePoints( t );

	final ArrayList< ViewId > missing = new ArrayList< ViewId >();
	final MissingViews missingViews = new MissingViews( missing );

	final ImgLoader imgLoader = new ImgLoader()
	{
		@Override
		public SetupImgLoader< ? > getSetupImgLoader( int setupId )
		{
			return new MySetupImgLoader( setupId );
		}
	};

	for ( final ViewSetup vs : setups )
	{
		final ViewRegistration vr = new ViewRegistration( t.get( 0 ).getId(), vs.getId() );

		final double minResolution = Math.min( Math.min( vs.getVoxelSize().dimension( 0 ), vs.getVoxelSize().dimension( 1 ) ), vs.getVoxelSize().dimension( 2 ) );
		
		final double calX = vs.getVoxelSize().dimension( 0 ) / minResolution;
		final double calY = vs.getVoxelSize().dimension( 1 ) / minResolution;
		final double calZ = vs.getVoxelSize().dimension( 2 ) / minResolution;
		
		final AffineTransform3D m = new AffineTransform3D();
		m.set( calX, 0.0f, 0.0f, 0.0f, 
			   0.0f, calY, 0.0f, 0.0f,
			   0.0f, 0.0f, calZ, 0.0f );
		final ViewTransform vt = new ViewTransformAffine( "Calibration", m );
		vr.preconcatenateTransform( vt );

		vr.updateModel();		
		
		registrations.add( vr );
	}

	final SequenceDescription sd = new SequenceDescription( timepoints, setups, imgLoader, missingViews );
	final SpimData data = new SpimData( new File( "" ), sd, new ViewRegistrations( registrations ) );

	return data;
}
 

@Test
public void testKernelDiffraction2D() {
	final Dimensions dims = new FinalDimensions(10, 10);
	final double NA = 1.4; // numerical aperture
	final double lambda = 610E-09; // wavelength
	final double ns = 1.33; // specimen refractive index
	final double ni = 1.5; // immersion refractive index, experimental
	final double resLateral = 100E-9; // lateral pixel size
	final double resAxial = 250E-9; // axial pixel size
	final double pZ = 2000E-9D; // position of particle
	final DoubleType type = new DoubleType(); // pixel type of created kernel

	final Img<DoubleType> kernel = //
		ops.create().kernelDiffraction(dims, NA, lambda, ns, ni, resLateral,
			resAxial, pZ, type);

	final double[] expected = { 0.03298495871588273, 0.04246786111102021,
		0.0543588031627261, 0.06650574371357207, 0.07370280610722534,
		0.07370280610722534, 0.06650574371357207, 0.0543588031627261,
		0.04246786111102021, 0.03298495871588273, 0.04246786111102021,
		0.05962205221267819, 0.08320071670150801, 0.10800022978800021,
		0.1247473245002288, 0.1247473245002288, 0.10800022978800021,
		0.08320071670150801, 0.05962205221267819, 0.04246786111102021,
		0.0543588031627261, 0.08320071670150801, 0.1247473245002288,
		0.1971468112729564, 0.2691722397359577, 0.2691722397359577,
		0.1971468112729564, 0.1247473245002288, 0.08320071670150801,
		0.0543588031627261, 0.06650574371357207, 0.10800022978800021,
		0.1971468112729564, 0.40090474481128285, 0.6227157103102976,
		0.6227157103102976, 0.40090474481128285, 0.1971468112729564,
		0.10800022978800021, 0.06650574371357207, 0.07370280610722534,
		0.1247473245002288, 0.2691722397359577, 0.6227157103102976, 1.0, 1.0,
		0.6227157103102976, 0.2691722397359577, 0.1247473245002288,
		0.07370280610722534, 0.07370280610722534, 0.1247473245002288,
		0.2691722397359577, 0.6227157103102976, 1.0, 1.0, 0.6227157103102976,
		0.2691722397359577, 0.1247473245002288, 0.07370280610722534,
		0.06650574371357207, 0.10800022978800021, 0.1971468112729564,
		0.40090474481128285, 0.6227157103102976, 0.6227157103102976,
		0.40090474481128285, 0.1971468112729564, 0.10800022978800021,
		0.06650574371357207, 0.0543588031627261, 0.08320071670150801,
		0.1247473245002288, 0.1971468112729564, 0.2691722397359577,
		0.2691722397359577, 0.1971468112729564, 0.1247473245002288,
		0.08320071670150801, 0.0543588031627261, 0.04246786111102021,
		0.05962205221267819, 0.08320071670150801, 0.10800022978800021,
		0.1247473245002288, 0.1247473245002288, 0.10800022978800021,
		0.08320071670150801, 0.05962205221267819, 0.04246786111102021,
		0.03298495871588273, 0.04246786111102021, 0.0543588031627261,
		0.06650574371357207, 0.07370280610722534, 0.07370280610722534,
		0.06650574371357207, 0.0543588031627261, 0.04246786111102021,
		0.03298495871588273,

	};

	assertArrayEquals(expected, asArray(kernel), 1e-4);
}
 

/**
 * create FFT memory, create FFT filter and run it
 */
@Override
public void compute(
	final RandomAccessibleInterval<I> input,
	final RandomAccessibleInterval<K> kernel,
	final RandomAccessibleInterval<O> output)
{

	//RandomAccessibleInterval<O> output = createOutput(input, kernel);

	final int numDimensions = input.numDimensions();

	// 1. Calculate desired extended size of the image

	final long[] paddedSize = new long[numDimensions];

	if (borderSize == null) {
		// if no border size was passed in, then extend based on kernel size
		for (int d = 0; d < numDimensions; ++d) {
			paddedSize[d] = (int) input.dimension(d) + (int) kernel.dimension(d) -
				1;
		}

	}
	else {
		// if borderSize was passed in
		for (int d = 0; d < numDimensions; ++d) {

			paddedSize[d] = Math.max(kernel.dimension(d) + 2 * borderSize[d], input
				.dimension(d) + 2 * borderSize[d]);
		}
	}

	RandomAccessibleInterval<I> paddedInput = padOp.calculate(input,
		new FinalDimensions(paddedSize));

	RandomAccessibleInterval<K> paddedKernel = padKernelOp.calculate(kernel,
		new FinalDimensions(paddedSize));

	RandomAccessibleInterval<C> fftInput = createOp.calculate(
		new FinalDimensions(paddedSize));

	RandomAccessibleInterval<C> fftKernel = createOp.calculate(
		new FinalDimensions(paddedSize));

	// TODO: in this case it is difficult to match the filter op in the
	// 'initialize' as we don't know the size yet, thus we can't create
	// memory
	// for the FFTs
	filter = createFilterComputer(paddedInput, paddedKernel, fftInput,
		fftKernel, output);

	filter.compute(paddedInput, paddedKernel, output);
}
 

@Override
public void compute(final RandomAccessibleInterval<T> input,
	final IterableInterval<DoubleType> tubeness)
{
	cancelReason = null;

	final int numDimensions = input.numDimensions();
	// Sigmas in pixel units.
	final double[] sigmas = new double[numDimensions];
	for (int d = 0; d < sigmas.length; d++) {
		final double cal = d < calibration.length ? calibration[d] : 1;
		sigmas[d] = sigma / cal;
	}

	/*
	 * Hessian.
	 */

	// Get a suitable image factory.
	final long[] gradientDims = new long[numDimensions + 1];
	final long[] hessianDims = new long[numDimensions + 1];
	for (int d = 0; d < numDimensions; d++) {
		hessianDims[d] = input.dimension(d);
		gradientDims[d] = input.dimension(d);
	}
	hessianDims[numDimensions] = numDimensions * (numDimensions + 1) / 2;
	gradientDims[numDimensions] = numDimensions;
	final Dimensions hessianDimensions = FinalDimensions.wrap(hessianDims);
	final FinalDimensions gradientDimensions = FinalDimensions.wrap(
		gradientDims);
	final ImgFactory<DoubleType> factory = ops().create().imgFactory(
		hessianDimensions);
	final Img<DoubleType> hessian = factory.create(hessianDimensions,
		new DoubleType());
	final Img<DoubleType> gradient = factory.create(gradientDimensions,
		new DoubleType());
	final Img<DoubleType> gaussian = factory.create(input, new DoubleType());

	// Handle multithreading.
	final int nThreads = Runtime.getRuntime().availableProcessors();
	final ExecutorService es = threadService.getExecutorService();

	try {
		// Hessian calculation.
		HessianMatrix.calculateMatrix(Views.extendBorder(input), gaussian,
			gradient, hessian, new OutOfBoundsBorderFactory<>(), nThreads, es,
			sigma);

		statusService.showProgress(1, 3);
		if (isCanceled()) return;

		// Hessian eigenvalues.
		final RandomAccessibleInterval<DoubleType> evs = TensorEigenValues
			.calculateEigenValuesSymmetric(hessian, TensorEigenValues
				.createAppropriateResultImg(hessian, factory, new DoubleType()),
				nThreads, es);

		statusService.showProgress(2, 3);
		if (isCanceled()) return;

		final AbstractUnaryComputerOp<Iterable<DoubleType>, DoubleType> method;
		switch (numDimensions) {
			case 2:
				method = new Tubeness2D(sigma);
				break;
			case 3:
				method = new Tubeness3D(sigma);
				break;
			default:
				System.err.println("Cannot compute tubeness for " + numDimensions +
					"D images.");
				return;
		}
		ops().transform().project(tubeness, evs, method, numDimensions);

		statusService.showProgress(3, 3);

		return;
	}
	catch (final IncompatibleTypeException | InterruptedException
			| ExecutionException e)
	{
		e.printStackTrace();
		return;
	}
}
 

protected void createNormalizationImageSemiNonCirculant(Interval fastFFTInterval) {

		// k is the window size (valid image region)
		final int length = k.numDimensions();

		final long[] n = new long[length];
		final long[] nFFT = new long[length];

		// n is the valid image size plus the extended region
		// also referred to as object space size
		for (int d = 0; d < length; d++) {
			n[d] = k.dimension(d) + l.dimension(d) - 1;
		}

		// nFFT is the size of n after (potentially) extending further
		// to a fast FFT size
		for (int d = 0; d < length; d++) {
			nFFT[d] = fastFFTInterval.dimension(d);
		}

		FinalDimensions fd = new FinalDimensions(nFFT);

		// create the normalization image
		normalization = create.calculate(fd);

		// size of the measurement window
		final Point size = new Point(length);
		final long[] sizel = new long[length];

		for (int d = 0; d < length; d++) {
			size.setPosition(k.dimension(d), d);
			sizel[d] = k.dimension(d);
		}

		// starting point of the measurement window when it is centered in fft space
		final Point start = new Point(length);
		final long[] startl = new long[length];
		final long[] endl = new long[length];

		for (int d = 0; d < length; d++) {
			start.setPosition((nFFT[d] - k.dimension(d)) / 2, d);
			startl[d] = (nFFT[d] - k.dimension(d)) / 2;
			endl[d] = startl[d] + sizel[d] - 1;
		}

		// size of the object space
		final Point maskSize = new Point(length);
		final long[] maskSizel = new long[length];

		for (int d = 0; d < length; d++) {
			maskSize.setPosition(Math.min(n[d], nFFT[d]), d);
			maskSizel[d] = Math.min(n[d], nFFT[d]);
		}

		// starting point of the object space within the fft space
		final Point maskStart = new Point(length);
		final long[] maskStartl = new long[length];

		for (int d = 0; d < length; d++) {
			maskStart.setPosition((Math.max(0, nFFT[d] - n[d]) / 2), d);
			maskStartl[d] = (Math.max(0, nFFT[d] - n[d]) / 2);
		}

		final RandomAccessibleInterval<O> temp = Views.interval(normalization,
			new FinalInterval(startl, endl));
		final Cursor<O> normCursor = Views.iterable(temp).cursor();

		// draw a cube the size of the measurement space
		while (normCursor.hasNext()) {
			normCursor.fwd();
			normCursor.get().setReal(1.0);
		}

		final Img<O> tempImg = create.calculate(fd);

		// 3. correlate psf with the output of step 2.
		correlater.compute(normalization, tempImg);

		normalization = tempImg;

		final Cursor<O> cursorN = normalization.cursor();

		while (cursorN.hasNext()) {
			cursorN.fwd();

			if (cursorN.get().getRealFloat() <= 1e-3f) {
				cursorN.get().setReal(1.0f);

			}
		}
	}
 

public static SpimData grid3x2()
{
	final ArrayList< ViewSetup > setups = new ArrayList< ViewSetup >();
	final ArrayList< ViewRegistration > registrations = new ArrayList< ViewRegistration >();

	final Channel c0 = new Channel( 0, "RFP" );
	final Channel c1 = new Channel( 1, "YFP" );
	final Channel c2 = new Channel( 2, "GFP" );

	final Angle a0 = new Angle( 0 );
	final Illumination i0 = new Illumination( 0 );

	final Tile t0 = new Tile( 0, "Tile0", new double[]{ 0.0, 0.0, 0.0 } );
	final Tile t1 = new Tile( 1, "Tile1", new double[]{ 450.0, 0.0, 0.0 } );
	final Tile t2 = new Tile( 2, "Tile2", new double[]{ 0.0, 450.0, 0.0 } );
	final Tile t3 = new Tile( 3, "Tile3", new double[]{ 450.0, 450.0, 0.0 } );

	final Dimensions d0 = new FinalDimensions( 512l, 512l, 86l );
	final VoxelDimensions vd0 = new FinalVoxelDimensions( "px", 0.4566360, 0.4566360, 2.0000000 );

	setups.add( new ViewSetup( 0, "setup 0", d0, vd0, t0, c0, a0, i0 ) );
	setups.add( new ViewSetup( 1, "setup 1", d0, vd0, t1, c0, a0, i0 ) );
	setups.add( new ViewSetup( 2, "setup 2", d0, vd0, t2, c0, a0, i0 ) );
	setups.add( new ViewSetup( 3, "setup 3", d0, vd0, t3, c0, a0, i0 ) );

	setups.add( new ViewSetup( 4, "setup 4", d0, vd0, t0, c1, a0, i0 ) );
	setups.add( new ViewSetup( 5, "setup 5", d0, vd0, t1, c1, a0, i0 ) );
	setups.add( new ViewSetup( 6, "setup 6", d0, vd0, t2, c1, a0, i0 ) );
	setups.add( new ViewSetup( 7, "setup 7", d0, vd0, t3, c1, a0, i0 ) );

	setups.add( new ViewSetup( 8, "setup 8", d0, vd0, t0, c2, a0, i0 ) );
	setups.add( new ViewSetup( 9, "setup 9", d0, vd0, t1, c2, a0, i0 ) );
	setups.add( new ViewSetup( 10, "setup 10", d0, vd0, t2, c2, a0, i0 ) );
	setups.add( new ViewSetup( 11, "setup 11", d0, vd0, t3, c2, a0, i0 ) );

	final ArrayList< TimePoint > t = new ArrayList< TimePoint >();
	t.add( new TimePoint( 0 ) );
	final TimePoints timepoints = new TimePoints( t );

	final ArrayList< ViewId > missing = new ArrayList< ViewId >();
	final MissingViews missingViews = new MissingViews( missing );

	final ImgLoader imgLoader = new ImgLoader()
	{
		@Override
		public SetupImgLoader< ? > getSetupImgLoader( int setupId )
		{
			return new MySetupImgLoader( setupId );
		}
	};

	for ( final ViewSetup vs : setups )
	{
		final ViewRegistration vr = new ViewRegistration( t.get( 0 ).getId(), vs.getId() );

		final Tile tile = vs.getTile();

		final double minResolution = Math.min( Math.min( vs.getVoxelSize().dimension( 0 ), vs.getVoxelSize().dimension( 1 ) ), vs.getVoxelSize().dimension( 2 ) );
		
		final double calX = vs.getVoxelSize().dimension( 0 ) / minResolution;
		final double calY = vs.getVoxelSize().dimension( 1 ) / minResolution;
		final double calZ = vs.getVoxelSize().dimension( 2 ) / minResolution;

		final AffineTransform3D m = new AffineTransform3D();
		m.set( calX, 0.0f, 0.0f, 0.0f, 
			   0.0f, calY, 0.0f, 0.0f,
			   0.0f, 0.0f, calZ, 0.0f );
		final ViewTransform vt = new ViewTransformAffine( "Calibration", m );
		vr.preconcatenateTransform( vt );

		final AffineTransform3D translation = new AffineTransform3D();

		if ( tile.hasLocation() )
		{
			translation.set( tile.getLocation()[ 0 ] / calX, 0, 3 );
			translation.set( tile.getLocation()[ 1 ] / calY, 1, 3 );
			translation.set( tile.getLocation()[ 2 ] / calZ, 2, 3 );
		}

		vr.preconcatenateTransform( new ViewTransformAffine( "Translation", translation ) );

		vr.updateModel();		
		
		registrations.add( vr );
	}

	final SequenceDescription sd = new SequenceDescription( timepoints, setups, imgLoader, missingViews );
	final SpimData data = new SpimData( new File( "" ), sd, new ViewRegistrations( registrations ) );

	return data;
}
 

/**
 * Creates an {@link Img} of type {@link DoubleType} with the given
 * dimensions.
 */
public Img<DoubleType> img(final long[] dims) {
	return img(new FinalDimensions(dims), new DoubleType());
}