net.imglib2.img.basictypeaccess.array.FloatArray Java Examples

@SuppressWarnings("unchecked")
final protected static void convolve2BlockCUDA(
		final Block blockStruct, final int deviceId, final Img< FloatType > image,
		final Img< FloatType > result, final Img< FloatType > block, final Img< FloatType > kernel2 )
{
	// ratio outside of the deconvolved space (psi) is 1
	blockStruct.copyBlock( Views.extendValue( image, new FloatType( 1.0f ) ), block );

	// convolve block with kernel2 using CUDA
	final float[] blockF = ((FloatArray)((ArrayImg< net.imglib2.type.numeric.real.FloatType, ? > )block).update( null ) ).getCurrentStorageArray();
	final float[] kernel2F = ((FloatArray)((ArrayImg< net.imglib2.type.numeric.real.FloatType, ? > )kernel2).update( null ) ).getCurrentStorageArray();

	MVDeconFFT.cuda.convolution3DfftCUDAInPlace(
			blockF, getCUDACoordinates( CUDAOutput.getImgSizeInt( block ) ),
			kernel2F, getCUDACoordinates( CUDAOutput.getImgSizeInt( kernel2 ) ),
			deviceId );

	blockStruct.pasteBlock( result, block );
}
 

@SuppressWarnings("unchecked")
final protected static void convolve1BlockCUDA(
		final Block blockStruct, final int deviceId, final Img< FloatType > image,
		final Img< FloatType > result, final Img< FloatType > block, final Img< FloatType > kernel1, final int i )
{
	long time = System.currentTimeMillis();
	blockStruct.copyBlock( Views.extendMirrorSingle( image ), block );
	System.out.println( " block " + i + "(CPU  " + deviceId + "): copy " + (System.currentTimeMillis() - time) );

	// convolve block with kernel1 using CUDA
	time = System.currentTimeMillis();
	final float[] blockF = ((FloatArray)((ArrayImg< net.imglib2.type.numeric.real.FloatType, ? > )block).update( null ) ).getCurrentStorageArray();
	final float[] kernel1F = ((FloatArray)((ArrayImg< net.imglib2.type.numeric.real.FloatType, ? > )kernel1).update( null ) ).getCurrentStorageArray();
	
	MVDeconFFT.cuda.convolution3DfftCUDAInPlace(
			blockF, getCUDACoordinates( CUDAOutput.getImgSizeInt( block ) ),
			kernel1F, getCUDACoordinates( CUDAOutput.getImgSizeInt( kernel1 ) ),
			deviceId );
	System.out.println( " block " + i + "(CUDA " + deviceId + "): compute " + (System.currentTimeMillis() - time) );

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

public static void main(String[] args)
{
	final ImagePlus imp = IJ.openImage( "/Users/david/Desktop/stable HelaK-GFP-H2A.Z20000.tif" );
	new ImageJ();

	RandomAccessibleInterval< ? extends RealType > img = ImageJFunctions.wrapReal( imp );
	ArrayImg< FloatType, FloatArray > f = ArrayImgs.floats( 1024, 1024 );
	ArrayImg< FloatType, FloatArray > w = ArrayImgs.floats( 1024, 1024 );
	RandomAccessibleInterval< FloatType > interp = (RandomAccessibleInterval< FloatType >) getLinearInterpolation( img, new FloatType(), new float[] {0.5f,0.5f}, Executors.newSingleThreadExecutor() ).getA();
	
	RandomAccessibleInterval< FloatType > weight = new ArrayImgFactory( new FloatType() ).create( interp );
	applyWeights( interp, weight, new float[] {0.5f,0.5f}, new float[] {0,0}, new float[] {20,20}, false, Executors.newSingleThreadExecutor() );
	addTranslated( Views.iterable( interp ), f, new int[] {500, 700}, Executors.newSingleThreadExecutor() );
	addTranslated( Views.iterable( interp ), f, new int[] {400, 500}, Executors.newSingleThreadExecutor() );
	addTranslated( Views.iterable( weight ), w, new int[] {500, 700}, Executors.newSingleThreadExecutor() );
	addTranslated( Views.iterable( weight ), w, new int[] {400, 500}, Executors.newSingleThreadExecutor() );

	normalizeWeights( f, w, Executors.newSingleThreadExecutor() );
	
	ImageJFunctions.show( f );
}
 

public ArrayImg<FloatType, FloatArray> generateFloatArrayTestImg(
	final boolean fill, final long... dims)
{
	final float[] array = new float[(int) Intervals.numElements(
		new FinalInterval(dims))];

	if (fill) {
		seed = 17;
		for (int i = 0; i < array.length; i++) {
			array[i] = (float) pseudoRandom() / (float) Integer.MAX_VALUE;
		}
	}

	return ArrayImgs.floats(array, dims);
}
 

public static void main( final String[] args )
{
	new ImageJ();

	final double[] coefficients = new double[]{
			0, 2, 4, 8,
			1, 1, 1, 1,
			1, 10, 5, 1,
			1, 1, 1, 1,

			0, 10, 20, 30,
			40, 50, 60, 70,
			80, 90, 100, 110,
			120, 130, 140, 150
	};

	final LinearIntensityMap< DoubleType > transform = new LinearIntensityMap< DoubleType >( ArrayImgs.doubles( coefficients, 4, 4, 2 ) );

	//final ImagePlus imp = new ImagePlus( "http://upload.wikimedia.org/wikipedia/en/2/24/Lenna.png" );
	final ImagePlus imp1 = new ImagePlus( "http://fly.mpi-cbg.de/~saalfeld/Pictures/norway.jpg");

	final ArrayImg< FloatType, FloatArray > image1 = ArrayImgs.floats( ( float[] )imp1.getProcessor().convertToFloatProcessor().getPixels(), imp1.getWidth(), imp1.getHeight() );
	final ArrayImg< UnsignedByteType, ByteArray > image2 = ArrayImgs.unsignedBytes( ( byte[] )imp1.getProcessor().convertToByteProcessor().getPixels(), imp1.getWidth(), imp1.getHeight() );
	final ArrayImg< UnsignedShortType, ShortArray > image3 = ArrayImgs.unsignedShorts( ( short[] )imp1.getProcessor().convertToShortProcessor().getPixels(), imp1.getWidth(), imp1.getHeight() );
	final ArrayImg< ARGBType, IntArray > image4 = ArrayImgs.argbs( ( int[] )imp1.getProcessor().getPixels(), imp1.getWidth(), imp1.getHeight() );

	ImageJFunctions.show( ArrayImgs.doubles( coefficients, 4, 4, 2 ) );

	transform.run( image1 );
	transform.run( image2 );
	transform.run( image3 );
	transform.run( image4 );

	ImageJFunctions.show( image1 );
	ImageJFunctions.show( image2 );
	ImageJFunctions.show( image3 );
	ImageJFunctions.show( image4 );
}
 

@SuppressWarnings({ "unchecked", "rawtypes" })
private <A extends ArrayDataAccess<A>> SCIFIOCellLoader<T, A>
	createCellLoader(final NativeTypeFactory<T, A> typeFactory)
{
	switch (typeFactory.getPrimitiveType()) {
		case BYTE:
			return new SCIFIOCellLoader(new ByteArrayLoader(reader, subregion),
				o -> new ByteArray((byte[]) o));
		case CHAR:
			return new SCIFIOCellLoader(new CharArrayLoader(reader, subregion),
				o -> new CharArray((char[]) o));
		case DOUBLE:
			return new SCIFIOCellLoader(new DoubleArrayLoader(reader, subregion),
				o -> new DoubleArray((double[]) o));
		case FLOAT:
			return new SCIFIOCellLoader(new FloatArrayLoader(reader, subregion),
				o -> new FloatArray((float[]) o));
		case INT:
			return new SCIFIOCellLoader(new IntArrayLoader(reader, subregion),
				o -> new IntArray((int[]) o));
		case LONG:
			return new SCIFIOCellLoader(new LongArrayLoader(reader, subregion),
				o -> new LongArray((long[]) o));
		case SHORT:
			return new SCIFIOCellLoader(new ShortArrayLoader(reader, subregion),
				o -> new ShortArray((short[]) o));
		default:
			throw new IllegalArgumentException();
	}
}
 

private static final void copy3dArray( final int threadIdx, final int numThreads, final RandomAccessible< FloatType > source, final ArrayImg< FloatType, ? > block, final long[] offset )
{
	final int w = (int)block.dimension( 0 );
	final int h = (int)block.dimension( 1 );
	final int d = (int)block.dimension( 2 );

	final long offsetX = offset[ 0 ];
	final long offsetY = offset[ 1 ];
	final long offsetZ = offset[ 2 ];
	final float[] blockArray = ((FloatArray)block.update( null ) ).getCurrentStorageArray();

	// define where we will query the RandomAccess on the source
	final FinalInterval interval = new FinalInterval( new long[] { offsetX, offsetY, offsetZ }, new long[] { offsetX + w - 1, offsetY + h - 1, offsetZ + d - 1 } );
	final RandomAccess< FloatType > randomAccess = source.randomAccess( interval );

	final long[] tmp = new long[]{ 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.get().get();
				randomAccess.fwd( 0 );
			}

			randomAccess.move( -w, 0 );
			randomAccess.fwd( 1 );
		}
	}
}
 

