mpicbg.imglib.cursor.LocalizableByDimCursor Java Examples

private static final void updateBinPositive(final int[] absFreq, final int histogramBins, final int x, final int y, final int z, final LocalizableByDimCursor<FloatType> cursor)
{
	// compute bins
	final int mx = x - cursor.getPosition( 0 );
	final int my = y - cursor.getPosition( 1 );
	final int mz = z - cursor.getPosition( 2 );
	
	cursor.move( mx, 0 );
	cursor.move( my, 1 );
	cursor.move( mz, 2 );
	
	int bin = (int) (cursor.getType().get() * histogramBins);

	cursor.move( -mx, 0 );
	cursor.move( -my, 1 );
	cursor.move( -mz, 2 );

	// for the case of value being exactly 1
	if (bin >= histogramBins) bin = histogramBins - 1;
	if (bin < 0) bin = 0;

	absFreq[bin]++;			
}
 

private static final void updateBinNegative(final int[] absFreq, final int histogramBins, final int x, final int y, final int z, final LocalizableByDimCursor<FloatType> cursor)
{
	// compute bins
	final int mx = x - cursor.getPosition( 0 );
	final int my = y - cursor.getPosition( 1 );
	final int mz = z - cursor.getPosition( 2 );
	
	cursor.move( mx, 0 );
	cursor.move( my, 1 );
	cursor.move( mz, 2 );
	
	int bin = (int) (cursor.getType().get() * histogramBins);

	cursor.move( -mx, 0 );
	cursor.move( -my, 1 );
	cursor.move( -mz, 2 );

	// for the case of value being exactly 1
	if (bin >= histogramBins) bin = histogramBins - 1;
	if (bin < 0) bin = 0;

	absFreq[bin]--;			
}
 

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 static final void paste( final long start, final long loopSize, final Image< FloatType > target, final Image< FloatType > block, 
		final int[] effectiveOffset, final int[] effectiveSize, final int[] effectiveLocalOffset )
{
	final int numDimensions = target.getNumDimensions();
	
	// iterate over effective size
	final LocalizableCursor<?> cursor = ArrayLocalizableCursor.createLinearCursor( effectiveSize );
	
	// read from block
	final LocalizableByDimCursor<FloatType> blockRandomAccess  = block.createLocalizableByDimCursor();
	
	// write to target		
	final LocalizableByDimCursor<FloatType> targetRandomAccess  = target.createLocalizableByDimCursor();
	
	cursor.fwd( start );
	
	final int[] tmp = new int[ numDimensions ];
	
	for ( long l = 0; l < loopSize; ++l )
	{
		cursor.fwd();
		cursor.getPosition( tmp );
		
		// move to the relative local offset where the real data starts
		for ( int d = 0; d < numDimensions; ++d )
			tmp[ d ] += effectiveLocalOffset[ d ];
		
		blockRandomAccess.setPosition( tmp );
		
		// move to the right position in the image
		for ( int d = 0; d < numDimensions; ++d )
			tmp[ d ] += effectiveOffset[ d ] - effectiveLocalOffset[ d ];

		targetRandomAccess.setPosition( tmp );

		// write the pixel
		targetRandomAccess.getType().set( blockRandomAccess.getType() );
	}
}
 

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

/**
 * Compute the average in the area, the RandomAccess needs to be positioned at the top left front corner:
 * 
 * fromX - start coordinate in x (exclusive in integral image coordinates, inclusive in image coordinates)
 * fromY - start coordinate in y (exclusive in integral image coordinates, inclusive in image coordinates)
 * fromZ - start coordinate in z (exclusive in integral image coordinates, inclusive in image coordinates)
 * 
 * @param vX - number of pixels in x
 * @param vY - number of pixels in y
 * @param vZ - number of pixels in z
 * @param randomAccess - randomAccess on the integral image
 * @return
 */
final public static long computeSum2( final int vX, final int vY, final int vZ, final LocalizableByDimCursor< LongType > randomAccess )
{
	long sum = -randomAccess.getType().get();
	
	randomAccess.move( vX, 0 );
	sum += randomAccess.getType().get();
	
	randomAccess.move( vY, 1 );
	sum += -randomAccess.getType().get();
	
	randomAccess.move( -vX, 0 );
	sum += randomAccess.getType().get();
	
	randomAccess.move( vZ, 2 );
	sum += -randomAccess.getType().get();
	
	randomAccess.move( vX, 0 );
	sum += randomAccess.getType().get();
	
	randomAccess.move( -vY, 1 );
	sum += -randomAccess.getType().get();
	
	randomAccess.move( -vX, 0 );
	sum += randomAccess.getType().get();
	
	return sum;
}
 

