net.imglib2.interpolation.randomaccess.NLinearInterpolatorFactory Java Examples

@Test
public void testScaling() {
	Img<ByteType> in = generateByteArrayTestImg(true, new long[] { 10, 10 });
	double[] scaleFactors = new double[] { 2, 2 };
	@SuppressWarnings("unchecked")
	RandomAccessibleInterval<ByteType> out = (RandomAccessibleInterval<ByteType>) ops.run(DefaultScaleView.class, in,
		scaleFactors, new NLinearInterpolatorFactory<ByteType>());

	assertEquals(out.dimension(0), 20);
	assertEquals(out.dimension(1), 20);

	RandomAccess<ByteType> inRA = in.randomAccess();
	RandomAccess<ByteType> outRA = out.randomAccess();
	inRA.setPosition(new long[] { 5, 5 });
	outRA.setPosition(new long[] { 10, 10 });
	assertEquals(inRA.get().get(), outRA.get().get());

}
 

/**
 *
 * @param reader container
 * @param dataset dataset
 * @param transform transforms voxel data into real world coordinates
 * @param priority in fetching queue
 * @param name initialize with this name
 * @param <T> data type
 * @param <V> viewer type
 * @return {@link DataSource}
 * @throws IOException if any N5 operation throws {@link IOException}
 */
public static <T extends NativeType<T> & RealType<T>, V extends Volatile<T> & NativeType<V> & RealType<V>>
DataSource<T, V> openRawAsSource(
		final N5Reader reader,
		final String dataset,
		final AffineTransform3D transform,
		final SharedQueue queue,
		final int priority,
		final String name) throws IOException, ReflectionException {
	return openScalarAsSource(
			reader,
			dataset,
			transform,
			queue,
			priority,
			i -> i == Interpolation.NLINEAR
			     ? new NLinearInterpolatorFactory<>()
			     : new NearestNeighborInterpolatorFactory<>(),
			i -> i == Interpolation.NLINEAR
			     ? new NLinearInterpolatorFactory<>()
			     : new NearestNeighborInterpolatorFactory<>(),
			name
	                         );
}
 

@Override
public String call() throws Exception 
{
	final NLinearInterpolatorFactory< FloatType > f = new NLinearInterpolatorFactory< FloatType >();
	
	// make the interpolators and get the transformations
	final RealRandomAccess< FloatType > ir = Views.interpolate( Views.extendMirrorSingle( img ), f ).realRandomAccess();
	final RealRandomAccess< FloatType > wr = blending.realRandomAccess();

	final Cursor< FloatType > cursor = Views.iterable( transformedImg ).localizingCursor();
	final Cursor< FloatType > cursorW = Views.iterable( weightImg ).cursor();

	final float[] s = new float[ 3 ];
	final float[] t = new float[ 3 ];

	cursor.jumpFwd( portion.getStartPosition() );
	cursorW.jumpFwd( portion.getStartPosition() );

	for ( int j = 0; j < portion.getLoopSize(); ++j )
		loop( cursor, cursorW, ir, wr, transform, s, t, offsetX, offsetY, offsetZ, imgSizeX, imgSizeY, imgSizeZ );

	return portion + " finished successfully (transform input & precompute weights).";
}
 

@Override
public String call() throws Exception 
{
	final NLinearInterpolatorFactory< FloatType > f = new NLinearInterpolatorFactory< FloatType >();
	
	// make the interpolators and get the transformations
	final RealRandomAccess< FloatType > ir = Views.interpolate( Views.extendMirrorSingle( img ), f ).realRandomAccess();
	final Cursor< FloatType > cursor = Views.iterable( transformedImg ).localizingCursor();

	final float[] s = new float[ 3 ];
	final float[] t = new float[ 3 ];
	
	cursor.jumpFwd( portion.getStartPosition() );
	
	for ( int j = 0; j < portion.getLoopSize(); ++j )
		loop( cursor, ir, transform, s, t, offsetX, offsetY, offsetZ, imgSizeX, imgSizeY, imgSizeZ );
	
	return portion + " finished successfully (transform input & no weights).";
}
 

