mpicbg.imglib.image.Image Java Examples

final static public void computeMinMax( final Image< FloatType > img, final FloatType min, final FloatType max )
{
	min.set( Float.MAX_VALUE );
	max.set( Float.MIN_VALUE );
	
	for ( final FloatType t : img )
	{
		final float value = t.get();
		
		if ( value > max.get() )
			max.set( value );
		
		if ( value < min.get() )
			min.set( value );
	}
}
 

final protected static void convolve1BlockCUDA( final Block blockStruct, final int i, final int deviceId, final Image<FloatType> image, final Image<FloatType> result, final Image<FloatType> block, 
		final Image<FloatType> kernel1, final int[] blockSize )
{
	long time = System.currentTimeMillis();
	blockStruct.copyBlock( image, block );
	System.out.println( " block " + i + "(CPU  " + deviceId + "): copy " + (System.currentTimeMillis() - time) );

	// convolve block with kernel1 using CUDA
	time = System.currentTimeMillis();				
	LRFFT.cuda.convolution3DfftCUDAInPlace( ((FloatArray)((Array)block.getContainer()).update( null )).getCurrentStorageArray(), getCUDACoordinates( blockSize ), 
			((FloatArray)((Array)kernel1.getContainer()).update( null )).getCurrentStorageArray(), getCUDACoordinates( kernel1.getDimensions() ), deviceId );
	System.out.println( " block " + i + "(CUDA " + deviceId + "): compute " + (System.currentTimeMillis() - time) );

	time = System.currentTimeMillis();
	blockStruct.pasteBlock( result, block );
	System.out.println( " block " + i + "(CPU  " + deviceId + "): paste " + (System.currentTimeMillis() - time) );
}
 

@SuppressWarnings("rawtypes")
public static final Img< net.imglib2.type.numeric.real.FloatType > wrap( final Image< FloatType > i )
{
	final ContainerFactory c = i.getContainerFactory();
	
	if ( c instanceof ImagePlusContainerFactory )
	{
		try
		{
			return net.imglib2.img.display.imagej.ImageJFunctions.wrapFloat( ((ImagePlusContainer)i.getContainer()).getImagePlus() );
		}
		catch (mpicbg.imglib.exception.ImgLibException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
			return null;
		}			
	}
	else
	{
		return ImgLib1.wrapFloatToImgLib2( i );
	}
}
 

@SuppressWarnings("rawtypes")
public static final Image< FloatType > wrap( final Img< net.imglib2.type.numeric.real.FloatType > i )
{
	if ( i instanceof ImagePlusImg )
	{
		try
		{
			return ImageJFunctions.wrapFloat( ((ImagePlusImg) i).getImagePlus() );
		}
		catch (ImgLibException e)
		{
			// TODO Auto-generated catch block
			e.printStackTrace();
			return null;
		}			
	}
	else
	{
		return ImgLib2.wrapFloatToImgLib1( i );
	}
}
 

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 private static BigDecimal sumImage( final Image<FloatType> img )
{
	BigDecimal sum = new BigDecimal( 0, MathContext.UNLIMITED );
	
	final Cursor<FloatType> cursorImg = img.createCursor();
	
	while ( cursorImg.hasNext() )
	{
		cursorImg.fwd();
		sum = sum.add( BigDecimal.valueOf( (double)cursorImg.getType().get() ) );
	}
	
	cursorImg.close();

	return sum;
}
 

public ExtractPSF( final SPIMConfiguration config )
{
	//
	// load the files
	//
	final ViewStructure viewStructure = ViewStructure.initViewStructure( config, 0, new AffineModel3D(), "ViewStructure Timepoint 0", config.debugLevelInt );						

	for ( ViewDataBeads view : viewStructure.getViews() )
	{
		view.loadDimensions();
		view.loadSegmentation();
		view.loadRegistration();

		BeadRegistration.concatenateAxialScaling( view, ViewStructure.DEBUG_MAIN );
	}
	
	this.viewStructure = viewStructure;
	this.pointSpreadFunctions = new ArrayList<Image<FloatType>>();
	this.originalPSFs = new ArrayList<Image<FloatType>>();
	this.conf = config;
	
	setPSFSize( Multi_View_Deconvolution.psfSize, Multi_View_Deconvolution.isotropic, Multi_View_Deconvolution.psfSize3d );
}
 