@OpMethod(op = net.imagej.ops.math.ConstantToPlanarImage.SubtractFloat.class)
public PlanarImg<FloatType, FloatArray> subtract(
	final PlanarImg<FloatType, FloatArray> image, final float value)
{
	@SuppressWarnings("unchecked")
	final PlanarImg<FloatType, FloatArray> result =
		(PlanarImg<FloatType, FloatArray>) ops().run(
			net.imagej.ops.Ops.Math.Subtract.class, image, value);
	return result;
}
 

@OpMethod(ops = {
	net.imagej.ops.math.ConstantToArrayImageP.SubtractFloat.class,
	net.imagej.ops.math.ConstantToArrayImage.SubtractFloat.class })
public ArrayImg<FloatType, FloatArray> subtract(
	final ArrayImg<FloatType, FloatArray> image, final float value)
{
	@SuppressWarnings("unchecked")
	final ArrayImg<FloatType, FloatArray> result =
		(ArrayImg<FloatType, FloatArray>) ops().run(Ops.Math.Subtract.NAME, image,
			value);
	return result;
}
 

@OpMethod(op = net.imagej.ops.math.ConstantToPlanarImage.MultiplyFloat.class)
public PlanarImg<FloatType, FloatArray> multiply(
	final PlanarImg<FloatType, FloatArray> image, final float value)
{
	@SuppressWarnings("unchecked")
	final PlanarImg<FloatType, FloatArray> result =
		(PlanarImg<FloatType, FloatArray>) ops().run(
			net.imagej.ops.Ops.Math.Multiply.class, image, value);
	return result;
}
 

