mpicbg.imglib.image.ImageFactory Java Examples

final protected static Thread getCPUThread1( final AtomicInteger ai, final Block[] blocks, final int[] blockSize, final ImageFactory< FloatType > factory,
		final Image<FloatType> image, final Image<FloatType> result, final FourierConvolution<FloatType, FloatType> fftConvolution1 )
{
	final Thread cpuThread1 = new Thread(new Runnable()
	{
		public void run()
		{
			final Image< FloatType > block = factory.createImage( blockSize );

			int i;

			while ( ( i = ai.getAndIncrement() ) < blocks.length )
			{
				convolve1BlockCPU( blocks[ i ], i, image, result, block, fftConvolution1 );					
			}
			
			block.close();
		}
	});
	
	return cpuThread1;
}
 

final protected static Thread getCPUThread2( final AtomicInteger ai, final Block[] blocks, final int[] blockSize, final ImageFactory< FloatType > factory,
		final Image<FloatType> image, final Image<FloatType> result, final FourierConvolution<FloatType, FloatType> fftConvolution2 )
{
	final Thread cpuThread2 = new Thread(new Runnable()
	{
		public void run()
		{
			final Image< FloatType > block = factory.createImage( blockSize );

			int i;

			while ( ( i = ai.getAndIncrement() ) < blocks.length )
			{
				convolve2BlockCPU( blocks[ i ], image, result, block, fftConvolution2 );					
			}
			
			block.close();
		}
	});
	
	return cpuThread2;
}
 

final protected static Thread getCUDAThread1( final AtomicInteger ai, final Block[] blocks, final int[] blockSize, final ImageFactory< FloatType > factory,
		final Image<FloatType> image, final Image<FloatType> result, final int deviceId, final Image<FloatType> kernel1 )
{
	final Thread cudaThread1 = new Thread(new Runnable()
	{
		public void run()
		{
			final Image< FloatType > block = factory.createImage( blockSize );

			int i;

			while ( ( i = ai.getAndIncrement() ) < blocks.length )
			{
				convolve1BlockCUDA( blocks[ i ], i, deviceId, image, result, block, kernel1, blockSize );					
			}
			
			block.close();
		}
	});
	
	return cudaThread1;
}
 

final protected static Thread getCUDAThread2( final AtomicInteger ai, final Block[] blocks, final int[] blockSize, final ImageFactory< FloatType > factory,
		final Image<FloatType> image, final Image<FloatType> result, final int deviceId, final Image<FloatType> kernel2 )
{
	final Thread cudaThread2 = new Thread(new Runnable()
	{
		public void run()
		{
			final Image< FloatType > block = factory.createImage( blockSize );

			int i;

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

public MappingFusionParalell( final ViewStructure viewStructure, final ViewStructure referenceViewStructure,
							  final ArrayList<IsolatedPixelWeightenerFactory<?>> isolatedWeightenerFactories,
							  final ArrayList<CombinedPixelWeightenerFactory<?>> combinedWeightenerFactories )
{
	super( viewStructure, referenceViewStructure, isolatedWeightenerFactories, combinedWeightenerFactories );


	if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_MAIN )
		IOFunctions.println("(" + new Date(System.currentTimeMillis()) + "): Reserving memory for fused image.");

	final ImageFactory<FloatType> fusedImageFactory = new ImageFactory<FloatType>( new FloatType(), conf.processImageFactory );
	fusedImage = fusedImageFactory.createImage( new int[]{ imgW, imgH, imgD }, "Fused image");

	if (fusedImage == null)
	{
		if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_ERRORONLY )
			IOFunctions.println("MappingFusionParalell.constructor: Cannot create output image: " + conf.processImageFactory.getErrorMessage());

		return;
	}
}
 

public MappingFusionParalellMaxWeight( final ViewStructure viewStructure, final ViewStructure referenceViewStructure,
							  final ArrayList<IsolatedPixelWeightenerFactory<?>> isolatedWeightenerFactories,
							  final ArrayList<CombinedPixelWeightenerFactory<?>> combinedWeightenerFactories )
{
	super( viewStructure, referenceViewStructure, isolatedWeightenerFactories, combinedWeightenerFactories );


	if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_MAIN )
		IOFunctions.println("(" + new Date(System.currentTimeMillis()) + "): Reserving memory for fused image.");

	final ImageFactory<FloatType> fusedImageFactory = new ImageFactory<FloatType>( new FloatType(), conf.processImageFactory );
	fusedImage = fusedImageFactory.createImage( new int[]{ imgW, imgH, imgD }, "Fused image");

	if (fusedImage == null)
	{
		if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_ERRORONLY )
			IOFunctions.println("MappingFusionParalell.constructor: Cannot create output image: " + conf.processImageFactory.getErrorMessage());

		return;
	}
}
 

