TrakEM2_/src/main/java/org/janelia/intensity/LinearIntensityMap.java

   1 /**
   2  * License: GPL
   3  *
   4  * This program is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU General Public License 2
   6  * as published by the Free Software Foundation.
   7  *
   8  * This program is distributed in the hope that it will be useful,
   9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11  * GNU General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU General Public License
  14  * along with this program; if not, write to the Free Software
  15  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  16  */
  17 package org.janelia.intensity;
  18
  19 import ij.ImageJ;
  20 import ij.ImagePlus;
  21 import net.imglib2.Cursor;
  22 import net.imglib2.Dimensions;
  23 import net.imglib2.FinalInterval;
  24 import net.imglib2.IterableInterval;
  25 import net.imglib2.RandomAccessible;
  26 import net.imglib2.RandomAccessibleInterval;
  27 import net.imglib2.RealRandomAccessible;
  28 import net.imglib2.img.array.ArrayImg;
  29 import net.imglib2.img.array.ArrayImgs;
  30 import net.imglib2.img.basictypeaccess.array.ByteArray;
  31 import net.imglib2.img.basictypeaccess.array.FloatArray;
  32 import net.imglib2.img.basictypeaccess.array.IntArray;
  33 import net.imglib2.img.basictypeaccess.array.ShortArray;
  34 import net.imglib2.img.display.imagej.ImageJFunctions;
  35 import net.imglib2.interpolation.InterpolatorFactory;
  36 import net.imglib2.interpolation.randomaccess.NLinearInterpolatorFactory;
  37 import net.imglib2.interpolation.randomaccess.NearestNeighborInterpolatorFactory;
  38 import net.imglib2.realtransform.RealViews;
  39 import net.imglib2.realtransform.Scale;
  40 import net.imglib2.realtransform.Translation;
  41 import net.imglib2.type.numeric.ARGBType;
  42 import net.imglib2.type.numeric.NumericType;
  43 import net.imglib2.type.numeric.RealType;
  44 import net.imglib2.type.numeric.integer.UnsignedByteType;
  45 import net.imglib2.type.numeric.integer.UnsignedShortType;
  46 import net.imglib2.type.numeric.real.DoubleType;
  47 import net.imglib2.type.numeric.real.FloatType;
  48 import net.imglib2.view.Views;
  49 import net.imglib2.view.composite.CompositeIntervalView;
  50 import net.imglib2.view.composite.RealComposite;
  51
  52 /**
  53  * Transfer intensities by a linear function <em>y</em>=<em>ax</em>+<em>b</em>
  54  * with the coefficients <em>a</em> and <em>b</em> being stored as a
  55  * {@link RealComposite RealComposite&lt;T&gt;} 2d-vector
  56  * (<em>a</em>, <em>b</em>) in a map.  The
  57  * map is a {@link RealRandomAccessible} and the defined interval as passed as
  58  * an independent parameter.  If the coefficients are passed as a raster, the
  59  * interval is that of a raster.  The map is scaled such that it applies to
  60  * the interval of the input image.
  61  *
  62  * @author Stephan Saalfeld <saalfelds@janelia.hhmi.org>
  63  */
  64 public class LinearIntensityMap< T extends RealType< T > >
  65 {
  66         static public enum Interpolation{ NN, NL };
  67
  68         final static private < T extends RealType< T > >InterpolatorFactory< RealComposite< T >, RandomAccessible< RealComposite< T > > > interpolatorFactory( final Interpolation interpolation )
  69         {
  70                 switch ( interpolation )
  71                 {
  72                 case NN:
  73                         return new NearestNeighborInterpolatorFactory< RealComposite< T > >();
  74                 default:
  75                         return new NLinearInterpolatorFactory< RealComposite< T > >();
  76                 }
  77         }
  78
  79         final protected Dimensions dimensions;
  80         final protected Translation translation;
  81         final protected RealRandomAccessible< RealComposite< T > > coefficients;
  82
  83         final protected InterpolatorFactory< RealComposite< T >, RandomAccessible< RealComposite< T > > > interpolatorFactory;
  84
  85         public LinearIntensityMap( final RandomAccessibleInterval< T > source, final InterpolatorFactory< RealComposite< T >, RandomAccessible< RealComposite< T > > > interpolatorFactory )
  86         {
  87                 this.interpolatorFactory = interpolatorFactory;
  88                 final CompositeIntervalView< T, RealComposite< T > > collapsedSource = Views.collapseReal( source );
  89                 dimensions = new FinalInterval( collapsedSource );
  90                 final double[] shift = new double[ dimensions.numDimensions() ];
  91                 for ( int d = 0; d < shift.length; ++d )
  92                         shift[ d ] = 0.5;
  93                 translation = new Translation( shift );
  94
  95                 final RandomAccessible< RealComposite< T > > extendedCollapsedSource = Views.extendBorder( collapsedSource );
  96                 coefficients = Views.interpolate( extendedCollapsedSource, interpolatorFactory );
  97         }
  98
  99         public LinearIntensityMap( final RandomAccessibleInterval< T > source )
 100         {
 101                 this( source, new NLinearInterpolatorFactory< RealComposite< T > >() );
 102         }
 103
 104         public LinearIntensityMap( final RandomAccessibleInterval< T > source, final Interpolation interpolation )
 105         {
 106                 this( source, LinearIntensityMap.< T >interpolatorFactory( interpolation ) );
 107         }
 108
 109         @SuppressWarnings( { "rawtypes", "unchecked" } )
 110         public < S extends NumericType< S > > void run( final RandomAccessibleInterval< S > image )
 111         {
 112                 assert image.numDimensions() == dimensions.numDimensions() : "Number of dimensions do not match.";
 113
 114                 final double[] s = new double[ dimensions.numDimensions() ];
 115                 for ( int d = 0; d < s.length; ++d )
 116                         s[ d ] = image.dimension( d ) / dimensions.dimension( d );
 117                 final Scale scale = new Scale( s );
 118
 119                 System.out.println( "translation-n " + translation.numDimensions() );
 120
 121                 final RandomAccessibleInterval< RealComposite< T > > stretchedCoefficients =
 122                                 Views.offsetInterval(
 123                                                 Views.raster(
 124                                                                 RealViews.transform(
 125                                                                                 RealViews.transform(
 126                                                                                                 coefficients,
 127                                                                                                 translation ),
 128                                                                                 scale ) ),
 129                                                 image );
 130
 131                 /* decide on type which mapping to use */
 132                 final S t = image.randomAccess().get();
 133
 134                 if ( ARGBType.class.isInstance( t ) )
 135                         mapARGB( Views.flatIterable( ( RandomAccessibleInterval< ARGBType > )image ), Views.flatIterable( stretchedCoefficients ) );
 136                 else if ( RealComposite.class.isInstance( t ) )
 137                         mapComposite( Views.flatIterable( ( RandomAccessibleInterval )image ), Views.flatIterable( stretchedCoefficients ) );
 138                 else if ( RealType.class.isInstance( t ) )
 139                 {
 140                         final RealType< ? > r = ( RealType )t;
 141                         if ( r.getMinValue() > -Double.MAX_VALUE || r.getMaxValue() < Double.MAX_VALUE )
 142 //                          TODO Bug in javac does not enable cast from RandomAccessibleInterval< S > to RandomAccessibleInterval< RealType >, remove when fixed
 143                                 mapCrop( Views.flatIterable( ( RandomAccessibleInterval< RealType > )( Object )image ), Views.flatIterable( stretchedCoefficients ) );
 144                         else
 145 //              TODO Bug in javac does not enable cast from RandomAccessibleInterval< S > to RandomAccessibleInterval< RealType >, remove when fixed
 146                                 map( Views.flatIterable( ( RandomAccessibleInterval< RealType > )( Object )image ), Views.flatIterable( stretchedCoefficients ) );
 147                 }
 148
 149         }
 150
 151         final static protected < S extends RealType< S >, T extends RealType< T > > void map(
 152                         final IterableInterval< S > image,
 153                         final IterableInterval< RealComposite< T > > coefficients )
 154         {
 155                 final Cursor< S > cs = image.cursor();
 156                 final Cursor< RealComposite< T > > ct = coefficients.cursor();
 157
 158                 while ( cs.hasNext() )
 159                 {
 160                         final S s = cs.next();
 161                         final RealComposite< T > t = ct.next();
 162                         s.setReal( s.getRealDouble() * t.get( 0 ).getRealDouble() + t.get( 1 ).getRealDouble() );
 163                 }
 164         }
 165
 166         final static protected < S extends RealType< S >, T extends RealType< T > > void mapCrop(
 167                         final IterableInterval< S > image,
 168                         final IterableInterval< RealComposite< T > > coefficients )
 169         {
 170                 final Cursor< S > cs = image.cursor();
 171                 final Cursor< RealComposite< T > > ct = coefficients.cursor();
 172                 final S firstValue = cs.next();
 173                 final double minS = firstValue.getMinValue();
 174                 final double maxS = firstValue.getMaxValue();
 175
 176                 while ( cs.hasNext() )
 177                 {
 178                         final S s = cs.next();
 179                         final RealComposite< T > t = ct.next();
 180
 181                         s.setReal( Math.max( minS, Math.min( maxS, s.getRealDouble() * t.get( 0 ).getRealDouble() + t.get( 1 ).getRealDouble() ) ) );
 182                 }
 183         }
 184
 185         final static protected < S extends RealType< S >, T extends RealType< T > > void mapComposite(
 186                         final IterableInterval< RealComposite< S > > image,
 187                         final IterableInterval< RealComposite< T > > coefficients )
 188         {
 189                 final Cursor< RealComposite< S > > cs = image.cursor();
 190                 final Cursor< RealComposite< T > > ct = coefficients.cursor();
 191
 192                 while ( cs.hasNext() )
 193                 {
 194                         final RealComposite< S > c = cs.next();
 195                         final RealComposite< T > t = ct.next();
 196
 197                         for ( final S s : c )
 198                                 s.setReal( s.getRealDouble() * t.get( 0 ).getRealDouble() + t.get( 1 ).getRealDouble() );
 199                 }
 200         }
 201
 202         final static protected < T extends RealType< T > > void mapARGB(
 203                         final IterableInterval< ARGBType > image,
 204                         final IterableInterval< RealComposite< T > > coefficients )
 205         {
 206                 final Cursor< ARGBType > cs = image.cursor();
 207                 final Cursor< RealComposite< T > > ct = coefficients.cursor();
 208
 209                 while ( cs.hasNext() )
 210                 {
 211                         final RealComposite< T > t = ct.next();
 212                         final double alpha = t.get( 0 ).getRealDouble();
 213                         final double beta = t.get( 1 ).getRealDouble();
 214
 215                         final ARGBType s = cs.next();
 216                         final int argb = s.get();
 217                         final int a = ( ( argb >> 24 ) & 0xff );
 218                         final double r = ( ( argb >> 16 ) & 0xff ) * alpha + beta;
 219                         final double g = ( ( argb >> 8 ) & 0xff ) * alpha + beta;
 220                         final double b = ( argb & 0xff ) * alpha + beta;
 221
 222                         s.set(
 223                                         ( a << 24 ) |
 224                                         ( ( r < 0 ? 0 : r > 255 ? 255 : ( int )( r + 0.5 ) ) << 16 ) |
 225                                         ( ( g < 0 ? 0 : g > 255 ? 255 : ( int )( g + 0.5 ) ) << 8 ) |
 226                                         ( b < 0 ? 0 : b > 255 ? 255 : ( int )( b + 0.5 ) ) );
 227                 }
 228         }
 229
 230         public static void main( final String[] args )
 231         {
 232                 new ImageJ();
 233
 234                 final double[] coefficients = new double[]{
 235                                 0, 2, 4, 8,
 236                                 1, 1, 1, 1,
 237                                 1, 10, 5, 1,
 238                                 1, 1, 1, 1,
 239
 240                                 0, 10, 20, 30,
 241                                 40, 50, 60, 70,
 242                                 80, 90, 100, 110,
 243                                 120, 130, 140, 150
 244                 };
 245
 246                 final LinearIntensityMap< DoubleType > transform = new LinearIntensityMap< DoubleType >( ArrayImgs.doubles( coefficients, 4, 4, 2 ) );
 247
 248                 //final ImagePlus imp = new ImagePlus( "http://upload.wikimedia.org/wikipedia/en/2/24/Lenna.png" );
 249                 final ImagePlus imp1 = new ImagePlus( "http://fly.mpi-cbg.de/~saalfeld/Pictures/norway.jpg");
 250
 251                 final ArrayImg< FloatType, FloatArray > image1 = ArrayImgs.floats( ( float[] )imp1.getProcessor().convertToFloatProcessor().getPixels(), imp1.getWidth(), imp1.getHeight() );
 252                 final ArrayImg< UnsignedByteType, ByteArray > image2 = ArrayImgs.unsignedBytes( ( byte[] )imp1.getProcessor().convertToByteProcessor().getPixels(), imp1.getWidth(), imp1.getHeight() );
 253                 final ArrayImg< UnsignedShortType, ShortArray > image3 = ArrayImgs.unsignedShorts( ( short[] )imp1.getProcessor().convertToShortProcessor().getPixels(), imp1.getWidth(), imp1.getHeight() );
 254                 final ArrayImg< ARGBType, IntArray > image4 = ArrayImgs.argbs( ( int[] )imp1.getProcessor().getPixels(), imp1.getWidth(), imp1.getHeight() );
 255
 256                 ImageJFunctions.show( ArrayImgs.doubles( coefficients, 4, 4, 2 ) );
 257
 258                 transform.run( image1 );
 259                 transform.run( image2 );
 260                 transform.run( image3 );
 261                 transform.run( image4 );
 262
 263                 ImageJFunctions.show( image1 );
 264                 ImageJFunctions.show( image2 );
 265                 ImageJFunctions.show( image3 );
 266                 ImageJFunctions.show( image4 );
 267         }
 268 }