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

public void copyBlock( final Image< FloatType > source, final Image< FloatType > block )
{	
	final AtomicInteger ai = new AtomicInteger(0);					
       final Thread[] threads = SimpleMultiThreading.newThreads( numThreads );
       for ( int ithread = 0; ithread < threads.length; ++ithread )
           threads[ithread] = new Thread(new Runnable()
           {
               public void run()
               {
               	final int threadIdx = ai.getAndIncrement();
               	
               	// get chunk of pixels to process
               	final Chunk myChunk = threadChunks.get( threadIdx );
               	
               	if ( source.getNumDimensions() == 3 && Array.class.isInstance( source.getContainer() ) && Array.class.isInstance( block.getContainer() ) )
               		copy3dArray( threadIdx, numThreads, source, block, offset, inside, factory );
               	else if ( source.getNumDimensions() == 3 && Array.class.isInstance( block.getContainer() ) )
               		copy3d( threadIdx, numThreads, source, block, offset, inside, factory );
               	else
               		copy( myChunk.getStartPosition(), myChunk.getLoopSize(), source, block, offset, inside, factory );
               }
           });
       
       SimpleMultiThreading.startAndJoin( threads );

}
 

public void pasteBlock( final Image< FloatType > target, final Image< FloatType > block )
{	
	final AtomicInteger ai = new AtomicInteger(0);					
       final Thread[] threads = SimpleMultiThreading.newThreads( numThreads );
       for ( int ithread = 0; ithread < threads.length; ++ithread )
           threads[ithread] = new Thread(new Runnable()
           {
               public void run()
               {
               	final int threadIdx = ai.getAndIncrement();
               	
               	// get chunk of pixels to process
               	final Chunk myChunk = threadChunks.get( threadIdx );
               	
               	if ( target.getNumDimensions() == 3 && Array.class.isInstance( target.getContainer() ) )
               		paste3d( threadIdx, numThreads, target, block, effectiveOffset, effectiveSize, effectiveLocalOffset );
               	else
               		paste( myChunk.getStartPosition(), myChunk.getLoopSize(), target, block, effectiveOffset, effectiveSize, effectiveLocalOffset );
               }
           });
       
       SimpleMultiThreading.startAndJoin( threads );

}
 

public static void translate( final Image< FloatType > img, final float[] vector )
{
	final Image< FloatType > tmp = img.clone();
	
	final LocalizableCursor< FloatType > cursor1 = img.createLocalizableCursor();		
	final Interpolator< FloatType > interpolator = tmp.createInterpolator( new LinearInterpolatorFactory<FloatType>( new OutOfBoundsStrategyMirrorFactory<FloatType>() ) );
	
	final int numDimensions = img.getNumDimensions();
	final float[] pos = new float[ numDimensions ];
	
	while ( cursor1.hasNext() )
	{
		cursor1.fwd();
		
		for( int d = 0; d < numDimensions; ++d )
			pos[ d ] = cursor1.getPosition( d ) - vector[ d ];
		
		interpolator.setPosition( pos );
		cursor1.getType().set( interpolator.getType() );
	}
	
	cursor1.close();
	interpolator.close();
}
 

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

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