IDT2D(final Image<BitType> img) {
	this.w = img.getDimension(0);
	this.h = img.getDimension(1);
	ImageFactory<IntType> f = new ImageFactory<IntType>(new IntType(), new ArrayContainerFactory());
	this.result = f.createImage(new int[]{w, h});

	// Set all result pixels to infinity
	final int infinity = (w + h) * 9;
	for (final IntType v : this.result) {
		v.set(infinity);
	}

	// init result pixels with those of the image:
	this.csrc = img.createLocalizableByDimCursor();
	this.cout = result.createLocalizableByDimCursor();

	int count = 0;
	for (int y = 0; y < h; y++) {
		for (int x = 0; x < w; x++) {
			if (isBoundary(x, y)) {
				setOutValueAt(x, y, 0);
				count++;
			} else if (isJustOutside(x, y)) {
				setOutValueAt(x, y, -1);
			}
		}
	}

	if (count > 0) {
		propagate();
	}

	csrc.close();
	cout.close();
}
 

public Image<FloatType> getFoundBeads( final ViewDataBeads view )
{
	// display found beads
	ImageFactory<FloatType> factory = new ImageFactory<FloatType>( new FloatType(), view.getViewStructure().getSPIMConfiguration().inputImageFactory );
	Image<FloatType> img = factory.createImage( view.getImageSize() );
	
	LocalizableByDimCursor3D<FloatType> cursor = (LocalizableByDimCursor3D<FloatType>) img.createLocalizableByDimCursor();		

   	final float[] tmp = new float[ img.getNumDimensions() ];

	for ( Bead bead : view.getBeadStructure().getBeadList())
	{
       	for ( int d = 0; d < tmp.length; ++d )
       		tmp[ d ] = (float)bead.getL()[ d ];

		final LocalizablePoint p = new LocalizablePoint( tmp );
		
		HyperSphereIterator<FloatType> it = new HyperSphereIterator<FloatType>( img, p, 1, new OutOfBoundsStrategyValueFactory<FloatType>() );
		
		for ( final FloatType f : it )
			f.setOne();			
	}
	
	cursor.close();

	return img;
}
 

public MappingFusionSequential( final ViewStructure viewStructure, final ViewStructure referenceViewStructure,
		  						final ArrayList<IsolatedPixelWeightenerFactory<?>> isolatedWeightenerFactories,
		  						final ArrayList<CombinedPixelWeightenerFactory<?>> combinedWeightenerFactories,
		  						final int numParalellStacks )
{
	super( viewStructure, referenceViewStructure, isolatedWeightenerFactories, combinedWeightenerFactories );

	if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_MAIN )
		IOFunctions.println("(" + new Date(System.currentTimeMillis()) + "): Reserving memory for fused image.");

	final ImageFactory<FloatType> fusedImageFactory = new ImageFactory<FloatType>( new FloatType(), conf.processImageFactory );

	fusedImage = fusedImageFactory.createImage( new int[]{ imgW, imgH, imgD }, "Fused image");
	weights = fusedImageFactory.createImage( new int[]{ imgW, imgH, imgD }, "Weights image");
	this.numParalellStacks = numParalellStacks;

	if ( fusedImage == null )
	{
		if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_ERRORONLY )
			IOFunctions.println("MappingFusionSequentialImages.constructor: Cannot create output image: " + conf.processImageFactory.getErrorMessage());

		if ( weights != null )
			weights.close();

		return;
	}

	if ( weights == null )
	{
		if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_ERRORONLY )
			IOFunctions.println("MappingFusionSequentialImages.constructor: Cannot create weights image: " + conf.processImageFactory.getErrorMessage());

		fusedImage.close();

		return;
	}
}
 

