public Align(final RandomAccessibleInterval< T > template, final ImgFactory< FloatType > factory, WarpFunction model)
{
	this.template = template;

	n = template.numDimensions();
	warpFunction = model;
	numParameters = warpFunction.numParameters();
	
	currentTransform = new AffineTransform( n );
	
	final long[] dim = new long[n + 1];
	for ( int d = 0; d < n; ++d )
		dim[d] = template.dimension( d );
	dim[n] = n;
	final Img< FloatType > gradients = factory.create( dim, new FloatType() );
	gradients( Views.extendBorder( template ), gradients );

	dim[n] = numParameters;
	descent = factory.create( dim, new FloatType() );
	computeSteepestDescents( gradients, warpFunction, descent );

	Hinv = computeInverseHessian( descent );

	error = factory.create( template, new FloatType() );
}
 

public static <T extends RealType<T>, S extends RealType<S>, R extends RealType<R>, C extends ComplexType<C>> RandomAccessibleInterval<R> calculatePCM(
		RandomAccessibleInterval<T> img1, RandomAccessibleInterval<S> img2, int[] extension,
		ImgFactory<R> factory, R type, ImgFactory<C> fftFactory, C fftType, ExecutorService service){

	
	// TODO: Extension absolute per dimension in pixels, i.e. int[] extension
	// TODO: not bigger than the image dimension because the second mirroring is identical to the image
	
	Dimensions extSize = PhaseCorrelation2Util.getExtendedSize(img1, img2, extension);
	long[] paddedDimensions = new long[extSize.numDimensions()];
	long[] fftSize = new long[extSize.numDimensions()];
	FFTMethods.dimensionsRealToComplexFast(extSize, paddedDimensions, fftSize);
	
	RandomAccessibleInterval<C> fft1 = fftFactory.create(fftSize, fftType);
	RandomAccessibleInterval<C> fft2 = fftFactory.create(fftSize, fftType);
	
	FFT.realToComplex(Views.interval(PhaseCorrelation2Util.extendImageByFactor(img1, extension), 
			FFTMethods.paddingIntervalCentered(img1, new FinalInterval(paddedDimensions))), fft1, service);
	FFT.realToComplex(Views.interval(PhaseCorrelation2Util.extendImageByFactor(img2, extension), 
			FFTMethods.paddingIntervalCentered(img2, new FinalInterval(paddedDimensions))), fft2, service);
	
	RandomAccessibleInterval<R> pcm = calculatePCMInPlace(fft1, fft2, factory, type, service);
	return pcm;
	
}
 

/**
 * @see ImgOpener#openImgs(Reader, ImgFactory, SCIFIOConfig)
 */
public static <T extends RealType<T>> List<SCIFIOImgPlus<T>> openAll(
		final Reader reader, final ImgFactory<T> imgFactory,
		final SCIFIOConfig config)
{
	final ImgOpener opener = opener();
	List<SCIFIOImgPlus<T>> imgPlus = null;
	try {
		imgPlus = opener.openImgs(reader, imgFactory, config);
		register(imgPlus, opener);
	}
	catch (final ImgIOException e) {
		openError(reader.getMetadata().getSourceLocation(), e);
	}
	return imgPlus;
}
 

/**
 * This method creates a noise image that has a specified mean. Every pixel
 * has a value uniformly distributed around mean with the maximum spread
 * specified.
 *
 * @return IllegalArgumentException if specified means and spreads are not
 *         valid
 */
public static <T extends RealType<T> & NativeType<T>>
	Img<T> produceMeanBasedNoiseImage(T type, int width,
		int height, double mean, double spread, double[] smoothingSigma,
		long seed) throws IllegalArgumentException
{
	if (mean < spread || (mean + spread) > type.getMaxValue()) {
		throw new IllegalArgumentException(
			"Mean must be larger than spread, and mean plus spread must be smaller than max of the type");
	}
	// create the new image
	ImgFactory<T> imgFactory = new ArrayImgFactory<>(type);
	Img<T> noiseImage = imgFactory.create(width, height);

	Random r = new Random(seed);
	for (T value : Views.iterable(noiseImage)) {
		value.setReal(mean + ((r.nextDouble() - 0.5) * spread));
	}

	// TODO: call Ops filter.gauss instead
	return gaussianSmooth(noiseImage, smoothingSigma);
}
 