@OpMethod(ops = { net.imagej.ops.math.ConstantToArrayImageP.AddFloat.class,
	net.imagej.ops.math.ConstantToArrayImage.AddFloat.class })
public ArrayImg<FloatType, FloatArray> add(
	final ArrayImg<FloatType, FloatArray> image, final float value)
{
	@SuppressWarnings("unchecked")
	final ArrayImg<FloatType, FloatArray> result =
		(ArrayImg<FloatType, FloatArray>) ops().run(Ops.Math.Add.NAME, image,
			value);
	return result;
}
 

@OpMethod(op = net.imagej.ops.math.ConstantToPlanarImage.DivideFloat.class)
public PlanarImg<FloatType, FloatArray> divide(
	final PlanarImg<FloatType, FloatArray> image, final float value)
{
	@SuppressWarnings("unchecked")
	final PlanarImg<FloatType, FloatArray> result =
		(PlanarImg<FloatType, FloatArray>) ops().run(
			net.imagej.ops.Ops.Math.Divide.class, image, value);
	return result;
}
 

@OpMethod(ops = { net.imagej.ops.math.ConstantToArrayImageP.DivideFloat.class,
	net.imagej.ops.math.ConstantToArrayImage.DivideFloat.class })
public ArrayImg<FloatType, FloatArray> divide(
	final ArrayImg<FloatType, FloatArray> image, final float value)
{
	@SuppressWarnings("unchecked")
	final ArrayImg<FloatType, FloatArray> result =
		(ArrayImg<FloatType, FloatArray>) ops().run(Ops.Math.Divide.NAME, image,
			value);
	return result;
}
 

@OpMethod(op = net.imagej.ops.math.ConstantToPlanarImage.AddFloat.class)
public PlanarImg<FloatType, FloatArray> add(
	final PlanarImg<FloatType, FloatArray> image, final float value)
{
	@SuppressWarnings("unchecked")
	final PlanarImg<FloatType, FloatArray> result =
		(PlanarImg<FloatType, FloatArray>) ops().run(
			net.imagej.ops.Ops.Math.Add.class, image, value);
	return result;
}
 