public PreDeconvolutionFusionSequential( final ViewStructure viewStructure, final ViewStructure referenceViewStructure, 
							  final ArrayList<IsolatedPixelWeightenerFactory<?>> isolatedWeightenerFactories, 
							  final ArrayList<CombinedPixelWeightenerFactory<?>> combinedWeightenerFactories )
{
	super( viewStructure, referenceViewStructure, isolatedWeightenerFactories, combinedWeightenerFactories );				
	
	// normalize the weights so the the sum for each pixel over all views is 1?
	this.normalize = true;
	
	if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_MAIN )
		IOFunctions.println("(" + new Date(System.currentTimeMillis()) + "): Reserving memory for fused image.");
	
	final ImageFactory<FloatType> imageFactory = new ImageFactory<FloatType>( new FloatType(), conf.processImageFactory );
	numViews = viewStructure.getNumViews();
	
	if ( conf.deconvolutionJustShowOverlap )
	{
		overlap = imageFactory.createImage( new int[]{ imgW, imgH, imgD }, "overlap" );
		images = null;
		weights = null;
		extractPSF = null;
	}
	else
	{
		overlap = null;
		
		if ( conf.extractPSF )
			extractPSF = new ExtractPSF( viewStructure );
		else
			extractPSF = ExtractPSF.loadAndTransformPSF( conf.psfFiles, conf.transformPSFs, viewStructure );
		
		images = new Image[ numViews ];
		weights = new Image[ numViews ];

		if ( extractPSF == null )
			return;
		
		for ( int view = 0; view < numViews; view++ )
		{
			weights[ view ] = imageFactory.createImage( new int[]{ imgW, imgH, imgD }, "weights_" + view );
			images[ view ] = imageFactory.createImage( new int[]{ imgW, imgH, imgD }, "view_" + view ); 
			
			if ( images[ view ] == null || weights[ view ] == null )
			{
				if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_ERRORONLY )
					IOFunctions.println("PreDeconvolutionFusion.constructor: Cannot create output image: " + conf.processImageFactory.getErrorMessage() );

				return;
			}
		}
	}
}
 

public MappingFusionSequentialDifferentOutput( final ViewStructure viewStructure, final ViewStructure referenceViewStructure,
		  									   final ArrayList<IsolatedPixelWeightenerFactory<?>> isolatedWeightenerFactories,
		  									   final ArrayList<CombinedPixelWeightenerFactory<?>> combinedWeightenerFactories,
		  									   final int numParalellStacks )
{
	super( viewStructure, referenceViewStructure, isolatedWeightenerFactories, combinedWeightenerFactories );

	numViews = viewStructure.getNumViews();

	if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_MAIN )
		IOFunctions.println("(" + new Date(System.currentTimeMillis()) + "): Reserving memory for fused images.");

	if ( angleIndiciesStatic == null )
	{
		angleIndicies = new int[ numViews ];

		for ( int view = 0; view < numViews; view++ )
			angleIndicies[ view ] = view;
	}
	else
	{
		if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_MAIN )
			IOFunctions.println("(" + new Date(System.currentTimeMillis()) + "): WARNING: Using statically defined angle-indices from class mpicbg.spim.fusion.MappingFusionSequentialDifferentOutput: " + Util.printCoordinates( angleIndicies ) );
		
		angleIndicies = angleIndiciesStatic.clone();
	}

	this.numParalellStacks = numParalellStacks;

	IJ.log( "nump = " + numParalellStacks );

	fusedImages = new Image[ angleIndicies.length ];
	final ImageFactory<FloatType> fusedImageFactory = new ImageFactory<FloatType>( new FloatType(), conf.processImageFactory );

	final long size = (4l * imgW * imgH * imgD)/(1000l*1000l);

	for (int i = 0; i < angleIndicies.length; i++)
	{
		if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_MAIN )
			IOFunctions.println( "(" + new Date(System.currentTimeMillis()) + "): Reserving " + size + " MiB for '" + viewStructure.getViews().get( angleIndicies[i] ).getName() + "'" );

		fusedImages[ i ] = fusedImageFactory.createImage( new int[]{ imgW, imgH, imgD }, "Fused image" );

		if ( fusedImages[i] == null && viewStructure.getDebugLevel() <= ViewStructure.DEBUG_ERRORONLY )
			IOFunctions.printErr("MappingFusionSequentialDifferentOutput.constructor: Cannot create output image: " + conf.processImageFactory.getErrorMessage());
	}
}
 