/**
 * 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 > computeInvertedKernel( final Image< FloatType > kernel )
{
	final Image< FloatType > invKernel = kernel.clone();
	
	for ( int d = 0; d < invKernel.getNumDimensions(); ++d )
		new MirrorImage< FloatType >( invKernel, d ).process();
	
	return invKernel;
}
 

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

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

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

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

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

public static Image<FloatType> computeContentBasedGauss( final Image< FloatType > img, final int fusionSigma1, final int fusionSigma2, final float zStretching )
{
	// get the kernels			
	final double[] k1 = new double[ img.getNumDimensions() ];
	final double[] k2 = new double[ img.getNumDimensions() ];
	
	for ( int d = 0; d < img.getNumDimensions() - 1; ++d )
	{
		k1[ d ] = fusionSigma1;
		k2[ d ] = fusionSigma2;
	}
	
	k1[ img.getNumDimensions() - 1 ] = fusionSigma1 / zStretching;
	k2[ img.getNumDimensions() - 1 ] = fusionSigma2 / zStretching;		
	
	final Image<FloatType> kernel1 = FourierConvolution.createGaussianKernel( new ArrayContainerFactory(), k1 );
	final Image<FloatType> kernel2 = FourierConvolution.createGaussianKernel( new ArrayContainerFactory(), k2 );

	System.out.println( new Date( System.currentTimeMillis() )  + " conv1" );
	
	// compute I*sigma1
	FourierConvolution<FloatType, FloatType> fftConv1 = new FourierConvolution<FloatType, FloatType>( img, kernel1 );
	
	fftConv1.process();		
	final Image<FloatType> conv1 = fftConv1.getResult();
	
	fftConv1.close();
	fftConv1 = null;

	System.out.println( new Date( System.currentTimeMillis() )  + " comp" );

	// compute ( I - I*sigma1 )^2
	final Cursor<FloatType> cursorImg = img.createCursor();
	final Cursor<FloatType> cursorConv = conv1.createCursor();
	
	while ( cursorImg.hasNext() )
	{
		cursorImg.fwd();
		cursorConv.fwd();
		
		final float diff = cursorImg.getType().get() - cursorConv.getType().get();
		
		cursorConv.getType().set( diff*diff );
	}

	System.out.println( new Date( System.currentTimeMillis() )  + " conv2" );

	// compute ( ( I - I*sigma1 )^2 ) * sigma2
	FourierConvolution<FloatType, FloatType> fftConv2 = new FourierConvolution<FloatType, FloatType>( conv1, kernel2 );
	fftConv2.process();	
	
	Image<FloatType> gaussContent = fftConv2.getResult();

	fftConv2.close();
	fftConv2 = null;
	
	// close the unnecessary image
	kernel1.close();
	kernel2.close();
	conv1.close();
	
	ViewDataBeads.normalizeImage( gaussContent );
	
	return gaussContent;
}
 

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

public static void gaussFit( final Image< FloatType > image, final ArrayList< Bead > beadList, final double[] typicalSigma )
{
	final AtomicInteger ai = new AtomicInteger( 0 );					
       final Thread[] threads = SimpleMultiThreading.newThreads();
       
	final Vector<Chunk> threadChunks = SimpleMultiThreading.divideIntoChunks( beadList.size(), threads.length );
	
	final int[] count = new int[ threads.length ];
	final double[][] sigmas = new double[ threads.length ][ 3 ];
	
       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 beads to process
               	final Chunk myChunk = threadChunks.get( myNumber );
               	final int loopSize = (int)myChunk.getLoopSize();
               	final int start = (int)myChunk.getStartPosition();
               	final int end = start + loopSize;
               	
               	final GaussianPeakFitterND<FloatType> fitter = new GaussianPeakFitterND<FloatType>( image );
               	
               	final float[] tmp = new float[ image.getNumDimensions() ];

           		// do as many pixels as wanted by this thread
                   for ( int j = start; j < end; ++j )
                   {
                   	final Bead bead = beadList.get( j );
                   	
                   	for ( int d = 0; d < tmp.length; ++d )
                   		tmp[ d ] = (float)bead.getL()[ d ];
                   	
           			final double[] results = fitter.process( new LocalizablePoint( tmp ), typicalSigma );
           			
           			//double a = results[ 0 ];
           			//double x = results[ 1 ];
           			//double y = results[ 2 ];
           			//double z = results[ 3 ];
           			double sx = 1/Math.sqrt( results[ 4 ] );
           			double sy = 1/Math.sqrt( results[ 5 ] );
           			double sz = 1/Math.sqrt( results[ 6 ] );
           			
           			bead.getL()[ 0 ] = bead.getW()[ 0 ] = results[ 1 ];
           			bead.getL()[ 1 ] = bead.getW()[ 1 ] = results[ 2 ];
           			bead.getL()[ 2 ] = bead.getW()[ 2 ] = results[ 3 ];
           			
           			bead.relocalized = true;
           			
           			if ( ! (Double.isNaN( sx ) || Double.isNaN( sy ) || Double.isNaN( sz ) ) )
           			{
           				sigmas[ myNumber ][ 0 ] += sx;
           				sigmas[ myNumber ][ 1 ] += sy;
           				sigmas[ myNumber ][ 2 ] += sz;
           				count[ myNumber ]++;
           			}

                   }
               }
           });
       
       SimpleMultiThreading.startAndJoin( threads );

       // compute final average sigma's
	for ( int i = 1; i < sigmas.length; ++i )
	{
		sigmas[ 0 ][ 0 ] += sigmas[ i ][ 0 ];
		sigmas[ 0 ][ 1 ] += sigmas[ i ][ 1 ];
		sigmas[ 0 ][ 2 ] += sigmas[ i ][ 2 ];
		count[ 0 ] += count[ i ];
	}
	
	IJ.log( "avg sigma: [" + ( sigmas[ 0 ][ 0 ] / count[ 0 ] ) + " px, " + ( sigmas[ 0 ][ 1 ] / count[ 0 ] ) + " px, " + ( sigmas[ 0 ][ 2 ] / count[ 0 ] ) + " px]" );
}
 

public static ArrayList< InterestPoint > computeQuadraticLocalization( final ArrayList< SimplePeak > peaks, final Image< FloatType > domImg, final boolean findMin, final boolean findMax, final float threshold, final boolean keepIntensity )
{
	IOFunctions.println("(" + new Date(System.currentTimeMillis()) + "): Subpixel localization using quadratic n-dimensional fit");

	final ArrayList< DifferenceOfGaussianPeak<FloatType> > peakList = new ArrayList<DifferenceOfGaussianPeak<FloatType>>();

	for ( final SimplePeak peak : peaks )
		if ( ( peak.isMax && findMax ) || ( peak.isMin && findMin ) )
			peakList.add( new DifferenceOfGaussianPeak<FloatType>( peak.location, new FloatType( peak.intensity ), SpecialPoint.MAX ) );
	

	final SubpixelLocalization<FloatType> spl = new SubpixelLocalization<FloatType>( domImg, peakList );
	spl.setAllowMaximaTolerance( true );
	spl.setMaxNumMoves( 10 );
	spl.setNumThreads( Threads.numThreads() );

	if ( !spl.checkInput() || !spl.process() )
		IOFunctions.println("(" + new Date(System.currentTimeMillis()) + "): Warning! Failed to compute subpixel localization " + spl.getErrorMessage() );
	
	final int n = domImg.getNumDimensions();

	final ArrayList< InterestPoint > peaks2 = new ArrayList< InterestPoint >();
	
	int id = 0;
	
	for ( DifferenceOfGaussianPeak<FloatType> detection : peakList )
	{
		if ( Math.abs( detection.getValue().get() ) > threshold )
		{
			final double[] tmp = new double[ n ];
			for ( int d = 0; d < n; ++d )
				tmp[ d ] = detection.getSubPixelPosition( d );

			if ( keepIntensity )
				peaks2.add( new InterestPointValue( id++, tmp, detection.getValue().get() ) );
			else
				peaks2.add( new InterestPoint( id++, tmp ) );
		}
	}

	return peaks2;
}
 

/**
 * 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 static < T extends RealType<T>, S extends RealType<S> > PairWiseStitchingResult computePhaseCorrelation( final Image<T> img1, final Image<S> img2, final int numPeaks, final boolean subpixelAccuracy )
{
	final PhaseCorrelation< T, S > phaseCorr = new PhaseCorrelation<T, S>( img1, img2 );
	phaseCorr.setInvestigateNumPeaks( numPeaks );
	
	if ( subpixelAccuracy )
		phaseCorr.setKeepPhaseCorrelationMatrix( true );
	
	phaseCorr.setComputeFFTinParalell( true );
	if ( !phaseCorr.process() )
	{
		Log.error( "Could not compute phase correlation: " + phaseCorr.getErrorMessage() );
		return null;
	}

	// result
	final PhaseCorrelationPeak pcp = phaseCorr.getShift();
	final float[] shift = new float[ img1.getNumDimensions() ];
	final PairWiseStitchingResult result;
	
	if ( subpixelAccuracy )
	{
		final Image<FloatType> pcm = phaseCorr.getPhaseCorrelationMatrix();		
	
		final ArrayList<DifferenceOfGaussianPeak<FloatType>> list = new ArrayList<DifferenceOfGaussianPeak<FloatType>>();		
		final Peak p = new Peak( pcp );
		list.add( p );
				
		final SubpixelLocalization<FloatType> spl = new SubpixelLocalization<FloatType>( pcm, list );
		final boolean move[] = new boolean[ pcm.getNumDimensions() ];
		for ( int i = 0; i < pcm.getNumDimensions(); ++i )
			move[ i ] = false;
		spl.setCanMoveOutside( true );
		spl.setAllowedToMoveInDim( move );
		spl.setMaxNumMoves( 0 );
		spl.setAllowMaximaTolerance( false );
		spl.process();
		
		final Peak peak = (Peak)list.get( 0 );
		
		for ( int d = 0; d < img1.getNumDimensions(); ++d )
			shift[ d ] = peak.getPCPeak().getPosition()[ d ] + peak.getSubPixelPositionOffset( d );
		
		pcm.close();
		
		result = new PairWiseStitchingResult( shift, pcp.getCrossCorrelationPeak(), p.getValue().get() );
	}
	else
	{
		for ( int d = 0; d < img1.getNumDimensions(); ++d )
			shift[ d ] = pcp.getPosition()[ d ];
		
		result = new PairWiseStitchingResult( shift, pcp.getCrossCorrelationPeak(), pcp.getPhaseCorrelationPeak() );
	}
	
	return result;
}