@OpMethod(ops = {
	net.imagej.ops.math.ConstantToArrayImageP.MultiplyFloat.class,
	net.imagej.ops.math.ConstantToArrayImage.MultiplyFloat.class })
public ArrayImg<FloatType, FloatArray> multiply(
	final ArrayImg<FloatType, FloatArray> image, final float value)
{
	@SuppressWarnings("unchecked")
	final ArrayImg<FloatType, FloatArray> result =
		(ArrayImg<FloatType, FloatArray>) ops().run(Ops.Math.Multiply.NAME, image,
			value);
	return result;
}
 

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

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

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

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

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

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

	final float[] blockArray = ((FloatArray)block.update( null ) ).getCurrentStorageArray();
	final float[] targetArray = ((FloatArray)target.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 )
		{
			copyX( blockArray, targetArray, sX, iTarget, iBlock );

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

@Override
public boolean process()
{
	// do not operate at the edge, 80% of the memory is a good idea I think
	final long memAvail = Math.round( cudaDevice.getFreeDeviceMemory() * ( percentGPUMem / 100.0 ) );
	final long imgBytes = numPixels() * 4 * 2; // float, two images on the card at once

	final long[] numBlocksDim = net.imglib2.util.Util.int2long( computeNumBlocksDim( memAvail, imgBytes, percentGPUMem, img.numDimensions(), "CUDA-Device " + cudaDevice.getDeviceId() ) );
	final BlockGenerator< Block > generator;

	if ( accurate )
		generator = new BlockGeneratorVariableSizePrecise( numBlocksDim );
	else
		generator = new BlockGeneratorVariableSizeSimple( numBlocksDim );

	final Block[] blocks = generator.divideIntoBlocks( getImgSize( img ), getKernelSize( sigma ) );

	if ( !accurate && blocks.length == 1 && ArrayImg.class.isInstance( img ) )
	{
		IOFunctions.println( "Conovlving image as one single block." );
		long time = System.currentTimeMillis();

		// copy the only directly into the result
		blocks[ 0 ].copyBlock( img, result );
		long copy = System.currentTimeMillis();
		IOFunctions.println( "Copying data took " + ( copy - time ) + "ms" );

		// convolve
		final float[] resultF = ((FloatArray)((ArrayImg< net.imglib2.type.numeric.real.FloatType, ? > )result).update( null ) ).getCurrentStorageArray();
		cudaconvolve.gauss( resultF, getImgSizeInt( result ), sigma, OutOfBounds.EXTEND_BORDER_PIXELS, 0 );
		IOFunctions.println( "Convolution took " + ( System.currentTimeMillis() - copy ) + "ms using device=" + cudaDevice.getDeviceName() + " (id=" + cudaDevice.getDeviceId() + ")" );

		// no copy back required
	}
	else
	{
		final RandomAccessible< net.imglib2.type.numeric.real.FloatType > input;
		
		if ( accurate )
			input = Views.extendMirrorSingle( img );
		else
			input = img;
		
		for( final Block block : blocks )
		{
			//long time = System.currentTimeMillis();
			final ArrayImg< net.imglib2.type.numeric.real.FloatType, FloatArray > imgBlock = ArrayImgs.floats( block.getBlockSize() );

			// copy the block
			block.copyBlock( input, imgBlock );
			//long copy = System.currentTimeMillis();
			//IOFunctions.println( "Copying block took " + ( copy - time ) + "ms" );

			// convolve
			final float[] imgBlockF = ((FloatArray)((ArrayImg< net.imglib2.type.numeric.real.FloatType, ? > )imgBlock).update( null ) ).getCurrentStorageArray();
			cudaconvolve.gauss( imgBlockF, getImgSizeInt( imgBlock ), sigma, OutOfBounds.EXTEND_BORDER_PIXELS, 0 );
			//long convolve = System.currentTimeMillis();
			//IOFunctions.println( "Convolution took " + ( convolve - copy ) + "ms using device=" + cudaDevice.getDeviceName() + " (id=" + cudaDevice.getDeviceId() + ")" );

			// no copy back required
			block.pasteBlock( result, imgBlock );
			//IOFunctions.println( "Pasting block took " + ( System.currentTimeMillis() - convolve ) + "ms" );
		}
	}

	return true;
}
 

@Override
public FloatArray emptyArray(final int entities) {
	return new FloatArray(entities);
}