public static Image<FloatType> computeEntropy(final Image<FloatType> img, final ContainerFactory entropyType, final int histogramBins, final int windowSizeX, final int windowSizeY, final int windowSizeZ)
{
	// check if we can use fast forward algorithm		
	if ( Array3D.class.isInstance( img.getContainer() ) )
	{
		IOFunctions.println("Input is instance of Image<Float> using an Array3D, fast forward algorithm --- Fast Forward Algorithm available.");
		return EntropyFloatArray3D.computeEntropy( img, entropyType, histogramBins, windowSizeX, windowSizeY, windowSizeZ); 
	}	
	
	final float maxEntropy = getMaxEntropy(histogramBins);
	final ImageFactory<FloatType> factory = new ImageFactory<FloatType>( new FloatType(), entropyType );
	final Image<FloatType> entropy = factory.createImage( img.getDimensions(), "Entropy of " + img.getName() );
	
	final LocalizableByDimCursor<FloatType> entropyIterator = entropy.createLocalizableByDimCursor( new OutOfBoundsStrategyMirrorFactory<FloatType>() );
	
	final Entropy ei = Entropy.initEntropy( img, histogramBins, windowSizeX, windowSizeY, windowSizeZ, 0, 0, 0);
	
	final int directionZ = +1; 
	int directionY = +1;
	int directionX = +1;	
	
	final int width = img.getDimension( 0 );
	final int height = img.getDimension( 1 );
	final int depth = img.getDimension( 2 );
	
	for (int z = 0; z < depth; z++)
	{    			
		for (int y = 0; y < height; y++)
   		{
   			for (int x = 0; x < width; x++)
   			{
   				if (x != 0)
   					ei.updateEntropyX(directionX);
   				    	
   				entropyIterator.move( ei.getX() - entropyIterator.getPosition(0), 0 );
   				entropyIterator.move( ei.getY() - entropyIterator.getPosition(1), 1 );
   				entropyIterator.move( ei.getZ() - entropyIterator.getPosition(2), 2 );
   				
   				entropyIterator.getType().set( ei.getEntropy() / maxEntropy );
			}    			
   			directionX *= -1;
   			
   			if (y != height - 1)
   				ei.updateEntropyY(directionY);
   		}                            		
   		directionY *= -1;
   		
		if (z != depth - 1)
			ei.updateEntropyZ(directionZ);                            		
	}
			
	entropyIterator.close();
	ei.close();
	
	return entropy;
}
 

public static Image<FloatType> computeEntropy(final Image<FloatType> image, final ContainerFactory entropyType, final int histogramBins, final int windowSizeX, final int windowSizeY, final int windowSizeZ)
{
	final float maxEntropy = getMaxEntropy(histogramBins);
	
	final ImageFactory<FloatType> factory = new ImageFactory<FloatType>( new FloatType(), entropyType );
	final Image<FloatType> entropy = factory.createImage( image.getDimensions(), "Entropy of " + image.getName() );

	final LocalizableByDimCursor3D<FloatType> it = (LocalizableByDimCursor3D<FloatType>)entropy.createLocalizableByDimCursor();
	
	final EntropyFloatArray3D entropyObject = EntropyFloatArray3D.initEntropy( image, histogramBins, windowSizeX, windowSizeY, windowSizeZ, 0, 0, 0 );
	
	final int directionZ = +1; 
	int directionY = +1;
	int directionX = +1;
	
	for (int z = 0; z < image.getDimension( 2 ); z++)
	{
		for (int y = 0; y < image.getDimension( 1 ); y++)
   		{
   			for (int x = 0; x < image.getDimension( 0 ); x++)
   			{
   				if (x != 0)
   					entropyObject.updateEntropyX(directionX);
   				    				
   				it.setPosition( entropyObject.getX(), entropyObject.getY(), entropyObject.getZ() );
   				it.getType().set( entropyObject.getEntropy() / maxEntropy );
			}    			
   			directionX *= -1;
   			
   			if (y != image.getDimension( 1 ) - 1)
   				entropyObject.updateEntropyY(directionY);
   		}                            		
   		directionY *= -1;
   		
		if (z != image.getDimension( 2 ) - 1)
			entropyObject.updateEntropyZ(directionZ);                            		
	}
	
	entropyObject.cIn.close();
	entropyObject.cOut.close();
	
	return entropy;
}
 

