Java Code Examples for mpicbg.imglib.image.Image#getDimension()

private static final void copy3dArray( final int threadIdx, final int numThreads, final Image< FloatType > source, final Image< FloatType > block, final int[] offset, 
		final boolean inside, final OutOfBoundsStrategyFactory< FloatType > strategyFactory )
{
	final int w = block.getDimension( 0 );
	final int h = block.getDimension( 1 );
	final int d = block.getDimension( 2 );
	
	final int offsetX = offset[ 0 ];
	final int offsetY = offset[ 1 ];
	final int offsetZ = offset[ 2 ];
	final float[] blockArray = ((FloatArray)((Array)block.getContainer()).update( null )).getCurrentStorageArray();
	
	final LocalizableByDimCursor3D<FloatType> randomAccess;
	
	if ( inside )
		randomAccess = (LocalizableByDimCursor3D<FloatType>)source.createLocalizableByDimCursor();
	else
		randomAccess = (LocalizableByDimCursor3D<FloatType>)source.createLocalizableByDimCursor( strategyFactory );		
	
	for ( int z = threadIdx; z < d; z += numThreads )
	{
		randomAccess.setPosition( offsetX, offsetY, z + offsetZ );
		
		int i = z * h * w;
		
		for ( int y = 0; y < h; ++y )
		{
			randomAccess.setPosition( offsetX, 0 );
			
			for ( int x = 0; x < w; ++x )
			{
				blockArray[ i++ ] = randomAccess.getType().get();
				randomAccess.fwdX();
			}
			
			randomAccess.move( -w, 0 );
			randomAccess.fwdY();
		}
	}
}
 

protected ArrayList<Integer> getNeighboringLabels( final Image<IntType> connectedComponents, final ArrayList<Point3i> neighbors, final int x, final int y, final int z )
{
	final ArrayList<Integer> labels = new ArrayList<Integer>();
	final Iterator<Point3i> iterateNeighbors = neighbors.iterator();
	
	final int w = connectedComponents.getDimension( 0 );
	final int h = connectedComponents.getDimension( 1 );
	final int d = connectedComponents.getDimension( 2 );
	
	final LocalizableByDimCursor3D<IntType> cursor = (LocalizableByDimCursor3D<IntType>) connectedComponents.createLocalizableByDimCursor();

	while (iterateNeighbors.hasNext())
	{
		Point3i neighbor = iterateNeighbors.next();

		int xp = x + neighbor.x;
		int yp = y + neighbor.y;
		int zp = z + neighbor.z;

		if (xp >= 0 && yp >= 0 && zp >= 0 && xp < w && yp < h && zp < d )
		{
			cursor.setPosition( xp, yp, zp );
			int label = cursor.getType().get();

			if (label != 0 && !labels.contains(neighbor)) 
				labels.add(label);
		}
	}
	
	cursor.close();

	return labels;
}
 

public void extract( final int viewID, final int[] maxSize )
{
	final ArrayList<ViewDataBeads > views = viewStructure.getViews();
	final int numDimensions = 3;
	
	final int[] size;
	
	if ( this.size3d == null )
	{
		size = Util.getArrayFromValue( this.size, numDimensions );
		if ( !this.isotropic )
		{
			size[ numDimensions - 1 ] *= Math.max( 1, 5.0/views.get( 0 ).getZStretching() );
			if ( size[ numDimensions - 1 ] % 2 == 0 )
				size[ numDimensions - 1 ]++;
		}
	}
	else
	{
		size = this.size3d.clone();
	}
	
	IJ.log ( "PSF size: " + Util.printCoordinates( size ) );
	
	final ViewDataBeads view = views.get( viewID );		
		
	final Image<FloatType> originalPSF = extractPSF( view, size );
	final Image<FloatType> psf = transformPSF( originalPSF, (AbstractAffineModel3D<?>)view.getTile().getModel() );
	psf.setName( "PSF_" + view.getName() );
	
	for ( int d = 0; d < numDimensions; ++d )
		if ( psf.getDimension( d ) > maxSize[ d ] )
			maxSize[ d ] = psf.getDimension( d );
	
	pointSpreadFunctions.add( psf );
	originalPSFs.add( originalPSF );
	
	psf.getDisplay().setMinMax();
}
 

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

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

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