private static final void computeQuotient( final long start, final long loopSize, final Image< FloatType > psiBlurred, final LRFFT processingData )
{
	final Cursor<FloatType> cursorImg = processingData.getImage().createCursor();
	final Cursor<FloatType> cursorPsiBlurred = psiBlurred.createCursor();
	
	cursorImg.fwd( start );
	cursorPsiBlurred.fwd( start );
	
	for ( long l = 0; l < loopSize; ++l )
	{
		cursorImg.fwd();
		cursorPsiBlurred.fwd();
		
		final float imgValue = cursorImg.getType().get();
		final float psiBlurredValue = cursorPsiBlurred.getType().get();
		
		cursorPsiBlurred.getType().set( imgValue / psiBlurredValue );
	}
	
	cursorImg.close();
	cursorPsiBlurred.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;
}
 

/**
 * Adds additive gaussian noise: i = i + gauss(x, sigma)
 * 
 * @return the signal-to-noise ratio (measured)
 */
public static double addGaussianNoise( final Image< FloatType > img, final Random rnd, final float sigma, boolean onlyPositive )
{
	for ( final FloatType f : img )
	{
		float newValue = f.get() + (float)( rnd.nextGaussian() * sigma );
		
		if ( onlyPositive )
			newValue = Math.max( 0, newValue );
		
		f.set( newValue );
	}
	
	return 1;
}
 

/**
 * Make image the same size as defined, center it
 * 
 * @param img
 * @return
 */
public static Image< FloatType > makeSameSize( final Image< FloatType > img, final int[] sizeIn )
{
	final int[] size = sizeIn.clone();

	float min = Float.MAX_VALUE;

	for ( final FloatType f : img )
		min = Math.min( min, f.get() );
	
	final Image< FloatType > square = img.createNewImage( size );
	
	final LocalizableCursor< FloatType > squareCursor = square.createLocalizableCursor();
	final LocalizableByDimCursor< FloatType > inputCursor = img.createLocalizableByDimCursor( new OutOfBoundsStrategyValueFactory<FloatType>( new FloatType( min ) ) );
	
	while ( squareCursor.hasNext() )
	{
		squareCursor.fwd();
		squareCursor.getPosition( size );
		
		for ( int d = 0; d < img.getNumDimensions(); ++d )
			size[ d ] =  size[ d ] - square.getDimension( d )/2 + img.getDimension( d )/2;

		inputCursor.setPosition( size );
		squareCursor.getType().set( inputCursor.getType().get() );
	}
	
	return square;
}
 

public static <T extends RealType<T>> ArrayList< DifferenceOfGaussianPeak<T> > extractBeadsLaPlaceImgLib( 
		final Image<T> img,
		final float initialSigma,
		float minPeakValue,
		float minInitialPeakValue )
{
	return extractBeadsLaPlaceImgLib(img, new OutOfBoundsStrategyMirrorFactory<T>(), 0.5f, initialSigma, minPeakValue, minInitialPeakValue, 4, true, false, ViewStructure.DEBUG_MAIN );
}
 

public static Image< FloatType > computeContentBasedWeighting( final Image< FloatType > img, final int fusionSigma1, final int fusionSigma2, final float zStretching )
{
	// compute the radii
	final int rxy1 = Math.round( fusionSigma1 );
	final int rxy2 = Math.round( fusionSigma2 );

	final int rz1 = Math.round( fusionSigma1 / zStretching );
	final int rz2 = Math.round( fusionSigma2 / zStretching );
				
	// compute I*sigma1, store in imgConv
	final Image< LongType > integralImg = IntegralImage3d.compute( img );
	final Image< FloatType > imgConv = img.createNewImage();			
	DOM.meanMirror( integralImg, imgConv, rxy1*2 + 1, rxy1*2 + 1, rz1*2 + 1 );
	
	// compute ( I - I*sigma1 )^2, store in imgConv
	final Cursor<FloatType> cursorImg = img.createCursor();
	final Cursor<FloatType> cursorConv = imgConv.createCursor();
	
	while ( cursorImg.hasNext() )
	{
		cursorImg.fwd();
		cursorConv.fwd();
		
		final float diff = cursorImg.getType().get() - cursorConv.getType().get();
		
		cursorConv.getType().set( diff*diff );
	}
	
	// compute ( ( I - I*sigma1 )^2 ) * sigma2, store in imgConv
	IntegralImage3d.computeIntegralImage( integralImg, imgConv );
	
	DOM.meanMirror( integralImg, imgConv, rxy2*2 + 1, rxy2*2 + 1, rz2*2 + 1 );
	
	integralImg.close();
	
	ViewDataBeads.normalizeImage( imgConv );
	
	return imgConv;
}
 