public PreDeconvolutionFusion( final ViewStructure viewStructure, final ViewStructure referenceViewStructure, 
							  final ArrayList<IsolatedPixelWeightenerFactory<?>> isolatedWeightenerFactories, 
							  final ArrayList<CombinedPixelWeightenerFactory<?>> combinedWeightenerFactories )
{
	super( viewStructure, referenceViewStructure, isolatedWeightenerFactories, combinedWeightenerFactories );				
	
	// normalize the weights so the the sum for each pixel over all views is 1?
	this.normalize = true;

	if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_MAIN )
		IOFunctions.println("(" + new Date(System.currentTimeMillis()) + "): Reserving memory for fused image.");
	
	final ImageFactory<FloatType> imageFactory = new ImageFactory<FloatType>( new FloatType(), conf.processImageFactory );
	numViews = viewStructure.getNumViews();
			
	if ( conf.deconvolutionJustShowOverlap )
	{
		overlap = imageFactory.createImage( new int[]{ imgW, imgH, imgD }, "overlap" );
		images = null;
		weights = null;
		extractPSF = null;
	}
	else
	{
		overlap = null;
		
		if ( conf.extractPSF )
			extractPSF = new ExtractPSF( viewStructure );
		else
			extractPSF = ExtractPSF.loadAndTransformPSF( conf.psfFiles, conf.transformPSFs, viewStructure );			
		
		images = new Image[ numViews ];
		weights = new Image[ numViews ];

		if ( extractPSF == null )
			return;

		for ( int view = 0; view < numViews; view++ )
		{
			weights[ view ] = imageFactory.createImage( new int[]{ imgW, imgH, imgD }, "weights_" + view );
			images[ view ] = imageFactory.createImage( new int[]{ imgW, imgH, imgD }, "view_" + view ); 
			
			if ( images[ view ] == null || weights[ view ] == null )
			{
				if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_ERRORONLY )
					IOFunctions.println("PreDeconvolutionFusion.constructor: Cannot create output image: " + conf.processImageFactory.getErrorMessage() );

				return;
			}
		}	
	}
}
 

/**
 * return an {@code Image<T>} as input for the PhaseCorrelation.
 * 
 * @param imp - the {@link ImagePlus}
 * @param imgFactory - the {@link ImageFactory} defining wher to put it into
 * @param channel - which channel (if channel=0 means average all channels)
 * @param timepoint - which timepoint
 * 
 * @return - the {@link Image} or null if it was not an ImagePlus.GRAY8, ImagePlus.GRAY16 or ImagePlus.GRAY32
 */