/**
 * Averages all channels into the target image. The size is given by the dimensions of the target image,
 * the offset (if applicable) is given by an extra field
 * 
 * @param target - the target Image
 * @param offset - the offset of the area (might be [0,0] or [0,0,0])
 * @param sources - a list of input Images
 */
protected static < T extends RealType< T >, S extends RealType< S > > void averageAllChannels( final Image< T > target, final ArrayList< Image< S > > sources, final int[] offset )
{
	// get the major numbers
	final int numDimensions = target.getNumDimensions();
	final float numImages = sources.size();
	long imageSize = target.getDimension( 0 );
	
	for ( int d = 1; d < target.getNumDimensions(); ++d )
		imageSize *= target.getDimension( d );

	// run multithreaded
	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()
           {
               @Override
               public void run()
               {
               	// Thread ID
               	final int myNumber = ai.getAndIncrement();
       
               	// get chunk of pixels to process
               	final Chunk myChunk = threadChunks.get( myNumber );
               	final long startPos = myChunk.getStartPosition();
               	final long loopSize = myChunk.getLoopSize();
               	
           		// the cursor for the output
           		final LocalizableCursor< T > targetCursor =  target.createLocalizableCursor();
           		
           		// the input cursors
           		final ArrayList< LocalizableByDimCursor< S > > sourceCursors = new ArrayList< LocalizableByDimCursor< S > > ();
           		
           		for ( final Image< S > source : sources )
           			sourceCursors.add( source.createLocalizableByDimCursor() );
           		
           		// temporary array
           		final int[] location = new int[ numDimensions ]; 

           		// move to the starting position of the current thread
           		targetCursor.fwd( startPos );
                   
           		// do as many pixels as wanted by this thread
                   for ( long j = 0; j < loopSize; ++j )
           		{
           			targetCursor.fwd();
           			targetCursor.getPosition( location );
           			
           			for ( int d = 0; d < numDimensions; ++d )
           				location[ d ] += offset[ d ];
           			
           			float sum = 0;
           			
           			for ( final LocalizableByDimCursor< S > sourceCursor : sourceCursors )
           			{
           				sourceCursor.setPosition( location );
           				sum += sourceCursor.getType().getRealFloat();
           			}
           			
           			targetCursor.getType().setReal( sum / numImages );
           		}                	
               }
           });
       
       SimpleMultiThreading.startAndJoin( threads );		
}
 

public Entropy(final float stepSize, final Image<FloatType> img, final LocalizableByDimCursor<FloatType> cursor,
		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.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;

	this.cursor = cursor;
	cursor.setPosition( new int[]{ x, y, z} );
}
 

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

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

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

final public static void meanMirror( final Image< LongType> integralImg, final Image< FloatType > domImg, final int sx, final int sy, final int sz )
{
	mean( integralImg, domImg, sx, sy, sz );
			
	final int sxHalf = sx / 2;
	final int syHalf = sy / 2;
	final int szHalf = sz / 2;
	
	final int sxHalf2 = sxHalf * 2;
	final int syHalf2 = syHalf * 2;
	final int szHalf2 = szHalf * 2;
	
	final int w = domImg.getDimension( 0 );
	final int h = domImg.getDimension( 1 );
	final int d = domImg.getDimension( 2 );
	
	final int w1 = w - sxHalf - 1;
	final int h1 = h - syHalf - 1;
	final int d1 = d - szHalf - 1;
	
	final AtomicInteger ai = new AtomicInteger(0);					
       final Thread[] threads = SimpleMultiThreading.newThreads();
	final int numThreads = 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 LocalizableByDimCursor< FloatType > c1 = domImg.createLocalizableByDimCursor();
           		final LocalizableByDimCursor< FloatType > c2 = domImg.createLocalizableByDimCursor();

				final int[] p1 = new int[ 3 ];
				final int[] p2 = new int[ 3 ];
		
				// fill the remaining pixels with a mirror strategy (they are mostly blended away anyways)
				for ( int z = 0; z < d; ++z )
				{
					if ( z % numThreads == myNumber )
           			{
						boolean zSmaller = z < szHalf;
						boolean zBigger = z > d1;
			
						if ( zSmaller )
							p1[ 2 ] =  szHalf2 - z;
						else if ( zBigger )
							p1[ 2 ] = 2*d1 - z;
						else
							p1[ 2 ] = z;
			
						p2[ 2 ] = z;
						
						for ( int y = 0; y < h; ++y )
						{
							boolean ySmaller = y < syHalf;
							boolean yBigger = y > h1;
			
							if ( ySmaller )
								p1[ 1 ] =  syHalf2 - y;
							else if ( yBigger )
								p1[ 1 ] = 2*h1 - y;
							else
								p1[ 1 ] = y;
							
							p2[ 1 ] = y;
							
							for ( int x = 0; x < w; ++x )
							{
								boolean xSmaller = x < sxHalf;
								boolean xBigger = x > w1;
								
								if ( xSmaller || ySmaller || zSmaller || xBigger || yBigger || zBigger )
								{
									p2[ 0 ] = x; 
									c1.setPosition( p2 );
									
									if ( xSmaller )
										p1[ 0 ] =  sxHalf2 - x;
									else if ( xBigger )
										p1[ 0 ] = 2*w1 - x;
									else
										p1[ 0 ] = x;
									
									c2.setPosition( p1 );
									c1.getType().set( c2.getType().get() );
								}
							}
						}
           			}
				}
               }
           });

       SimpleMultiThreading.startAndJoin( threads );
	}
 