private static final void copy( final long start, final long loopSize, final Image< FloatType > source, final Image< FloatType > block, final int[] offset, final boolean inside, final OutOfBoundsStrategyFactory< FloatType > strategyFactory )
{
	final int numDimensions = source.getNumDimensions();
	final LocalizableCursor<FloatType> cursor = block.createLocalizableCursor();
	final LocalizableByDimCursor<FloatType> randomAccess;
	
	if ( inside )
		randomAccess = source.createLocalizableByDimCursor();
	else
		randomAccess = source.createLocalizableByDimCursor( strategyFactory );
	
	cursor.fwd( start );
	
	final int[] tmp = new int[ numDimensions ];
	
	for ( long l = 0; l < loopSize; ++l )
	{
		cursor.fwd();
		cursor.getPosition( tmp );
		
		for ( int d = 0; d < numDimensions; ++d )
			tmp[ d ] += offset[ d ];
		
		randomAccess.setPosition( tmp );
		cursor.getType().set( randomAccess.getType() );
	}
}
 

private final static <T extends Type<T>> float[] extractSliceFloat( final Image<T> img, final Display<T> display, final int dimX, final int dimY, final int[] dimensionPositions )
  {
final int sizeX = img.getDimension( dimX );
final int sizeY = img.getDimension( dimY );
  	
  	final LocalizablePlaneCursor<T> cursor = img.createLocalizablePlaneCursor();		
cursor.reset( dimX, dimY, dimensionPositions );   	

// store the slice image
  	float[] sliceImg = new float[ sizeX * sizeY ];
  	
  	if ( dimY < img.getNumDimensions() )
  	{
   	while ( cursor.hasNext() )
   	{
   		cursor.fwd();
   		sliceImg[ cursor.getPosition( dimX ) + cursor.getPosition( dimY ) * sizeX ] = display.get32Bit( cursor.getType() );    		
   	}
  	}
  	else // only a 1D image
  	{
   	while ( cursor.hasNext() )
   	{
   		cursor.fwd();
   		sliceImg[ cursor.getPosition( dimX ) ] = display.get32Bit( cursor.getType() );    		
   	}    		
  	}
  	
  	cursor.close();

  	return sliceImg;
  }
 

public static Image<FloatType> computeExponentialKernel( final Image<FloatType> kernel, final int numViews )
{
	final Image<FloatType> exponentialKernel = kernel.clone();
	
	for ( final FloatType f : exponentialKernel )
		f.set( pow( f.get(), numViews ) );
	
	//IJ.log("Jusrt using numViews/2 as exponent" );
	
	return exponentialKernel;
}
 

private static final Image<UnsignedByteType> wrap(final byte[] p, final int[] dims) {
	final Array<UnsignedByteType,ByteArray> c = new Array<UnsignedByteType,ByteArray>(
		new ArrayContainerFactory(),
		new ByteArray(p),
		dims, 1);
	final UnsignedByteType t = new UnsignedByteType(c);
	c.setLinkedType(t);
	return new Image<UnsignedByteType>(c, t);
}
 

public void setKernel( final Image<FloatType> kernel ) 
{
	this.kernel1 = kernel;
	
	init( iterationType, views );

	setCurrentIteration( -1 );
}
 

public ExtractPSF( final ViewStructure viewStructure )
{		
	this.viewStructure = viewStructure;
	this.pointSpreadFunctions = new ArrayList<Image<FloatType>>();
	this.originalPSFs = new ArrayList<Image<FloatType>>();
	this.conf = viewStructure.getSPIMConfiguration();
	
	setPSFSize( Multi_View_Deconvolution.psfSize, Multi_View_Deconvolution.isotropic, Multi_View_Deconvolution.psfSize3d );
}
 