@Override
public RealRandomAccess<FloatType> realRandomAccess( final RealInterval interval )
{
	return Views.interpolate(
			Views.extendZero( this.contentBasedImg ),
			new NLinearInterpolatorFactory< FloatType >()
			).realRandomAccess( interval );
}
 

@SuppressWarnings({ "unused", "unchecked" })
@Test
public void testOutOfBoundsFactoryIsNull() {
	Img<ByteType> in = generateByteArrayTestImg(true, new long[] { 10, 10 });
	double[] scaleFactors = new double[] { 2, 2 };
	NLinearInterpolatorFactory<ByteType> nLinearInterpolatorFactory =
		new NLinearInterpolatorFactory<ByteType>();
	RandomAccessibleInterval<ByteType> out = (RandomAccessibleInterval<ByteType>) ops.run(DefaultScaleView.class, in,
		scaleFactors, nLinearInterpolatorFactory, null);
}
 

@Test(expected = IllegalArgumentException.class)
public void testContingency() {
	Img<ByteType> in = generateByteArrayTestImg(true, new long[] { 10, 10 });
	double[] scaleFactors = new double[] { 2, 2, 2 };
	ops.run(DefaultScaleView.class, in, scaleFactors,
		new NLinearInterpolatorFactory<ByteType>());
}
 

protected static < T extends RealType< T > & NativeType< T > > CompositeImage createOverlay( final T targetType, final ImagePlus imp1, final ImagePlus imp2, final InvertibleBoundable finalModel1, final InvertibleBoundable finalModel2, final int dimensionality ) 
{
	final ArrayList< ImagePlus > images = new ArrayList<ImagePlus>();
	images.add( imp1 );
	images.add( imp2 );
	
	final ArrayList< InvertibleBoundable > models = new ArrayList<InvertibleBoundable>();
	models.add( finalModel1 );
	models.add( finalModel2 );
	
	return createOverlay( targetType, images, models, dimensionality, 1, new NLinearInterpolatorFactory<FloatType>() );
}
 

@Override
public RealRandomAccess<FloatType> realRandomAccess()
{ 
	return Views.interpolate(
		Views.extendZero( this.contentBasedImg ),
		new NLinearInterpolatorFactory< FloatType >()
		).realRandomAccess();
}
 

public < T extends RealType< T > > InterpolatorFactory< T, RandomAccessible< T > > getInterpolatorFactory( final T type )
{
	if ( getInterpolation() == 0 )
		return new NearestNeighborInterpolatorFactory<T>();
	else
		return new NLinearInterpolatorFactory< T >();
}
 

final static private < T extends RealType< T > >InterpolatorFactory< RealComposite< T >, RandomAccessible< RealComposite< T > > > interpolatorFactory( final Interpolation interpolation )
{
	switch ( interpolation )
	{
	case NN:
		return new NearestNeighborInterpolatorFactory< RealComposite< T > >();
	default:
		return new NLinearInterpolatorFactory< RealComposite< T > >();
	}
}
 

private static <R extends RealType<R>, T extends NativeType<T> & RealType<T>> RealRandomAccessible<T> getInterpolatedDistanceTransformMask(
		final RandomAccessibleInterval<R> dt1,
		final RandomAccessibleInterval<R> dt2,
		final double distance,
		final T targetValue,
		final AffineTransform3D transformToSource)
{
	final RandomAccessibleInterval<R> distanceTransformStack = Views.stack(dt1, dt2);

	final R extendValue = Util.getTypeFromInterval(distanceTransformStack).createVariable();
	extendValue.setReal(extendValue.getMaxValue());
	final RealRandomAccessible<R> interpolatedDistanceTransform = Views.interpolate(
			Views.extendValue(distanceTransformStack, extendValue),
			new NLinearInterpolatorFactory<>()
		);

	final RealRandomAccessible<R> scaledInterpolatedDistanceTransform = RealViews.affineReal(
			interpolatedDistanceTransform,
			new Scale3D(1, 1, -distance)
		);

	final T emptyValue = targetValue.createVariable();
	final RealRandomAccessible<T> interpolatedShape = Converters.convert(
			scaledInterpolatedDistanceTransform,
			(in, out) -> out.set(in.getRealDouble() <= 0 ? targetValue : emptyValue),
			emptyValue.createVariable()
		);

	return RealViews.affineReal(interpolatedShape, transformToSource);
}
 