/**
 * Compute the average in the area
 * 
 * @param fromX - start coordinate in x (exclusive in integral image coordinates, inclusive in image coordinates)
 * @param fromY - start coordinate in y (exclusive in integral image coordinates, inclusive in image coordinates)
 * @param fromZ - start coordinate in z (exclusive in integral image coordinates, inclusive in image coordinates)
 * @param vX - number of pixels in x
 * @param vY - number of pixels in y
 * @param vZ - number of pixels in z
 * @param randomAccess - randomAccess on the integral image
 * @return
 */
final public static long computeSum( final int fromX, final int fromY, final int fromZ, final int vX, final int vY, final int vZ, final LocalizableByDimCursor< LongType > randomAccess )
{
	randomAccess.setPosition( fromX, 0 );
	randomAccess.setPosition( fromY, 1 );
	randomAccess.setPosition( fromZ, 2 );
	
	long sum = -randomAccess.getType().get();
	
	randomAccess.move( vX, 0 );
	sum += randomAccess.getType().get();
	
	randomAccess.move( vY, 1 );
	sum += -randomAccess.getType().get();
	
	randomAccess.move( -vX, 0 );
	sum += randomAccess.getType().get();
	
	randomAccess.move( vZ, 2 );
	sum += -randomAccess.getType().get();
	
	randomAccess.move( vX, 0 );
	sum += randomAccess.getType().get();
	
	randomAccess.move( -vY, 1 );
	sum += -randomAccess.getType().get();
	
	randomAccess.move( -vX, 0 );
	sum += randomAccess.getType().get();
	
	return sum;
}
 

private final void setPosition(final LocalizableByDimCursor<?> c, final int x, final int y) {
	position[0] = x;
	position[1] = y;
	c.setPosition(position);
}
 

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

final public static double[] computeDerivativeVector3( final LocalizableByDimCursor<FloatType> cursor )
{
	final double[] derivativeVector = new double[3];

	// x
	// derivativeVector[0] = (accessor.getRelative(1, 0, 0) - accessor.getRelative(-1, 0, 0)) / 2;
	//
	cursor.fwd( 0 );
	// we are now at (1, 0, 0)
	derivativeVector[0] = cursor.getType().get();
	cursor.bck( 0 );
	cursor.bck( 0 );
	// we are now at (-1, 0, 0)
	derivativeVector[0] -= cursor.getType().get();
	derivativeVector[0] /= 2.0f;
	
	// y
	// derivativeVector[1] = (accessor.getRelative(0, 1, 0) - accessor.getRelative(0, -1, 0)) / 2;
	//
	cursor.fwd( 0 );
	cursor.fwd( 1 );
	// we are now at (0, 1, 0)
	derivativeVector[1] = cursor.getType().get();
	cursor.bck( 1 );
	cursor.bck( 1 );
	// we are now at (0, -1, 0)
	derivativeVector[1] -= cursor.getType().get();
	derivativeVector[1] /= 2.0f;

	// z
	// derivativeVector[2] = (accessor.getRelative(0, 0, 1) - accessor.getRelative(0, 0, -1)) / 2;
	//
	cursor.fwd( 1 );
	cursor.fwd( 2 );
	// we are now at (0, 0, 1)
	derivativeVector[2] = cursor.getType().get();
	cursor.bck( 2 );
	cursor.bck( 2 );
	// we are now at (0, 0, -1)
	derivativeVector[2] -= cursor.getType().get();
	derivativeVector[2] /= 2.0f;
	
	// we are now at (0, 0, 0)
	cursor.fwd( 2 );

	return derivativeVector;

}
 

@Override
public LocalizableByDimCursor<FloatType> getResultIterator( OutOfBoundsStrategyFactory<FloatType> factory )
{
       // the iterator we need to get values from the weightening image
	return gaussContent.createLocalizableByDimCursor( factory );
}
 

@Override
public LocalizableByDimCursor<FloatType> getResultIterator()
{
       // the iterator we need to get values from the weightening image
	return gaussContent.createLocalizableByDimCursor();
}
 