final public static void mean( final Image< LongType> integralImg, final Image< FloatType > domImg, final int sx, final int sy, final int sz )
{
	final float sumPixels = sx * sy * sz;
	
	final int sxHalf = sx / 2;
	final int syHalf = sy / 2;
	final int szHalf = sz / 2;

	final int w = domImg.getDimension( 0 ) - ( sx / 2 ) * 2; // this makes sense as sx is odd
	final int h = domImg.getDimension( 1 ) - ( sy / 2 ) * 2;
	final int d = domImg.getDimension( 2 ) - ( sz / 2 ) * 2;

	final AtomicInteger ai = new AtomicInteger(0);					
       final Thread[] threads = SimpleMultiThreading.newThreads();
	final int numThreads = 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();

           		// for each computation we need 8 randomaccesses, so 16 all together
           		final LocalizableByDimCursor< LongType > r1 = integralImg.createLocalizableByDimCursor();
           		final LocalizableByDimCursor< LongType > r2 = integralImg.createLocalizableByDimCursor();
           		final LocalizableByDimCursor< LongType > r3 = integralImg.createLocalizableByDimCursor();
           		final LocalizableByDimCursor< LongType > r4 = integralImg.createLocalizableByDimCursor();
           		final LocalizableByDimCursor< LongType > r5 = integralImg.createLocalizableByDimCursor();
           		final LocalizableByDimCursor< LongType > r6 = integralImg.createLocalizableByDimCursor();
           		final LocalizableByDimCursor< LongType > r7 = integralImg.createLocalizableByDimCursor();
           		final LocalizableByDimCursor< LongType > r8 = integralImg.createLocalizableByDimCursor();
           		
           		final LocalizableByDimCursor< FloatType > result = domImg.createLocalizableByDimCursor();
           		
           		final int[] p = new int[ 3 ];

           		for ( int z = 0; z < d; ++z )
           		{
           			if ( z % numThreads == myNumber )
           			{
            			for ( int y = 0; y < h; ++y )
            			{
            				// set the result randomaccess
            				p[ 0 ] = sxHalf; p[ 1 ] = y + syHalf; p[ 2 ] = z + szHalf;
            				result.setPosition( p );
            				
            				// set all randomaccess for the first box accordingly
            				p[ 0 ] = 0; p[ 1 ] = y; p[ 2 ] = z;
            				r1.setPosition( p ); // negative

            				p[ 0 ] += sx;
            				r2.setPosition( p ); // positive
            				
            				p[ 1 ] += sy;
            				r3.setPosition( p ); // negative
            				
            				p[ 0 ] -= sx;
            				r4.setPosition( p ); // positive

            				p[ 2 ] += sz;
            				r5.setPosition( p ); // negative

            				p[ 0 ] += sx;
            				r6.setPosition( p ); // positive

            				p[ 1 ] -= sy;
            				r7.setPosition( p ); // negative

            				p[ 0 ] -= sx;
            				r8.setPosition( p ); // positive

            				for ( int x = 0; x < w; ++x )
            				{
            					final long s = -r1.getType().get() + r2.getType().get() - r3.getType().get() + r4.getType().get() - r5.getType().get() + r6.getType().get() - r7.getType().get() + r8.getType().get();

            					result.getType().set( (float)s/sumPixels );
            					
            					r1.fwd( 0 ); r2.fwd( 0 ); r3.fwd( 0 ); r4.fwd( 0 ); r5.fwd( 0 ); r6.fwd( 0 ); r7.fwd( 0 ); r8.fwd( 0 );
            					result.fwd( 0 );
            				}
            			}
           			}
           		}            		
               }
           });

       SimpleMultiThreading.startAndJoin( threads );
	}
 