/**
 *
 * @param meta
 * @param transform
 * @param dataExtension
 * @param extension
 * @param name
 * @param priority
 * @param channelDimension
 * @param channels
 * @throws IOException
 * @throws DataTypeNotSupported
 */
private N5ChannelDataSource(
		final N5Meta meta,
		final AffineTransform3D transform,
		final D dataExtension,
		final T extension,
		final String name,
		final SharedQueue queue,
		final int priority,
		final int channelDimension,
		final long[] channels) throws
		IOException, DataTypeNotSupported {

	final ImagesWithTransform<D, T>[] data = getData(
			meta.reader(),
			meta.dataset(),
			transform,
			queue,
			priority);
	final RandomAccessibleIntervalDataSource.DataWithInvalidate<D, T> dataWithInvalidate = RandomAccessibleIntervalDataSource.asDataWithInvalidate(data);
	this.meta = meta;
	this.channelDimension = channelDimension;
	this.name = name;
	this.transforms = dataWithInvalidate.transforms;
	this.invalidate = dataWithInvalidate.invalidate;

	this.channels = channels == null ? range((int) dataWithInvalidate.data[0].dimension(channelDimension)) : channels;
	this.numChannels = this.channels.length;

	this.intervals = dataWithInvalidate.data;
	extension.setValid(true);
	this.data = collapseDimension(dataWithInvalidate.data, this.channelDimension, this.channels, dataExtension);
	this.viewerData = collapseDimension(dataWithInvalidate.viewData, this.channelDimension, this.channels, extension);

	this.interpolation = ipol -> new NearestNeighborInterpolatorFactory<>();
	this.viewerInterpolation = ipol -> Interpolation.NLINEAR.equals(ipol) ? new NLinearInterpolatorFactory<>() : new NearestNeighborInterpolatorFactory<>();

	LOG.debug("Channel dimension {} has {} channels", channelDimension, numChannels);
}
 

private static <T extends RealType<T>> Function<Interpolation, InterpolatorFactory<T, RandomAccessible<T>>>
realTypeInterpolation()
{
	return i -> i.equals(Interpolation.NLINEAR)
	            ? new NLinearInterpolatorFactory<>()
	            : new NearestNeighborInterpolatorFactory<>();
}
 

@Override
public InterpolatorFactory<T, RandomAccessible<T>> apply(final Interpolation t)
{
	return t.equals(Interpolation.NLINEAR)
	       ? new NLinearInterpolatorFactory<>()
	       : new NearestNeighborInterpolatorFactory<>();
}
 

public LinearIntensityMap( final RandomAccessibleInterval< T > source )
{
	this( source, new NLinearInterpolatorFactory< RealComposite< T > >() );
}
 

/**
 * Extracts the PSF by averaging the local neighborhood RANSAC correspondences
 * @param size - the size in which the psf is extracted (in pixel units, z-scaling is ignored)
 * @return - the psf, NOT z-scaling corrected
 */
protected static < T extends RealType< T > & NativeType< T > > ArrayImg< T, ? > extractPSFLocal(
		final RandomAccessibleInterval< T > img,
		final ArrayList< double[] > locations,
		final long[] size )
{
	final int numDimensions = size.length;
	
	final ArrayImg< T, ? > psf = new ArrayImgFactory< T >().create( size, Views.iterable( img ).firstElement() );
	
	// Mirror produces some artifacts ... so we use periodic
	final RealRandomAccess< T > interpolator =
			Views.interpolate( Views.extendPeriodic( img ), new NLinearInterpolatorFactory< T >() ).realRandomAccess();
	
	final ArrayLocalizingCursor< T > psfCursor = psf.localizingCursor();
	
	final long[] sizeHalf = size.clone();
	for ( int d = 0; d < numDimensions; ++d )
		sizeHalf[ d ] /= 2;
	
	final int[] tmpI = new int[ size.length ];
	final double[] tmpD = new double[ size.length ];

	for ( final double[] position : locations )
	{
		psfCursor.reset();
		
		while ( psfCursor.hasNext() )
		{
			psfCursor.fwd();
			psfCursor.localize( tmpI );

			for ( int d = 0; d < numDimensions; ++d )
				tmpD[ d ] = tmpI[ d ] - sizeHalf[ d ] + position[ d ];
			
			interpolator.setPosition( tmpD );
			
			psfCursor.get().add( interpolator.get() );
		}
	}

	return psf;
}
 