/**
 * Gaussian Smooth of the input image using intermediate float format.
 * 
 * @param <T>
 * @param img
 * @param sigma
 * @return
 */
public static <T extends RealType<T> & NativeType<T>>
	Img<T> gaussianSmooth(RandomAccessibleInterval<T> img,
		double[] sigma)
{
	Interval interval = Views.iterable(img);

	ImgFactory<T> outputFactory = new ArrayImgFactory<>(Util.getTypeFromInterval(img));
	final long[] dim = new long[img.numDimensions()];
	img.dimensions(dim);
	Img<T> output = outputFactory.create(dim);

	final long[] pos = new long[img.numDimensions()];
	Arrays.fill(pos, 0);
	Localizable origin = new Point(pos);

	ImgFactory<FloatType> tempFactory = new ArrayImgFactory<>(new FloatType());
	RandomAccessible<T> input = Views.extendMirrorSingle(img);
	Gauss.inFloat(sigma, input, interval, output, origin, tempFactory);

	return output;
}
 

/**
 * @see ImgOpener#openImgs(Reader, ImgFactory, SCIFIOConfig)
 */
public static <T extends RealType<T>> List<SCIFIOImgPlus<T>> openAll(
	final Reader reader, final T type, final ImgFactory<T> imgFactory,
	final SCIFIOConfig config)
{
	final ImgOpener opener = opener();
	List<SCIFIOImgPlus<T>> imgPlus = null;
	try {
		imgPlus = opener.openImgs(reader, imgFactory.imgFactory(type), config);
		register(imgPlus, opener);
	}
	catch (final ImgIOException e) {
		openError(reader.getMetadata().getSourceLocation(), e);
	}
	return imgPlus;
}
 

protected ImgFactory<? extends NativeType<?>> selectImgFactory(final SlideBook6MetaData meta) {
    int[] dims = meta.imageSize(0);
    long maxNumPixels = dims[0];
    maxNumPixels *= dims[1];
    maxNumPixels *= dims[2];

    String s = "Maximum number of pixels in any view: n=" + Long.toString(maxNumPixels) +
            " px ";

    if (maxNumPixels < Integer.MAX_VALUE) {
        IOFunctions.println(s + "< " + Integer.MAX_VALUE + ", using ArrayImg.");
        return new ArrayImgFactory<FloatType>();
    } else {
        IOFunctions.println(s + ">= " + Integer.MAX_VALUE + ", using CellImg.");
        return new CellImgFactory<FloatType>(256);
    }
}
 

/**
 * For a local initialization without the XML
 * 
 * @param path
 * @param fileNamePattern
 * @param imgFactory
 * @param layoutTP - 0 == one, 1 == one per file, 2 == all in one file
 * @param layoutChannels - 0 == one, 1 == one per file, 2 == all in one file
 * @param layoutIllum - 0 == one, 1 == one per file, 2 == all in one file
 * @param layoutAngles - 0 == one, 1 == one per file, 2 == all in one file
 */
public LegacyStackImgLoader(
		final File path, final String fileNamePattern, final ImgFactory< ? extends NativeType< ? > > imgFactory,
		final int layoutTP, final int layoutChannels, final int layoutIllum, final int layoutAngles,
		final AbstractSequenceDescription< ?, ?, ? > sequenceDescription )
{
	super();
	this.path = path;
	this.fileNamePattern = fileNamePattern;
	this.layoutTP = layoutTP;
	this.layoutChannels = layoutChannels;
	this.layoutIllum = layoutIllum;
	this.layoutAngles = layoutAngles;
	this.sequenceDescription = sequenceDescription;
	
	this.init( imgFactory );
}
 