@Override
public ArrayList<DifferenceOfGaussianPeak< FloatType>> findPeaks( final Image< FloatType > laPlace )
{
	IOFunctions.println( new Date( System.currentTimeMillis() ) + ": Detecting peaks." );
	simplePeaks = InteractiveIntegral.findPeaks( laPlace, (float)min );

	return new ArrayList<DifferenceOfGaussianPeak< FloatType>>();
}
 

public DifferenceOfGaussianNewPeakFinder(
		final Image< FloatType> img, OutOfBoundsStrategyFactory< FloatType> outOfBoundsFactory,
		final double[] sigma1, final double[] sigma2, double minPeakValue, double normalizationFactor)
{
	super( img, outOfBoundsFactory, sigma1, sigma2, minPeakValue, normalizationFactor );

	this.s1 = sigma1;
	this.s2 = sigma2;
	this.min = minPeakValue;
}
 

public ImgLibVolume(final Image<T> img, final float[] origin) throws Exception {
	super();
	if (img.getNumDimensions() < 3) throw new Exception("Image does not support at least 3 dimensions.");

	this.img = img;
	this.xDim = img.getDimension(0);
	this.yDim = img.getDimension(1);
	this.zDim = img.getDimension(2);

	this.pw = img.getCalibration(0);
	this.ph = img.getCalibration(1);
	this.pd = img.getCalibration(2);

	System.out.println("dims: " + xDim + ", " + yDim + ", " + zDim + " :: " + pw +", " + ph + ", " + pd);

	float xSpace = (float)pw;
	float ySpace = (float)ph;
	float zSpace = (float)pd;

	// real coords
	minCoord.x = origin[0];
	minCoord.y = origin[1];
	minCoord.z = origin[2];

	maxCoord.x = minCoord.x + xDim * xSpace;
	maxCoord.y = minCoord.y + yDim * ySpace;
	maxCoord.z = minCoord.z + zDim * zSpace;

	initLoader();
}
 

final public static Image<FloatType> createImageFromArray( final float[] data, final int[] dim )
{
    final FloatAccess access = new FloatArray( data );
    final Array<FloatType, FloatAccess> array = 
        new Array<FloatType, FloatAccess>(new ArrayContainerFactory(), access, dim, 1 );
        
    // create a Type that is linked to the container
    final FloatType linkedType = new FloatType( array );
    
    // pass it to the DirectAccessContainer
    array.setLinkedType( linkedType );
    
    return new Image<FloatType>(array, new FloatType());
}
 

@Override
public Image<FloatType> getResultImage() 
{
	return entropy;
}
 