/**
 * @return the default interpolator factory for warp field instances.
 */
public static InterpolatorFactory<RealComposite<DoubleType>, RandomAccessible<RealComposite<DoubleType>>> getDefaultInterpolatorFactory() {
    return new NLinearInterpolatorFactory<>();
}
 

public static < T extends RealType< T > & NativeType< T > > ImagePlus createReRegisteredSeries( final T targetType, final ImagePlus imp, final ArrayList<InvertibleBoundable> models, final int dimensionality )
{
	final int numImages = imp.getNFrames();

	// the size of the new image
	final int[] size = new int[ dimensionality ];
	// the offset relative to the output image which starts with its local coordinates (0,0,0)
	final double[] offset = new double[ dimensionality ];

	final int[][] imgSizes = new int[ numImages ][ dimensionality ];
	
	for ( int i = 0; i < numImages; ++i )
	{
		imgSizes[ i ][ 0 ] = imp.getWidth();
		imgSizes[ i ][ 1 ] = imp.getHeight();
		if ( dimensionality == 3 )
			imgSizes[ i ][ 2 ] = imp.getNSlices();
	}
	
	// estimate the boundaries of the output image and the offset for fusion (negative coordinates after transform have to be shifted to 0,0,0)
	Fusion.estimateBounds( offset, size, imgSizes, models, dimensionality );
			
	// for output
	final ImgFactory< T > f = new ImagePlusImgFactory< T >();
	// the composite
	final ImageStack stack = new ImageStack( size[ 0 ], size[ 1 ] );

	for ( int t = 1; t <= numImages; ++t )
	{
		for ( int c = 1; c <= imp.getNChannels(); ++c )
		{
			final Img<T> out = f.create( size, targetType );
			final Img< FloatType > in = ImageJFunctions.convertFloat( Hyperstack_rearranger.getImageChunk( imp, c, t ) );

			fuseChannel( out, Views.interpolate( Views.extendZero( in ), new NLinearInterpolatorFactory< FloatType >() ), offset, models.get( t - 1 ) );

			try
			{
				final ImagePlus outImp = ((ImagePlusImg<?,?>)out).getImagePlus();
				for ( int z = 1; z <= out.dimension( 2 ); ++z )
					stack.addSlice( imp.getTitle(), outImp.getStack().getProcessor( z ) );
			} 
			catch (ImgLibException e) 
			{
				Log.error( "Output image has no ImageJ type: " + e );
			}				
		}
	}
	
	//convertXYZCT ...
	ImagePlus result = new ImagePlus( "registered " + imp.getTitle(), stack );
	
	// numchannels, z-slices, timepoints (but right now the order is still XYZCT)
	if ( dimensionality == 3 )
	{
		result.setDimensions( size[ 2 ], imp.getNChannels(), imp.getNFrames() );
		result = OverlayFusion.switchZCinXYCZT( result );
		return CompositeImageFixer.makeComposite( result, CompositeImage.COMPOSITE );
	}
	//Log.info( "ch: " + imp.getNChannels() );
	//Log.info( "slices: " + imp.getNSlices() );
	//Log.info( "frames: " + imp.getNFrames() );
	result.setDimensions( imp.getNChannels(), 1, imp.getNFrames() );
	
	if ( imp.getNChannels() > 1 )
		return CompositeImageFixer.makeComposite( result, CompositeImage.COMPOSITE );
	return result;
}
 