public static Entropy initEntropy( final Image<FloatType> img, final int histogramBins, final int windowSizeX, final int windowSizeY, final int windowSizeZ, final int x, final int y, final int z)
{
	final LocalizableByDimCursor<FloatType> cursor = img.createLocalizableByDimCursor( new OutOfBoundsStrategyMirrorFactory<FloatType>() );
	
	// arrays for guessing the probabilities
	final Entropy ei = new Entropy(0.001f, img, cursor, histogramBins, windowSizeX, windowSizeY, windowSizeZ, x, y, z);
	
	//
	// fill arrays
	//
	for (int zs = z - ei.windowSizeZHalf; zs <= z + ei.windowSizeZHalf; zs++)
		for (int ys = y - ei.windowSizeYHalf; ys <= y + ei.windowSizeYHalf; ys++)
			for (int xs = x - ei.windowSizeXHalf; xs <= x + ei.windowSizeXHalf; xs++)
			{
				// compute bin
				cursor.move( xs - cursor.getPosition( 0 ), 0 );
				cursor.move( ys - cursor.getPosition( 1 ), 1 );
				cursor.move( zs - cursor.getPosition( 2 ), 2 );
				
				int bin = (int) (cursor.getType().get() * histogramBins);

				// for the case of value being exactly 1
				if (bin >= histogramBins) bin = histogramBins - 1;
				if (bin < 0) bin = 0;

				ei.absFreq[bin]++;					
			}
			
	//
	// make probablities and compute the entropy
	//		
	ei.entropy = 0;
	for (int bin = 0; bin < histogramBins; bin++)
	{
		if (ei.absFreq[bin] > 0)
		{
			final float prob = ei.absFreq[bin] / ei.size;				
			ei.entropy -= prob * (Math.log(prob) / Math.log(2));
		}
	}
	
	IOFunctions.println(ei.entropy);
	
	return ei;
}
 

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

@Override
public LocalizableByDimCursor<FloatType> getResultIterator( OutOfBoundsStrategyFactory<FloatType> factory )
{
       // the iterator we need to get values from the weightening image
	return gaussContent.createLocalizableByDimCursor( factory );
}
 

@Override
public LocalizableByDimCursor<FloatType> getResultIterator()
{
       // the iterator we need to get values from the weightening image
	return gaussContent.createLocalizableByDimCursor();
}
 

@Override
public LocalizableByDimCursor<FloatType> getResultIterator( OutOfBoundsStrategyFactory<FloatType> factory )
{
       // the iterator we need to get values from the weightening image
	return entropy.createLocalizableByDimCursor( factory );
}
 

@Override
public LocalizableByDimCursor<FloatType> getResultIterator()
{
       // the iterator we need to get values from the weightening image
	return entropy.createLocalizableByDimCursor();
}
 

public static ArrayList<SimplePeak> findPeaks( final Image<FloatType> laPlace, final float minValue )
	{
	    final AtomicInteger ai = new AtomicInteger( 0 );
	    final Thread[] threads = SimpleMultiThreading.newThreads( Threads.numThreads() );
	    final int nThreads = threads.length;
	    final int numDimensions = laPlace.getNumDimensions();
	    
	    final Vector< ArrayList<SimplePeak> > threadPeaksList = new Vector< ArrayList<SimplePeak> >();
	    
	    for ( int i = 0; i < nThreads; ++i )
	    	threadPeaksList.add( new ArrayList<SimplePeak>() );

		for (int ithread = 0; ithread < threads.length; ++ithread)
	        threads[ithread] = new Thread(new Runnable()
	        {
	            public void run()
	            {
	            	final int myNumber = ai.getAndIncrement();
            	
	            	final ArrayList<SimplePeak> myPeaks = threadPeaksList.get( myNumber );	
	            	final LocalizableByDimCursor<FloatType> cursor = laPlace.createLocalizableByDimCursor();	            	
	            	final LocalNeighborhoodCursor<FloatType> neighborhoodCursor = LocalNeighborhoodCursorFactory.createLocalNeighborhoodCursor( cursor );
	            	
	            	final int[] position = new int[ numDimensions ];
	            	final int[] dimensionsMinus2 = laPlace.getDimensions();

            		for ( int d = 0; d < numDimensions; ++d )
            			dimensionsMinus2[ d ] -= 2;
	            	
MainLoop:           while ( cursor.hasNext() )
	                {
	                	cursor.fwd();
	                	cursor.getPosition( position );
	                	
	                	if ( position[ 0 ] % nThreads == myNumber )
	                	{
	                		for ( int d = 0; d < numDimensions; ++d )
	                		{
	                			final int pos = position[ d ];
	                			
	                			if ( pos < 1 || pos > dimensionsMinus2[ d ] )
	                				continue MainLoop;
	                		}

	                		// if we do not clone it here, it might be moved along with the cursor
	                		// depending on the container type used
	                		final float currentValue = cursor.getType().get();
	                		
	                		// it can never be a desired peak as it is too low
	                		if ( Math.abs( currentValue ) < minValue )
                				continue;

                			// update to the current position
                			neighborhoodCursor.update();

                			// we have to compare for example 26 neighbors in the 3d case (3^3 - 1) relative to the current position
                			final SpecialPoint specialPoint = isSpecialPoint( neighborhoodCursor, currentValue ); 
                			
                			if ( specialPoint == SpecialPoint.MIN )
                				myPeaks.add( new SimplePeak( position, Math.abs( currentValue ), true, false ) ); //( position, currentValue, specialPoint ) );
                			else if ( specialPoint == SpecialPoint.MAX )
                				myPeaks.add( new SimplePeak( position, Math.abs( currentValue ), false, true ) ); //( position, currentValue, specialPoint ) );
                			
                			// reset the position of the parent cursor
                			neighborhoodCursor.reset();	                				                		
	                	}
	                }
                
	                cursor.close();
            }
        });
	
		SimpleMultiThreading.startAndJoin( threads );		

		// put together the list from the various threads	
		final ArrayList<SimplePeak> dogPeaks = new ArrayList<SimplePeak>();
		
		for ( final ArrayList<SimplePeak> peakList : threadPeaksList )
			dogPeaks.addAll( peakList );		
		
		return dogPeaks;
	}
 