public static < T extends RealType<T> > Image<T> getImage( final ImagePlus imp, Roi roi, final ImageFactory<T> imgFactory, final int channel, final int timepoint )
{
	// first test the roi
	roi = getOnlyRectangularRoi( roi );
	
	// how many dimensions?
	final int numDimensions;		
	if ( imp.getNSlices() > 1 )
		numDimensions = 3;
	else
		numDimensions = 2;
	
	// the size of the image
	final int[] size = new int[ numDimensions ];
	final int[] offset = new int[ numDimensions ];
	
	if ( roi == null )
	{
		size[ 0 ] = imp.getWidth();
		size[ 1 ] = imp.getHeight();
		
		if ( numDimensions == 3 )
			size[ 2 ] = imp.getNSlices();
	}
	else
	{
		size[ 0 ] = roi.getBounds().width;
		size[ 1 ] = roi.getBounds().height;

		offset[ 0 ] = roi.getBounds().x;
		offset[ 1 ] = roi.getBounds().y;
		
		if ( numDimensions == 3 )
			size[ 2 ] = imp.getNSlices();
	}
	
	// create the Image
	final Image<T> img = imgFactory.createImage( size );
	final boolean success;
	
	// copy the content
	if ( channel == 0 )
	{
		// we need to average all channels
		success = averageAllChannels( img, offset, imp, timepoint );
	}
	else
	{
		// otherwise only copy one channel
		success = fillInChannel( img, offset, imp, channel, timepoint );
	}
	
	if ( success )
	{
		return img;
	}
	img.close();
	return null;
}
 

/**
 * Extract the current 2d region of interest from the souce image
 * 
 * @param source - the source image, a {@link Image} which is a copy of the {@link ImagePlus}
 * @param rectangle - the area of interest
 * @param extraSize - the extra size around so that detections at the border of the roi are not messed up
 * @return
 */
protected Image<FloatType> extractImage( final FloatImagePlus< net.imglib2.type.numeric.real.FloatType > source, final Rectangle rectangle, final int extraSize )
{
	final Image<FloatType> img = new ImageFactory<FloatType>( new FloatType(), new ArrayContainerFactory() ).createImage( new int[]{ rectangle.width+extraSize, rectangle.height+extraSize } );
	
	final int offsetX = rectangle.x - extraSize/2;
	final int offsetY = rectangle.y - extraSize/2;

	final int[] location = new int[ source.numDimensions() ];
	
	if ( location.length > 2 )
		location[ 2 ] = (imp.getCurrentSlice()-1)/imp.getNChannels();
			
	final LocalizableCursor<FloatType> cursor = img.createLocalizableCursor();
	final RandomAccess<net.imglib2.type.numeric.real.FloatType> positionable;
	
	if ( offsetX >= 0 && offsetY >= 0 && 
		 offsetX + img.getDimension( 0 ) < source.dimension( 0 ) && 
		 offsetY + img.getDimension( 1 ) < source.dimension( 1 ) )
	{
		// it is completely inside so we need no outofbounds for copying
		positionable = source.randomAccess();
	}
	else
	{
		positionable = Views.extendMirrorSingle( source ).randomAccess();
	}
		
	while ( cursor.hasNext() )
	{
		cursor.fwd();
		cursor.getPosition( location );
		
		location[ 0 ] += offsetX;
		location[ 1 ] += offsetY;
		
		positionable.setPosition( location );
		
		cursor.getType().set( positionable.get().get() );
	}
	
	return img;
}
 