protected static < T extends RealType< T > & NativeType< T > > ImagePlus fuse( final T targetType, final ImagePlus imp1, final ImagePlus imp2, final ArrayList<InvertibleBoundable> models, final StitchingParameters params )
{
	final ArrayList<ImagePlus> images = new ArrayList< ImagePlus >();
	images.add( imp1 );
	images.add( imp2 );
	
	if ( params.fusionMethod < 6 )
	{
		ImagePlus imp = Fusion.fuse( targetType, images, models, params.dimensionality, params.subpixelAccuracy, params.fusionMethod, null, false, params.ignoreZeroValuesFusion, params.displayFusion );
		return imp;
	}
	else if ( params.fusionMethod == 6 ) // overlay
	{
		// images are always the same, we just trigger different timepoints
		final InterpolatorFactory< FloatType, RandomAccessible< FloatType > > factory;
		
		if ( params.subpixelAccuracy )
			factory  = new NLinearInterpolatorFactory<FloatType>();
		else
			factory  = new NearestNeighborInterpolatorFactory< FloatType >();
	
		// fuses the first timepoint but estimates the boundaries for all timepoints as it gets all models
		final CompositeImage timepoint0 = OverlayFusion.createOverlay( targetType, images, models, params.dimensionality, 1, factory );
		
		if ( imp1.getNFrames() > 1 )
		{
			final ImageStack stack = new ImageStack( timepoint0.getWidth(), timepoint0.getHeight() );
			
			// add all slices of the first timepoint
			for ( int c = 1; c <= timepoint0.getStackSize(); ++c )
				stack.addSlice( "", timepoint0.getStack().getProcessor( c ) );
			
			//"Overlay into composite image"
			for ( int f = 2; f <= imp1.getNFrames(); ++f )
			{
				final CompositeImage tmp = OverlayFusion.createOverlay( targetType, images, models, params.dimensionality, f, factory );
				
				// add all slices of the first timepoint
				for ( int c = 1; c <= tmp.getStackSize(); ++c )
					stack.addSlice( "", tmp.getStack().getProcessor( c ) );					
			}
			
			//convertXYZCT ...
			ImagePlus result = new ImagePlus( params.fusedName, stack );
			
			// numchannels, z-slices, timepoints (but right now the order is still XYZCT)
			result.setDimensions( timepoint0.getNChannels(), timepoint0.getNSlices(), imp1.getNFrames() );
			return CompositeImageFixer.makeComposite( result, CompositeImage.COMPOSITE );
		}
		else
		{
			timepoint0.setTitle( params.fusedName );
			return timepoint0;
		}
	}
	else
	{
		//"Do not fuse images"
		return null;
	}
}
 

public <T extends RealType<T>, S extends RealType<S>> void calculateCrossCorr(RandomAccessibleInterval<T> img1, RandomAccessibleInterval<S> img2, 
		long minOverlapPx, boolean interpolateSubpixel)
{
	Pair<Interval, Interval> intervals = PhaseCorrelation2Util.getOverlapIntervals(img1, img2, shift);
	
	// no overlap found
	if (intervals == null) {
		crossCorr = Double.NEGATIVE_INFINITY;
		nPixel = 0;
		return;
	}
	
	nPixel = 1;
	for (int i = 0; i< intervals.getA().numDimensions(); i++){
		nPixel *= intervals.getA().dimension(i);
	}
	
	if (nPixel < minOverlapPx){
		crossCorr = Double.NEGATIVE_INFINITY;
		nPixel = 0;
		return;
	}

	// for subpixel move the underlying Img2 by the subpixel offset
	if ( subpixelShift != null && interpolateSubpixel )
	{
		RealRandomAccessible< S > rra = Views.interpolate( Views.extendMirrorSingle( img2 ), new NLinearInterpolatorFactory< S >() );

		InvertibleRealTransform transform = null;

		// e.g. subpixel = (-0.4, 0.1, -0.145)
		final double tx = subpixelShift.getDoublePosition( 0 ) - shift.getDoublePosition( 0 );
		final double ty = subpixelShift.getDoublePosition( 1 ) - shift.getDoublePosition( 1 );

		if ( rra.numDimensions() == 2 )
			transform = new Translation2D( -tx, -ty ); // -relative subpixel shift only
		else if ( rra.numDimensions() == 3 )
			transform = new Translation3D( -tx, -ty, shift.getDoublePosition( 2 ) - subpixelShift.getDoublePosition( 2 ) ); // -relative subpixel shift only

		img2 = Views.interval( Views.raster( RealViews.transform( rra, transform ) ), img2 );
	}

	// calculate cross correlation.
	// note that the overlap we calculate assumes zero-min input
	crossCorr = PhaseCorrelation2Util.getCorrelation(
			Views.zeroMin( Views.interval(Views.zeroMin(img1), intervals.getA())),
			Views.zeroMin( Views.interval(Views.zeroMin(img2), intervals.getB()))
		);
	
}
 