public static ExtractPSF loadAndTransformPSF( final ArrayList< String > fileName, final boolean transformPSFs, final ViewStructure viewStructure )
{
	ExtractPSF extractPSF = new ExtractPSF( viewStructure );
	
	final ArrayList<ViewDataBeads > views = viewStructure.getViews();
	final int numDimensions = 3;
			
	final int[] maxSize = new int[ numDimensions ];
	
	for ( int d = 0; d < numDimensions; ++d )
		maxSize[ d ] = 0;

	int i = 0;
	for ( final ViewDataBeads view : views )		
	{
        // extract the PSF for this one	        
   		if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_MAIN )
   			IOFunctions.println( "Loading PSF file '" + fileName.get( i ) + "' for " + view.getName() );

		final Image< FloatType > psfImage = LOCI.openLOCIFloatType( fileName.get( i ), viewStructure.getSPIMConfiguration().inputImageFactory );
		
		if ( psfImage == null )
		{
			IJ.log( "Could not find PSF file '" + fileName.get( i ) + "' - quitting." );
			return null;
		}
		
		++i;

		final Image<FloatType> psf;
		
		if ( transformPSFs )
		{
			if ( viewStructure.getDebugLevel() <= ViewStructure.DEBUG_MAIN )
				IOFunctions.println( "Transforming PSF for " + view.getName() );

			psf = transformPSF( psfImage, (AbstractAffineModel3D<?>)view.getTile().getModel() );
		}
		else
		{
			psf = psfImage.clone();
		}
		
		psf.setName( "PSF_" + view.getName() );
		
		for ( int d = 0; d < numDimensions; ++d )
			if ( psf.getDimension( d ) > maxSize[ d ] )
				maxSize[ d ] = psf.getDimension( d );
		
		extractPSF.pointSpreadFunctions.add( psf );
		extractPSF.originalPSFs.add( psfImage );
		
		psf.getDisplay().setMinMax();
	}
	
	extractPSF.computeAveragePSF( maxSize );
	
	return extractPSF;
}
 

/**
 * Transforms the extracted PSF using the affine transformation of the corresponding view
 * 
 * @param psf - the extracted psf (NOT z-scaling corrected)
 * @param model - the transformation model
 * @return the transformed psf which has odd sizes and where the center of the psf is also the center of the transformed psf
 */
protected static Image<FloatType> transformPSF( final Image<FloatType> psf, final AbstractAffineModel3D<?> model )
{
	// here we compute a slightly different transformation than the ImageTransform does
	// two things are necessary:
	// a) the center pixel stays the center pixel
	// b) the transformed psf has a odd size in all dimensions
	
	final int numDimensions = psf.getNumDimensions();
	
	final float[][] minMaxDim = ExtractPSF.getMinMaxDim( psf.getDimensions(), model );
	final float[] size = new float[ numDimensions ];		
	final int[] newSize = new int[ numDimensions ];		
	final float[] offset = new float[ numDimensions ];
	
	// the center of the psf has to be the center of the transformed psf as well
	// this is important!
	final double[] center = new double[ numDimensions ]; 
	
	for ( int d = 0; d < numDimensions; ++d )
		center[ d ] = psf.getDimension( d ) / 2;
	
	model.applyInPlace( center );

	for ( int d = 0; d < numDimensions; ++d )
	{						
		size[ d ] = minMaxDim[ d ][ 1 ] - minMaxDim[ d ][ 0 ];
		
		newSize[ d ] = (int)size[ d ] + 3;
		if ( newSize[ d ] % 2 == 0 )
			++newSize[ d ];
			
		// the offset is defined like this:
		// the transformed coordinates of the center of the psf
		// are the center of the transformed psf
		offset[ d ] = (float)center[ d ] - newSize[ d ]/2;
		
		//System.out.println( MathLib.printCoordinates( minMaxDim[ d ] ) + " size " + size[ d ] + " newSize " + newSize[ d ] );
	}
	
	final ImageTransform<FloatType> transform = new ImageTransform<FloatType>( psf, model, new LinearInterpolatorFactory<FloatType>( new OutOfBoundsStrategyValueFactory<FloatType>()));
	transform.setOffset( offset );
	transform.setNewImageSize( newSize );
	
	if ( !transform.checkInput() || !transform.process() )
	{
		System.out.println( "Error transforming psf: " + transform.getErrorMessage() );
		return null;
	}
	
	final Image<FloatType> transformedPSF = transform.getResult();
	
	ViewDataBeads.normalizeImage( transformedPSF );
	
	return transformedPSF;
}
 

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