final protected static Thread getCUDAThread2(
		final AtomicInteger ai, final ImgFactory< FloatType > blockFactory, final Block[] blocks, final int[] blockSize,
		final Img< FloatType > image, final Img< FloatType > result, final int deviceId, final Img< FloatType > kernel2 )
{
	final Thread cudaThread2 = new Thread( new Runnable()
	{
		public void run()
		{
			final Img< FloatType > block = blockFactory.create( Util.int2long( blockSize ), new FloatType() );

			int i;

			while ( ( i = ai.getAndIncrement() ) < blocks.length )
				convolve2BlockCUDA( blocks[ i ], deviceId, image, result, block, kernel2 );
		}
	});
	
	return cudaThread2;
}
 

public static < T extends RealType< T > > RandomAccessibleInterval< T > simple2x( final RandomAccessibleInterval<T> input, final ImgFactory< T > imgFactory, final boolean[] downsampleInDim )
{
	RandomAccessibleInterval< T > src = input;

	for ( int d = 0; d < input.numDimensions(); ++d )
		if ( downsampleInDim[ d ] )
		{
			final long dim[] = new long[ input.numDimensions() ];

			for ( int e = 0; e < input.numDimensions(); ++e )
			{
				if ( e == d )
					dim[ e ] = src.dimension( e ) / 2;
				else
					dim[ e ] = src.dimension( e );
			}

			final Img< T > img = imgFactory.create( dim, Views.iterable( input ).firstElement() );
			simple2x( src, img, d );
			src = img;
		}

	return src;
}
 

/**
 * This method creates a noise image that has a specified mean.
 * Every pixel has a value uniformly distributed around mean with
 * the maximum spread specified.
 *
 * @return a new noise image
 * @throws MissingPreconditionException if specified means and spreads are not valid
 */
public static <T extends RealType<T> & NativeType<T>> RandomAccessibleInterval<T> produceMeanBasedNoiseImage(T type, int width,
		int height, double mean, double spread, double[] smoothingSigma, long seed) throws MissingPreconditionException {
	if (mean < spread || (mean + spread) > type.getMaxValue()) {
		throw new MissingPreconditionException("Mean must be larger than spread, and mean plus spread must be smaller than max of the type");
	}
	// create the new image
	ImgFactory<T> imgFactory = new ArrayImgFactory<T>();
	RandomAccessibleInterval<T> noiseImage = imgFactory.create( new int[] {width, height}, type); // "Noise image");

	Random r = new Random(seed);
	for (T value : Views.iterable(noiseImage)) {
		value.setReal( mean + ( (r.nextDouble() - 0.5) * spread ) );
	}

	return gaussianSmooth(noiseImage, smoothingSigma);
}
 

/**
 * Gaussian Smooth of the input image using intermediate float format.
 * @param <T>
 * @param img
 * @param sigma
 * @return
 */
public static <T extends RealType<T> & NativeType<T>> RandomAccessibleInterval<T> gaussianSmooth(
		RandomAccessibleInterval<T> img, double[] sigma) {
	Interval interval = Views.iterable(img);

	ImgFactory<T> outputFactory = new ArrayImgFactory<T>(Util.getTypeFromInterval(img));
	final long[] dim = new long[ img.numDimensions() ];
	img.dimensions(dim);
	RandomAccessibleInterval<T> output = outputFactory.create( dim );

	final long[] pos = new long[ img.numDimensions() ];
	Arrays.fill(pos, 0);
	Localizable origin = new Point(pos);

	ImgFactory<FloatType> tempFactory = new ArrayImgFactory<FloatType>(new FloatType());
	RandomAccessible<T> input = Views.extendMirrorSingle(img);
	Gauss.inFloat(sigma, input, interval, output, origin, tempFactory);

	return output;
}
 