public static void main(String[] args)
{
	Img< FloatType > a = ImgLib2Util.openAs32Bit( new File( "73.tif.zip" ) );
	Img< FloatType > b = ImgLib2Util.openAs32Bit( new File( "74.tif.zip" ) );

	TranslationGet t1 = new Translation3D();
	TranslationGet t2 = new Translation3D(460, 0, 0);
	ArrayList< Pair< RealInterval, AffineGet > > views = new ArrayList<Pair<RealInterval, AffineGet>>();
	views.add( new ValuePair< RealInterval, AffineGet >( a, t1 ) );
	views.add( new ValuePair< RealInterval, AffineGet >( b, t2 ) );

	RealInterval overlap = BoundingBoxMaximalGroupOverlap.getMinBoundingIntervalSingle( views );

	final RealInterval transformed1 = TransformTools.applyTranslation( a, t1, new boolean[] {false, false, false} );
	final RealInterval transformed2 = TransformTools.applyTranslation( b, t2, new boolean[] {false, false, false} );

	// get overlap in images' coordinates
	final RealInterval localOverlap1 = TransformTools.getLocalOverlap( transformed1, overlap );
	final RealInterval localOverlap2 = TransformTools.getLocalOverlap( transformed2, overlap );

	// round to integer interval
	final Interval interval1 = TransformTools.getLocalRasterOverlap( localOverlap1 );
	final Interval interval2 = TransformTools.getLocalRasterOverlap( localOverlap2 );

	//final WarpFunction warp = new TranslationWarp(3);
	final WarpFunction warp = new RigidWarp(3);
	//final WarpFunction warp = new AffineWarp( 3 );

	// rotate second image
	AffineTransform3D rot = new AffineTransform3D();
	rot.rotate( 2, 2 * Math.PI / 180 );
	RandomAccessibleInterval< FloatType > rotated = Views.interval(
			RealViews.affine( 
					Views.interpolate( Views.extendBorder( Views.zeroMin( Views.interval( b, interval2 ) ) ), new NLinearInterpolatorFactory<>() ),
					rot.copy() ),
			interval2);

	// show input
	new ImageJ();
	ImageJFunctions.show( Views.interval( a,  interval1 ), "target" );
	ImageJFunctions.show( rotated, "in");

	// downsample input
	RandomAccessibleInterval< FloatType > simple2x1 = Downsample.simple2x( Views.zeroMin( Views.interval( a, interval1 ) ), new ArrayImgFactory<>(), new boolean[] {true, true, false} );
	RandomAccessibleInterval< FloatType > simple2x2 = Downsample.simple2x( Views.zeroMin( Views.interval( rotated, interval2 ) ), new ArrayImgFactory<>(), new boolean[] {true, true, false} );

	// align

	//Align< FloatType > lk = new Align<>( Views.zeroMin( Views.interval( a, interval1 ) ), new ArrayImgFactory<>(), warp );
	Align< FloatType > lk = new Align<>( simple2x1, new ArrayImgFactory<>(), warp );
	//System.out.println( Util.printCoordinates( lk.align( Views.zeroMin( Views.interval( b, interval2 ) ), 100, 0.01 ).getRowPackedCopy() ) );
	//final AffineTransform transform = lk.align( Views.zeroMin( rotated ), 100, 0.01 );
	final AffineTransform transform = lk.align( simple2x2, 100, 0.01 );

	final AffineTransform scale = new AffineTransform( 3 );
	scale.set( 2, 0, 0 );
	scale.set( 1, 1, 1 );

	transform.preConcatenate( scale );

	// transformation matrix
	System.out.println( Util.printCoordinates( transform.getRowPackedCopy() ) );

	// correct input and show
	RandomAccessibleInterval< FloatType > backRotated = Views.interval(
			RealViews.affine( 
					Views.interpolate( Views.extendBorder( Views.zeroMin( Views.interval( b, interval2 ) ) ), new NLinearInterpolatorFactory<>() ),
					rot.copy().preConcatenate( transform ).copy() ),
			interval2);

	ImageJFunctions.show( backRotated, "out" );

	// constructor needs column packed matrix, therefore the transpose
	Matrix mt = new Matrix( transform.getRowPackedCopy(), 4).transpose();
	Matrix rigid = mt.getMatrix( 0, 2, 0, 2 );

	// check whether result is rotation matrix (det == +-1, orthogonal)
	System.out.println( rigid.det() );
	System.out.println( Util.printCoordinates( rigid.times( rigid.transpose() ).getRowPackedCopy() ) );
}
 