final public static void computeDifferencOfMeanSlice( final Image< LongType> integralImg, final Image< FloatType > sliceImg, final int z, 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 sumPixels1 = sx1 * sy1 * sz1;
	final float sumPixels2 = sx2 * sy2 * sz2;
	
	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 = sliceImg.getDimension( 0 ) - ( Math.max( sx1, sx2 ) / 2 ) * 2;
	final int h = sliceImg.getDimension( 1 ) - ( Math.max( sy1, sy2 ) / 2 ) * 2;
	final int d = (integralImg.getDimension( 2 ) - 1) - ( Math.max( sz1, sz2 ) / 2 ) * 2;
			
	final long imageSize = w * h;

	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();
       
               	// get chunk of pixels to process
               	final Chunk myChunk = threadChunks.get( myNumber );
               	final long loopSize = myChunk.getLoopSize();
               	
           		final LocalizableCursor< FloatType > cursor = sliceImg.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 )
           			{
           				final float s1 = DOM.computeSum( x - sx1Half, y - sy1Half, z - sz1Half, sx1, sy1, sz1, randomAccess ) / sumPixels1;
           				final float s2 = DOM.computeSum( x - sx2Half, y - sy2Half, z - sz2Half, sx2, sy2, sz2, randomAccess ) / sumPixels2;
           			
           				result.set( ( s2 - s1 ) / ( max - min ) );
           			}
                   }
               }
           });
       
       SimpleMultiThreading.startAndJoin( threads );
	}
 

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

/**
 * Get projection along the smallest dimension (which is usually the rotation axis)
 * 
 * @return - the averaged, projected PSF
 */
public Image< FloatType > getMaxProjectionAveragePSF()
{
	final int[] dimensions = avgPSF.getDimensions();
	
	int minSize = dimensions[ 0 ];
	int minDim = 0;
	
	for ( int d = 0; d < dimensions.length; ++d )
	{
		if ( avgPSF.getDimension( d ) < minSize )
		{
			minSize = avgPSF.getDimension( d );
			minDim = d;
		}
	}
	
	final int[] projDim = new int[ dimensions.length - 1 ];
	
	int dim = 0;
	int sizeProjection = 0;
	
	// the new dimensions
	for ( int d = 0; d < dimensions.length; ++d )
		if ( d != minDim )
			projDim[ dim++ ] = dimensions[ d ];
		else
			sizeProjection = dimensions[ d ];
	
	final Image< FloatType > proj = avgPSF.getImageFactory().createImage( projDim );
	
	final LocalizableByDimCursor< FloatType > psfIterator = avgPSF.createLocalizableByDimCursor();
	final LocalizableCursor< FloatType > projIterator = proj.createLocalizableCursor();
	
	final int[] tmp = new int[ avgPSF.getNumDimensions() ];
	
	while ( projIterator.hasNext() )
	{
		projIterator.fwd();

		dim = 0;
		for ( int d = 0; d < dimensions.length; ++d )
			if ( d != minDim )
				tmp[ d ] = projIterator.getPosition( dim++ );

		tmp[ minDim ] = -1;
		
		float maxValue = -Float.MAX_VALUE;
		
		psfIterator.setPosition( tmp );
		for ( int i = 0; i < sizeProjection; ++i )
		{
			psfIterator.fwd( minDim );
			final float value = psfIterator.getType().get();
			
			if ( value > maxValue )
				maxValue = value;
		}
		
		projIterator.getType().set( maxValue );
	}
	
	proj.setName( "MIP of PSF's of " + viewStructure.getID() );
	
	return proj;
}
 

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

public abstract LocalizableByDimCursor<FloatType> getResultIterator(); 

public abstract LocalizableByDimCursor<FloatType> getResultIterator( OutOfBoundsStrategyFactory<FloatType> factory);