/**
 * @deprecated Use {@link #open(Location, ImgFactory)}.
 */
@Deprecated
public static <T extends RealType<T> & NativeType<T>> SCIFIOImgPlus<T> open(
	final Location source, final ImgFactory<T> imgFactory,
	@SuppressWarnings("unused") final T type)
{
	return open(opener(), source, imgFactory);
}
 

public static <T extends ComplexType<T> & NativeType<T>, S extends ComplexType<S> & NativeType <S>, R extends RealType<R>> RandomAccessibleInterval<R> calculatePCM(
		RandomAccessibleInterval<T> fft1, RandomAccessibleInterval<S> fft2, ImgFactory<R> factory, R type, ExecutorService service){
	
	long[] paddedDimensions = new long[fft1.numDimensions()];
	long[] realSize = new long[fft1.numDimensions()];
	
	FFTMethods.dimensionsComplexToRealFast(fft1, paddedDimensions, realSize);
	RandomAccessibleInterval<R> res = factory.create(realSize, type);
	
	final T typeT = Views.iterable(fft1).firstElement().createVariable();
	final S typeS = Views.iterable(fft2).firstElement().createVariable();
	RandomAccessibleInterval< T > fft1Copy;
	RandomAccessibleInterval< S > fft2Copy;

	try
	{
		fft1Copy = factory.imgFactory( typeT ).create(fft1, typeT );
		fft2Copy = factory.imgFactory( typeS ).create(fft2, typeS );
	}
	catch ( IncompatibleTypeException e )
	{
		throw new RuntimeException( "Cannot instantiate Img for type " + typeS.getClass().getSimpleName() + " or " + typeT.getClass().getSimpleName() );
	}
	
	
	calculatePCM(fft1, fft1Copy, fft2, fft2Copy, res, service);
	
	return res;
}
 

private static <T extends RealType<T> & NativeType<T>> List<SCIFIOImgPlus<T>>
openAll(final ImgOpener opener, final Location source, final ImgFactory<T> imgFactory,
        final SCIFIOConfig config)
{
	List<SCIFIOImgPlus<T>> imgPlus = null;
	try {
		imgPlus = opener.openImgs(source, imgFactory, config);
		register(imgPlus, opener);
	}
	catch (final ImgIOException e) {
		openError(source, e);
	}
	return imgPlus;
}
 

public static void main( String[] args )
{
	new ImageJ();
	
	// test blending
	ImgFactory< FloatType > f = new ArrayImgFactory< FloatType >();
	Img< FloatType > img = f.create( new int[] { 400, 400 }, new FloatType() ); 
	
	Cursor< FloatType > c = img.localizingCursor();
	final int numDimensions = img.numDimensions();
	final double[] tmp = new double[ numDimensions ];
	
	// for blending
	final long[] dimensions = new long[ numDimensions ];
	img.dimensions( dimensions );
	final float percentScaling = 0.2f;
	final double[] border = new double[ numDimensions ];
				
	while ( c.hasNext() )
	{
		c.fwd();
		
		for ( int d = 0; d < numDimensions; ++d )
			tmp[ d ] = c.getFloatPosition( d );
		
		c.get().set( (float)BlendingPixelFusion.computeWeight( tmp, dimensions, border, percentScaling ) );
	}
	
	ImageJFunctions.show( img );
	Log.debug( "done" );
}
 

/**
 * @deprecated Use {@link #openAll(String, ImgFactory)}.
 */
@Deprecated
public static <T extends RealType<T> & NativeType<T>> List<SCIFIOImgPlus<T>>
	openImgs(final String source, final ImgFactory<T> imgFactory)
{
	return openAll(source, imgFactory);
}
 