/**
 * Extracts the PSF by averaging the local neighborhood RANSAC correspondences
 * @param view - the SPIM view
 * @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 Image<FloatType> extractPSF( final ViewDataBeads view, final int[] size )
{
	final int numDimensions = size.length;
	
	// Mirror produces some artifacts ...
	final OutOfBoundsStrategyFactory<FloatType> outside = new OutOfBoundsStrategyPeriodicFactory<FloatType>();
	final InterpolatorFactory<FloatType> interpolatorFactory = new LinearInterpolatorFactory<FloatType>( outside );
	
	final ImageFactory<FloatType> imageFactory = new ImageFactory<FloatType>( new FloatType(), view.getViewStructure().getSPIMConfiguration().processImageFactory );
	final Image<FloatType> img = view.getImage();
	final Image<FloatType> psf = imageFactory.createImage( size );
	
	final Interpolator<FloatType> interpolator = img.createInterpolator( interpolatorFactory );
	final LocalizableCursor<FloatType> psfCursor = psf.createLocalizableCursor();
	
	final int[] sizeHalf = size.clone();		
	for ( int d = 0; d < numDimensions; ++d )
		sizeHalf[ d ] /= 2;
	
	int numRANSACBeads = 0;
	
	for ( final Bead bead : view.getBeadStructure().getBeadList() )
	{			
		final double[] position = bead.getL().clone();
		final int[] tmpI = new int[ position.length ];
		final double[] tmpF = new double[ position.length ];
		
		// check if it is a true correspondence
		if ( bead.getRANSACCorrespondence().size() > 0 ) 
		{
			++numRANSACBeads;
			psfCursor.reset();
			
			while ( psfCursor.hasNext() )
			{
				psfCursor.fwd();
				psfCursor.getPosition( tmpI );

				for ( int d = 0; d < numDimensions; ++d )
					tmpF[ d ] = tmpI[ d ] - sizeHalf[ d ] + position[ d ];
				
				interpolator.moveTo( tmpF );
				
				psfCursor.getType().add( interpolator.getType() );
			}
		}
	}

	// compute the average		
	final FloatType n = new FloatType( numRANSACBeads );
	
	psfCursor.reset();
	while ( psfCursor.hasNext() )
	{
		psfCursor.fwd();
		psfCursor.getType().div( n );			
	}	

	ViewDataBeads.normalizeImage( psf );
	
	// TODO: Remove
	//ImageJFunctions.show( psf );
	
	return psf;
}
 

protected static Image< FloatType > loadInitialImage( final String fileName, final boolean checkNumbers, final float minValue, final int[] dimensions, final ImageFactory< FloatType > imageFactory )
{
	IOFunctions.println( "Loading image '" + fileName + "' as start for iteration." );
	Image< FloatType > psi = LOCI.openLOCIFloatType( fileName, imageFactory );
	
	if ( psi == null )
	{
		IOFunctions.println( "Could not load image '" + fileName + "'." );
		return null;
	}
	else
	{
		boolean dimensionsMatch = true;
		
		for ( int d = 0; d < psi.getNumDimensions(); ++d )
			if ( psi.getDimension( d ) != dimensions[ d ] )
				dimensionsMatch = false;
		
		if ( !dimensionsMatch )
		{
			IOFunctions.println( "Dimensions of '" + fileName + "' do not match: " + Util.printCoordinates( psi.getDimensions() ) + " != " + Util.printCoordinates( dimensions ) );
			psi.close();
			return null;
		}
		
		if ( checkNumbers )
		{
			IOFunctions.println( "Checking values of '" + fileName + "' you can disable this check by setting mpicbg.spim.postprocessing.deconvolution2.BayesMVDeconvolution.checkNumbers = false;" );
			
			boolean smaller = false;
			boolean hasZerosOrNeg = false;
			
			for ( final FloatType v : psi )
			{
				if ( v.get() < minValue )
					smaller = true;

				if ( v.get() <= 0 )
				{
					hasZerosOrNeg = true;
					v.set( minValue );
				}
			}

			if ( smaller )
				IOFunctions.println( "Some values '" + fileName + "' are smaller than the minimal value of " + minValue + ", this can lead to instabilities." );

			if ( hasZerosOrNeg )
				IOFunctions.println( "Some values '" + fileName + "' were smaller or equal to zero, they have been replaced with the min value of " + minValue );
		}
	}
	
	return psi;
}
 

public LRFFT(
		final Image<FloatType> image,
		final Image<FloatType> weight,
		final Image<FloatType> kernel,
		final int[] deviceList, final boolean useBlocks, final int[] blockSize )
{
	this.image = image;
	this.kernel1 = kernel;
	this.weight = weight;
	
	this.deviceList = deviceList;
	this.device0 = deviceList[ 0 ];
	this.numDevices = deviceList.length;

	// figure out if we need GPU and/or CPU
	boolean anyGPU = false;
	boolean anyCPU = false;
	
	for ( final int i : deviceList )
	{
		if ( i >= 0 )
			anyGPU = true;
		else if ( i == -1 )
			anyCPU = true;
	}
	
	this.useCUDA = anyGPU;
	this.useCPU = anyCPU;
			
	if ( useBlocks )
	{
		this.useBlocks = true;
		
		// define the blocksize so that it is one single block
		this.blockSize = new int[ image.getNumDimensions() ];
					
		for ( int d = 0; d < this.blockSize.length; ++d )
			this.blockSize[ d ] = blockSize[ d ];
					
		this.blocks = Block.divideIntoBlocks( image.getDimensions(), this.blockSize, kernel.getDimensions() );
		
		// blocksize might change during division if they were too small
		 //this.blockSize = blockSize.clone();
		
		IOFunctions.println( "Number of blocks: " + this.blocks.length );
		
		this.factory = new ImageFactory< FloatType >( new FloatType(), new ArrayContainerFactory() );
	}
	else if ( this.useCUDA ) // and no blocks, i.e. one big block
	{
		this.useBlocks = true;
		
		// define the blocksize so that it is one single block
		this.blockSize = new int[ image.getNumDimensions() ];
		
		for ( int d = 0; d < this.blockSize.length; ++d )
			this.blockSize[ d ] = image.getDimension( d ) + kernel.getDimension( d ) - 1;
		
		this.blocks = Block.divideIntoBlocks( image.getDimensions(), this.blockSize, kernel.getDimensions() );
		this.factory = new ImageFactory< FloatType >( new FloatType(), new ArrayContainerFactory() );			
	}
	else
	{
		this.blocks = null;
		this.blockSize = null;
		this.factory = null;
		this.useBlocks = false;
	}		
}
 

public static void main( String[] args )
{
	new ImageJ();
	//final Image< FloatType > img = LOCI.openLOCIFloatType( "/Users/preibischs/Desktop/Geburtstagsfaust.jpg", new ArrayContainerFactory() );
	final Image< FloatType > img = LOCI.openLOCIFloatType( "/Users/preibischs/Desktop/spim.tif", new ArrayContainerFactory() );
			
	final Image< FloatType > kernel = FourierConvolution.createGaussianKernel( img.getContainerFactory(), new double[]{ 10, 10, 10 } );
	//ImageJFunctions.show( img );
	//final FourierConvolution< FloatType, FloatType > t = new FourierConvolution<FloatType, FloatType>( img, kernel );
	//t.process();
	//ImageJFunctions.show( t.getResult() );		
	
	final int[] imgSize = img.getDimensions();//new int[]{ 256, 384 };
	final int[] blockSize = new int[]{ 256, 256, 256 };
	
	//for ( int d = 0; d < blockSize.length; ++d )
	//	blockSize[ d ] = img.getDimension( d ) + kernel.getDimension( d ) - 1;
	
	final int[] kernelSize = kernel.getDimensions();//new int[]{ 51, 25 };
	
	final Block[] blocks = Block.divideIntoBlocks( imgSize, blockSize, kernelSize );
	
	final ArrayList< Image< FloatType > > blockImgs = new ArrayList< Image< FloatType > >();
	final ImageFactory< FloatType > factory = new ImageFactory< FloatType >( new FloatType(), new ArrayContainerFactory() );
	
	ImageJFunctions.show( img );	
	
	for ( int i = 0; i < blocks.length; ++i )
		blockImgs.add( factory.createImage( blockSize ) );
	
	long time = 0;
	
	//while ( time >= 0 )
	{
		time = System.currentTimeMillis();
		for ( int i = 0; i < blocks.length; ++i )
		{
			blocks[ i ].copyBlock( img, blockImgs.get( i ) );
			//ImageJFunctions.show( block );
		}

		System.out.println( System.currentTimeMillis() - time );
		//System.exit( 0 );
	}
	
	
	/*
	final FourierConvolution< FloatType, FloatType > t2 = new FourierConvolution<FloatType, FloatType>( blockImgs.get( 0 ), kernel );
	//t2.setExtendImageByKernelSize( false );

		//for ( final FloatType t : block )
		//	t.set( t.get() + 20*(i+1) );

	*/

	final Image< FloatType > img2 = img.createNewImage();
	
	//while ( time > 0 )
	{
		time = System.currentTimeMillis();

		for ( int i = 0; i < blocks.length; ++i )
		{
			//t2.replaceImage( blockImgs.get( i ) );
			//t2.process();
			
			blocks[ i ].pasteBlock( img2,  blockImgs.get( i ) );
			//ImageJFunctions.show( t.getResult() );
		}
		
		System.out.println( System.currentTimeMillis() - time );
	}
	ImageJFunctions.show( img2 );
	
}