final public static void computeDifferencOfMean( final Image< LongType> integralImg, final Image< FloatType > domImg, final int sx1, final int sy1, final int sz1, final int sx2, final int sy2, final int sz2, final float min, final float max  )
{
	final float diff = max - min;
	
	final float sumPixels1 = sx1 * sy1 * sz1;
	final float sumPixels2 = sx2 * sy2 * sz2;
	
	final float d1 = sumPixels1 * diff;
	final float d2 = sumPixels2 * diff;
	
	final int sx1Half = sx1 / 2;
	final int sy1Half = sy1 / 2;
	final int sz1Half = sz1 / 2;

	final int sx2Half = sx2 / 2;
	final int sy2Half = sy2 / 2;
	final int sz2Half = sz2 / 2;
	
	final int sxHalfMax = Math.max( sx1Half, sx2Half );
	final int syHalfMax = Math.max( sy1Half, sy2Half );
	final int szHalfMax = Math.max( sz1Half, sz2Half );

	final int w = domImg.getDimension( 0 ) - ( Math.max( sx1, sx2 ) / 2 ) * 2;
	final int h = domImg.getDimension( 1 ) - ( Math.max( sy1, sy2 ) / 2 ) * 2;
	final int d = domImg.getDimension( 2 ) - ( Math.max( sz1, sz2 ) / 2 ) * 2;
			
	final long imageSize = domImg.getNumPixels();

	final AtomicInteger ai = new AtomicInteger(0);					
       final Thread[] threads = SimpleMultiThreading.newThreads();

       final Vector<Chunk> threadChunks = SimpleMultiThreading.divideIntoChunks( imageSize, threads.length );
	
       for (int ithread = 0; ithread < threads.length; ++ithread)
           threads[ithread] = new Thread(new Runnable()
           {
               public void run()
               {
               	// Thread ID
               	final int myNumber = ai.getAndIncrement();
       
               	final int[] position = new int[ 3 ];
               	
               	// get chunk of pixels to process
               	final Chunk myChunk = threadChunks.get( myNumber );
               	final long loopSize = myChunk.getLoopSize();
               	
           		final LocalizableCursor< FloatType > cursor = domImg.createLocalizableCursor();
           		final LocalizableByDimCursor< LongType > randomAccess = integralImg.createLocalizableByDimCursor();

           		cursor.fwd( myChunk.getStartPosition() );
           		
           		// do as many pixels as wanted by this thread
                   for ( long j = 0; j < loopSize; ++j )
                   {                    	
           			final FloatType result = cursor.next();
           			
           			final int x = cursor.getPosition( 0 );
           			final int y = cursor.getPosition( 1 );
           			final int z = cursor.getPosition( 2 );
           			
           			final int xt = x - sxHalfMax;
           			final int yt = y - syHalfMax;
           			final int zt = z - szHalfMax;
           			
           			if ( xt >= 0 && yt >= 0 && zt >= 0 && xt < w && yt < h && zt < d )
           			{
           				position[ 0 ] = x - sx1Half;
           				position[ 1 ] = y - sy1Half;
           				position[ 2 ] = z - sz1Half;
           				randomAccess.setPosition( position );
           				final float s1 = computeSum2( sx1, sy1, sz1, randomAccess ) / d1;
           				
           				position[ 0 ] = x - sx2Half;
           				position[ 1 ] = y - sy2Half;
           				position[ 2 ] = z - sz2Half;
           				randomAccess.setPosition( position );
           				final float s2 = computeSum2( sx2, sy2, sz2, randomAccess ) / d2;

           				result.set( ( s2 - s1 ) );
             			}
           			else
           			{
           				result.set( 0 );
           			}
                   }
               }
           });
       
       SimpleMultiThreading.startAndJoin( threads );
       
       //ImageJFunctions.show( tmp1  ).setTitle( "s2" );
       //ImageJFunctions.show( tmp2  ).setTitle( "s1" );
       //ImageJFunctions.show( tmp3  ).setTitle( "1" );
	}
 

public EntropyFloatArray3D( final float stepSize, final Image<FloatType> img, final LocalizableByDimCursor3D<FloatType> cIn, final LocalizableByDimCursor3D<FloatType> cOut,
							final int histogramBins, final int windowSizeX, final int windowSizeY, final int windowSizeZ, int x, int y, int z)
{
	absFreq = new int[histogramBins];

	this.img = img;
	this.cIn = cIn;
	this.cOut = cOut;
	this.histogramBins = histogramBins;
	
	if (windowSizeX %2 == 0)
		this.windowSizeX = windowSizeX + 1;
	else
		this.windowSizeX = windowSizeX;
		
	if (windowSizeY %2 == 0)
		this.windowSizeY = windowSizeY + 1;
	else
		this.windowSizeY = windowSizeY;

	if (windowSizeZ %2 == 0)
		this.windowSizeZ = windowSizeZ + 1;
	else
		this.windowSizeZ = windowSizeZ;
	
	size = this.windowSizeX * this.windowSizeY * this.windowSizeZ;
	
	windowSizeXHalf = this.windowSizeX/2;
	windowSizeYHalf = this.windowSizeY/2;
	windowSizeZHalf = this.windowSizeZ/2;
	
	if (z - windowSizeZHalf >= 0 && z + windowSizeZHalf < img.getDimension(2) )
		outOfImageZ = false;
	else
		outOfImageZ = true;
	
	if (y - windowSizeYHalf >= 0 && y + windowSizeYHalf < img.getDimension(1) )
		outOfImageY = false;
	else 
		outOfImageY = true;
	
	if (x - windowSizeXHalf >= 0 && x + windowSizeXHalf < img.getDimension(0))
		outOfImageX = false;
	else
		outOfImageX = true;		
	
	preComputed = preComputeProbabilities(stepSize);
	
	this.x = x;
	this.y = y;		
	this.z = z;
}
 

@Override
public Image<FloatType> getResultImage() {
	return gaussContent;
}
 

public static <T extends Type<T>> boolean saveAsTiffs( final Image<T> img, String directory, final int type )
{
	return saveAsTiffs( img, directory, img.getName(), type );
}
 

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

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