@OpMethod(op = net.imagej.ops.create.imgFactory.DefaultCreateImgFactory.class)
public <T extends NativeType<T>> ImgFactory<T> imgFactory() {
	// NB: The generic typing of ImgFactory is broken; see:
	// https://github.com/imglib/imglib2/issues/91
	@SuppressWarnings("unchecked")
	final ImgFactory<T> result = (ImgFactory<T>) ops().run(
		Ops.Create.ImgFactory.class);
	return result;
}
 

@OpMethod(
	ops = net.imagej.ops.create.imgFactory.DefaultCreateImgFactory.class)
public <T extends NativeType<T>> ImgFactory<T> imgFactory(
	final Dimensions dims)
{
	// NB: The generic typing of ImgFactory is broken; see:
	// https://github.com/imglib/imglib2/issues/91
	@SuppressWarnings("unchecked")
	final ImgFactory<T> result = (ImgFactory<T>) ops().run(
		Ops.Create.ImgFactory.class, dims);
	return result;
}
 

/**
 * @see #open(Location, ImgFactory, SCIFIOConfig)
 */
public static <T extends RealType<T> & NativeType<T>> SCIFIOImgPlus<T>
open(final String source, final ImgFactory<T> imgFactory,
     final SCIFIOConfig config)
{
	return first(openAll(source, imgFactory, config));
}
 

@OpMethod(
	op = net.imagej.ops.create.imgLabeling.DefaultCreateImgLabeling.class)
public <L, T extends IntegerType<T>> ImgLabeling<L, T> imgLabeling(
	final Dimensions dims, final T outType, final ImgFactory<T> fac)
{
	@SuppressWarnings("unchecked")
	final ImgLabeling<L, T> result =
		(ImgLabeling<L, T>) ops().run(
			Ops.Create.ImgLabeling.class, dims,
			outType, fac);
	return result;
}
 

/**
 * @deprecated Use {@link #openAll(String, ImgFactory)}.
 */
@Deprecated
public static <T extends RealType<T> & NativeType<T>> List<SCIFIOImgPlus<T>>
	openImgs(final String source, final ImgFactory<T> imgFactory, final T type)
{
	return openAll(source, imgFactory);
}
 

@OpMethod(
	op = net.imagej.ops.create.imgLabeling.CreateImgLabelingFromInterval.class)
public <L, T extends IntegerType<T>> ImgLabeling<L, T> imgLabeling(
	final Interval interval, final T outType, final ImgFactory<T> fac)
{
	@SuppressWarnings("unchecked")
	final ImgLabeling<L, T> result =
		(ImgLabeling<L, T>) ops().run(
			Ops.Create.ImgLabeling.class,
			interval, outType, fac);
	return result;
}
 

@OpMethod(
	op = net.imagej.ops.create.imgLabeling.CreateImgLabelingFromInterval.class)
public <L, T extends IntegerType<T>> ImgLabeling<L, T> imgLabeling(
	final Interval interval, final T outType, final ImgFactory<T> fac,
	final int maxNumLabelSets)
{
	@SuppressWarnings("unchecked")
	final ImgLabeling<L, T> result =
		(ImgLabeling<L, T>) ops().run(
			Ops.Create.ImgLabeling.class,
			interval, outType, fac, maxNumLabelSets);
	return result;
}
 

/**
 * @deprecated Use {@link #open(String, ImgFactory, SCIFIOConfig)}.
 */
@Deprecated
public static <T extends RealType<T> & NativeType<T>> SCIFIOImgPlus<T>
	openImg(final String source, final ImgFactory<T> imgFactory,
		final SCIFIOConfig config) throws ImgIOException
{
	return open(source, imgFactory, config);
}
 