/**
 * Compute the pixel-wise difference between an affine-transformed source
 * image and a target image.
 *
 * @param source
 *            The source image.
 * @param transform
 *            A coordinate transformation to apply to the source image.
 * @param target
 *            The target image.
 * @param difference
 *            Output image. The pixel-wise difference between the
 *            transformed source image and the target image is stored here.
 * @param service
 *            thread pool for difference calculation
 * @param nTasks
 *            number of image parts that are processed in parallel
 * @param <T> pixel type source
 * @param <S> pixel type target
 */
public static < T extends RealType< T >,  S extends RealType< S > > void computeDifference(
		final RandomAccessible< T > source,
		final AffineTransform transform,
		final RandomAccessible< T > target,
		final RandomAccessibleInterval< S > difference,
		final ExecutorService service,
		final int nTasks)
{
	final RealRandomAccessible< T > interpolated = Views.interpolate( source, new NLinearInterpolatorFactory< T >() );
	final RandomAccessible< T > warped = RealViews.affine( interpolated, transform );

	final long stepSize = Views.iterable( difference ).size() / nTasks;

	final List<Callable< Void >> tasks = new ArrayList<>();
	final AtomicInteger ai = new AtomicInteger( 0 );
	for (int iO = 0; iO<nTasks; iO++)
	{
		tasks.add( new Callable< Void >()
		{
			@Override
			public Void call() throws Exception
			{
				final int i = ai.getAndIncrement();
				final Cursor< T > cw = Views.flatIterable( Views.interval( warped, difference ) ).cursor();
				final Cursor< T > ct = Views.flatIterable( Views.interval( target, difference ) ).cursor();
				final Cursor< S > cd = Views.flatIterable( difference ).cursor();

				cw.jumpFwd( stepSize * i );
				ct.jumpFwd( stepSize * i );
				cd.jumpFwd( stepSize * i );

				final long end = i == nTasks - 1 ? Views.iterable( difference ).size() - stepSize * i : stepSize;
				int count = 0;
				while (count++ < end)
				{
					cd.next().setReal( ( cw.next().getRealDouble() - ct.next().getRealDouble() ));
				}
				return null;
			}
		} );
	}

	try
	{
		List< Future< Void > > futures = service.invokeAll( tasks );
		for (Future< Void > f: futures)
			f.get();
	}
	catch ( InterruptedException | ExecutionException e )
	{
		e.printStackTrace();
	}
}
 

public double getCurrentCorrelation(final RandomAccessibleInterval< T > image)
{
	final RealRandomAccessible< T > interpolated = Views.interpolate( Views.extendBorder( image ), new NLinearInterpolatorFactory< T >() );
	final RandomAccessible< T > warped = RealViews.affine( interpolated, currentTransform );
	return PhaseCorrelation2Util.getCorrelation( Views.interval( warped, template ), template );
}
 