private static final void paste3d( final int threadIdx, final int numThreads, final Image< FloatType > target, final Image< FloatType > block, 
		final int[] effectiveOffset, final int[] effectiveSize, final int[] effectiveLocalOffset )
{
	// min position in the output
	final int minX = effectiveOffset[ 0 ];
	final int minY = effectiveOffset[ 1 ];
	final int minZ = effectiveOffset[ 2 ];

	// max+1 of the output area
	final int maxX = effectiveSize[ 0 ] + minX;
	final int maxY = effectiveSize[ 1 ] + minY;
	final int maxZ = effectiveSize[ 2 ] + minZ;

	// size of the output area
	final int sX = effectiveSize[ 0 ];

	// min position in the output
	final int minXb = effectiveLocalOffset[ 0 ];
	final int minYb = effectiveLocalOffset[ 1 ];
	final int minZb = effectiveLocalOffset[ 2 ];

	// size of the target image
	final int w = target.getDimension( 0 );
	final int h = target.getDimension( 1 );

	// size of the block image
	final int wb = block.getDimension( 0 );
	final int hb = block.getDimension( 1 );

	final float[] blockArray = ((FloatArray)((Array)block.getContainer()).update( null )).getCurrentStorageArray();
	final float[] targetArray = ((FloatArray)((Array)target.getContainer()).update( null )).getCurrentStorageArray();
			
	for ( int z = minZ + threadIdx; z < maxZ; z += numThreads )
	{
		final int zBlock = z - minZ + minZb;
		
		int iTarget = z * h * w + minY * w + minX;
		int iBlock = zBlock * hb * wb + minYb * wb + minXb;
		
		for ( int y = minY; y < maxY; ++y )
		{			
			for ( int x = minX; x < maxX; ++x )
				targetArray[ iTarget++ ] = blockArray[ iBlock++ ];

			iTarget -= sX;
			iTarget += w;
			
			iBlock -= sX;
			iBlock += wb;
		}
	}
}
 

private static final void copy3d( final int threadIdx, final int numThreads, final Image< FloatType > source, final Image< FloatType > block, final int[] offset, 
		final boolean inside, final OutOfBoundsStrategyFactory< FloatType > strategyFactory )
{
	final int w = block.getDimension( 0 );
	final int h = block.getDimension( 1 );
	final int d = block.getDimension( 2 );
	
	final int offsetX = offset[ 0 ];
	final int offsetY = offset[ 1 ];
	final int offsetZ = offset[ 2 ];
	final float[] blockArray = ((FloatArray)((Array)block.getContainer()).update( null )).getCurrentStorageArray();
	
	final LocalizableByDimCursor<FloatType> randomAccess;
	
	if ( inside )
		randomAccess = source.createLocalizableByDimCursor();
	else
		randomAccess = source.createLocalizableByDimCursor( strategyFactory );		
	
	final int[] tmp = new int[]{ offsetX, offsetY, 0 };
	
	for ( int z = threadIdx; z < d; z += numThreads )
	{
		tmp[ 2 ] = z + offsetZ;
		randomAccess.setPosition( tmp );
		
		int i = z * h * w;
		
		for ( int y = 0; y < h; ++y )
		{
			randomAccess.setPosition( offsetX, 0 );
			
			for ( int x = 0; x < w; ++x )
			{
				blockArray[ i++ ] = randomAccess.getType().get();
				randomAccess.fwd( 0 );
			}
			
			randomAccess.move( -w, 0 );
			randomAccess.fwd( 1 );
		}
	}
}