public static void display( final AbstractSpimData< ? > spimData, final ViewId viewId, final int pixelType, final String name, final boolean virtual )
{
	final ImgLoader imgLoader = (ImgLoader)spimData.getSequenceDescription().getImgLoader();
	final ImgFactory< ? extends NativeType< ? > > factory;
	final AbstractImgFactoryImgLoader il;

	// load as ImagePlus directly if possible
	if ( AbstractImgFactoryImgLoader.class.isInstance( imgLoader ) )
	{
		il = (AbstractImgFactoryImgLoader)imgLoader;
		factory = il.getImgFactory();
		il.setImgFactory( new ImagePlusImgFactory< FloatType >());
	}
	else
	{
		il = null;
		factory = null;
	}

	// display it
	DisplayImage export = new DisplayImage( virtual );
	
	if ( pixelType == 0 )
		export.exportImage( ((ImgLoader)spimData.getSequenceDescription().getImgLoader()).getSetupImgLoader( viewId.getViewSetupId() ).getFloatImage( viewId.getTimePointId(), false ), name );
	else
	{
		@SuppressWarnings( "unchecked" )
		RandomAccessibleInterval< UnsignedShortType > img =
			( RandomAccessibleInterval< UnsignedShortType > ) ((ImgLoader)spimData.getSequenceDescription().getImgLoader()).getSetupImgLoader( viewId.getViewSetupId() ).getImage( viewId.getTimePointId() );
		export.exportImage( img, name );
	}

	if ( factory != null && il != null )
		il.setImgFactory( factory );
}
 

private static <T extends RealType<T> & NativeType<T>> List<SCIFIOImgPlus<T>>
openAll(final ImgOpener opener, final Location source, final ImgFactory<T> imgFactory)
{
	List<SCIFIOImgPlus<T>> imgPlus = null;
	try {
		imgPlus = opener.openImgs(source, imgFactory);
		register(imgPlus, opener);
	}
	catch (final ImgIOException e) {
		openError(source, e);
	}
	return imgPlus;
}
 

protected ImgFactory< ? extends NativeType< ? > > selectImgFactory( final LightSheetZ1MetaData meta )
{
	long maxNumPixels = 0;

	for ( int a = 0; a < meta.numAngles(); ++a )
	{
		final int[] dim = meta.imageSizes().get( a );

		long n = 1;

		for ( int d = 0; d < dim.length; ++d )
		{
			n *= (long)dim[ d ];

			if ( dim[ d ] <= 0 )
			{
				IOFunctions.println( "Dimensions couldn't be read from metadata, using CellImg(256)." );
				return new CellImgFactory< FloatType >( 256 );
			}
		}

		maxNumPixels = Math.max( n, maxNumPixels );
	}

	int smallerLog2 = (int)Math.ceil( Math.log( maxNumPixels ) / Math.log( 2 ) );

	String s = "Maximum number of pixels in any view: n=" + maxNumPixels + 
			" (2^" + (smallerLog2-1) + " < n < 2^" + smallerLog2 + " px), ";

	if ( smallerLog2 <= 31 )
	{
		IOFunctions.println( s + "using ArrayImg." );
		return new ArrayImgFactory< FloatType >();
	}
	else
	{
		IOFunctions.println( s + "using CellImg(256)." );
		return new CellImgFactory< FloatType >( 256 );
	}
}
 

@Override
protected StackImgLoader createAndInitImgLoader( final String path, final File basePath, final ImgFactory< ? extends NativeType< ? > > imgFactory, SequenceDescription sequenceDescription )
{
	return new StackImgLoaderIJ(
			new File( basePath.getAbsolutePath(), path ),
			fileNamePattern, imgFactory,
			hasMultipleTimePoints, hasMultipleChannels, hasMultipleIlluminations, hasMultipleAngles,
			sequenceDescription );
}
 

@Override
protected StackImgLoader< ? > createAndInitImgLoader( final String path, final File basePath, final ImgFactory< ? extends NativeType< ? > > imgFactory, SequenceDescription sequenceDescription )
{
	return new StackImgLoaderLOCI(
			new File( basePath.getAbsolutePath(), path ),
			fileNamePattern, imgFactory,
			hasMultipleTimePoints, hasMultipleChannels, hasMultipleIlluminations, hasMultipleAngles,
			sequenceDescription );
}