public static void main(String[] args)
{
	RandomAccessibleInterval< FloatType > a = ImgLib2Util.openAs32Bit( new File( "73.tif.zip" ) );
	RandomAccessibleInterval< FloatType > b = ImgLib2Util.openAs32Bit( new File( "74.tif.zip" ) );
	
	long slice = 40;
	ImageJFunctions.show( a );
	
	a = Views.zeroMin( Views.hyperSlice( a, 2, slice ));
	b = Views.zeroMin( Views.hyperSlice( b, 2, slice ));

	TranslationGet t1 = new Translation2D();
	TranslationGet t2 = new Translation2D(460, 0);
	ArrayList< Pair< RealInterval, AffineGet > > views = new ArrayList<Pair<RealInterval, AffineGet>>();
	views.add( new ValuePair< RealInterval, AffineGet >( a, t1 ) );
	views.add( new ValuePair< RealInterval, AffineGet >( b, t2 ) );

	RealInterval overlap = BoundingBoxMaximalGroupOverlap.getMinBoundingIntervalSingle( views );

	final RealInterval transformed1 = TransformTools.applyTranslation( a, t1, new boolean[] {false, false} );
	final RealInterval transformed2 = TransformTools.applyTranslation( b, t2, new boolean[] {false, false} );

	// get overlap in images' coordinates
	final RealInterval localOverlap1 = TransformTools.getLocalOverlap( transformed1, overlap );
	final RealInterval localOverlap2 = TransformTools.getLocalOverlap( transformed2, overlap );

	// round to integer interval
	final Interval interval1 = TransformTools.getLocalRasterOverlap( localOverlap1 );
	final Interval interval2 = TransformTools.getLocalRasterOverlap( localOverlap2 );

	//final WarpFunction warp = new TranslationWarp(3);
	final WarpFunction warp = new RigidWarp(2);
	//final WarpFunction warp = new AffineWarp( 3 );

	// rotate second image
	AffineTransform2D rot = new AffineTransform2D();
	rot.rotate( 1.4 * Math.PI / 180 );
	RandomAccessibleInterval< FloatType > rotated = Views.interval(
			RealViews.affine( 
					Views.interpolate( Views.extendMirrorSingle( Views.zeroMin( Views.interval( b, interval2 ) ) ), new NLinearInterpolatorFactory<>() ),
					rot.copy() ),
			interval2);

	// show input
	new ImageJ();
	ImageJFunctions.show( Views.interval( a,  interval1 ) );
	ImageJFunctions.show( rotated );

	// downsample input
	RandomAccessibleInterval< FloatType > simple2x1 = Downsample.simple2x( Views.zeroMin( Views.interval( a, interval1 ) ), new ArrayImgFactory<>(), new boolean[] {false, false} );
	RandomAccessibleInterval< FloatType > simple2x2 = Downsample.simple2x( Views.zeroMin( Views.interval( rotated, interval2 ) ), new ArrayImgFactory<>(), new boolean[] {false, false} );

	// align

	//Align< FloatType > lk = new Align<>( Views.zeroMin( Views.interval( a, interval1 ) ), new ArrayImgFactory<>(), warp );
	Align< FloatType > lk = new Align<>( simple2x1, new ArrayImgFactory<>(), warp );
	//System.out.println( Util.printCoordinates( lk.align( Views.zeroMin( Views.interval( b, interval2 ) ), 100, 0.01 ).getRowPackedCopy() ) );
	//final AffineTransform transform = lk.align( Views.zeroMin( rotated ), 100, 0.01 );
	final AffineTransform transform = lk.align( simple2x2, 100, 0.1 );

	// transformation matrix
	System.out.println( Util.printCoordinates( transform.getRowPackedCopy() ) );

	// correct input and show
	RandomAccessibleInterval< FloatType > backRotated = Views.interval(
			RealViews.affine( 
					Views.interpolate( Views.extendMirrorSingle( Views.zeroMin( Views.interval( b, interval2 ) ) ), new NLinearInterpolatorFactory<>() ),
					rot.copy().preConcatenate( transform ).copy() ),
			interval2);

	ImageJFunctions.show( backRotated );

	// constructor needs column packed matrix, therefore the transpose
	Matrix mt = new Matrix( transform.getRowPackedCopy(), 3).transpose();
	Matrix rigid = mt.getMatrix( 0, 1, 0, 1 );

	// check whether result is rotation matrix (det == +-1, orthogonal)
	System.out.println( rigid.det() );
	System.out.println( Util.printCoordinates( rigid.times( rigid.transpose() ).getRowPackedCopy() ) );
}
 

@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 );
	
}