mpi/fruitfly/registration/FloatArray2DScaleOctave.java

   1 package mpi.fruitfly.registration;
   2
   3 /**
   4  * single octave of a discrete {@link FloatArray2DScaleSpace}
   5  *
   6  * This class is optimized for the Difference Of Gaussian detector used in
   7  * David Lowe's SIFT-algorithm \citep{Loew04}.
   8  *
   9  * The scale space itself consists of an arbitrary number of octaves.  This
  10  * number is implicitly defined by the minimal image size {@link #MIN_SIZE}.
  11  * Octaves contain overlapping scales of the scalespace.  Thus it is possible
  12  * to execute several operations that depend on adjacent scales within one
  13  * octave.
  14  *
  15  * BibTeX:
  16  * <pre>
  17  * &#64;article{Lowe04,
  18  *   author  = {David G. Lowe},
  19  *   title   = {Distinctive Image Features from Scale-Invariant Keypoints},
  20  *   journal = {International Journal of Computer Vision},
  21  *   year    = {2004},
  22  *   volume  = {60},
  23  *   number  = {2},
  24  *   pages   = {91--110},
  25  * }
  26  * </pre>
  27  *
  28  * License: GPL
  29  *
  30  * This program is free software; you can redistribute it and/or
  31  * modify it under the terms of the GNU General Public License 2
  32  * as published by the Free Software Foundation.
  33  *
  34  * This program is distributed in the hope that it will be useful,
  35  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  36  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  37  * GNU General Public License for more details.
  38  *
  39  * You should have received a copy of the GNU General Public License
  40  * along with this program; if not, write to the Free Software
  41  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  42  *
  43  * @author Stephan Saalfeld <stephan.saalfeld@inf.tu-dresden.de>
  44  * @version 0.1b
  45  */
  46
  47 import mpi.fruitfly.math.datastructures.FloatArray2D;
  48 import mpi.fruitfly.registration.ImageFilter;
  49
  50 public class FloatArray2DScaleOctave
  51 {
  52         public enum State { EMPTY, STUB, COMPLETE }
  53
  54         public State state = State.EMPTY;
  55
  56         public int width = 0;
  57         public int height = 0;
  58
  59         private float K = 2.0f;
  60         private float K_MIN1_INV = 1.0f / ( K - 1.0f );
  61
  62         /**
  63          * steps per octave
  64          *
  65          * an octave consists of STEPS + 3 images to be
  66          */
  67         public int STEPS = 1;
  68
  69         /**
  70          * sigma of gaussian kernels corresponding to the steps of the octave
  71          *
  72          * the first member is the sigma of the gaussian kernel that is assumed to
  73          * be the generating kernel of the first gaussian image instance of the
  74          * octave
  75          */
  76         public float[] SIGMA;
  77 //      public float[] getSigma()
  78 //      {
  79 //              return SIGMA;
  80 //      }
  81
  82         /**
  83          * sigma of gaussian kernels required to create the corresponding gaussian
  84          * image instances from the first one
  85          */
  86         private float[] SIGMA_DIFF;
  87
  88         /**
  89          * 1D gaussian kernels required to create the corresponding gaussian
  90          * image instances from the first one
  91          */
  92         private float[][] KERNEL_DIFF;
  93
  94         /**
  95          * gaussian smoothed images
  96          */
  97         private FloatArray2D[] l;
  98         public FloatArray2D[] getL()
  99         {
 100                 return l;
 101         }
 102         public FloatArray2D getL( int i )
 103         {
 104                 return l[ i ];
 105         }
 106
 107         /**
 108          * scale normalised difference of gaussian images
 109          */
 110         private FloatArray2D[] d;
 111         public FloatArray2D[] getD()
 112         {
 113                 return d;
 114         }
 115         public FloatArray2D getD( int i )
 116         {
 117                 return d[ i ];
 118         }
 119
 120         /**
 121          * gradients of the gaussian smoothed images; 0=>amplitudes; 1=>orientations
 122          */
 123         private FloatArray2D[][] l1;
 124         /**
 125          * get the gradients of the corresponding gaussian image, generates it on
 126          * demand, if not yet available.
 127          *
 128          * @param i index will not be checked for efficiency reasons, so take care
 129          *   that it is within a valid range
 130          *
 131          * @returns reference to the gradients
 132          */
 133         public FloatArray2D[] getL1( int i )
 134         {
 135                 if ( l1[ i ] == null )
 136                 {
 137                         l1[ i ] = ImageFilter.createGradients( l[ i ] );
 138                 }
 139                 return l1[ i ];
 140         }
 141
 142         /**
 143          * Constructor
 144          *
 145          * @param img image being the first gaussian instance of the scale octave
 146          *   img must be a 2d-array of float values in range [0.0f, ..., 1.0f]
 147          * @param initial_sigma inital gaussian sigma
 148          */
 149         public FloatArray2DScaleOctave(
 150                         FloatArray2D img,
 151                         int steps,
 152                         float initial_sigma )
 153         {
 154                 state = State.EMPTY;
 155
 156                 width = img.width;
 157                 height = img.height;
 158
 159                 STEPS = steps;
 160
 161                 K = ( float )Math.pow( 2.0, 1.0 / ( float )STEPS );
 162                 K_MIN1_INV = 1.0f / ( K - 1.0f );
 163
 164                 SIGMA = new float[ STEPS + 3 ];
 165                 SIGMA[ 0 ] = initial_sigma;
 166                 SIGMA_DIFF = new float[ STEPS + 3 ];
 167                 SIGMA_DIFF[ 0 ] = 0.0f;
 168                 KERNEL_DIFF = new float[ STEPS + 3 ][];
 169
 170                 //System.out.println( "sigma[0] = " + SIGMA[ 0 ] + "; sigma_diff[0] = " + SIGMA_DIFF[ 0 ] );
 171
 172                 for ( int i = 1; i < STEPS + 3; ++i )
 173                 {
 174                         SIGMA[ i ] = initial_sigma * ( float )Math.pow( 2.0f, ( float )i / ( float )STEPS );
 175                         SIGMA_DIFF[ i ] = ( float )Math.sqrt( SIGMA[ i ] * SIGMA[ i ] - initial_sigma * initial_sigma );
 176
 177                         //System.out.println( "sigma[" + i + "] = " + SIGMA[ i ] + "; sigma_diff[" + i + "] = " + SIGMA_DIFF[ i ] );
 178
 179                         KERNEL_DIFF[ i ] = ImageFilter.createGaussianKernel1D(
 180                                         SIGMA_DIFF[ i ],
 181                                         true );
 182                 }
 183                 l = new FloatArray2D[ 1 ];
 184                 l[ 0 ] = img;
 185                 d = null;
 186                 l1 = null;
 187         }
 188
 189         /**
 190          * Constructor
 191          *
 192          * faster initialisation with precomputed gaussian kernels
 193          *
 194          * @param img image being the first gaussian instance of the scale octave
 195          * @param initial_sigma inital gaussian sigma
 196          *
 197          */
 198         public FloatArray2DScaleOctave(
 199                         FloatArray2D img,
 200                         float[] sigma,
 201                         float[] sigma_diff,
 202                         float[][] kernel_diff )
 203         {
 204                 state = State.EMPTY;
 205
 206                 width = img.width;
 207                 height = img.height;
 208
 209                 STEPS = sigma.length - 3;
 210
 211                 K = ( float )Math.pow( 2.0, 1.0 / ( float )STEPS );
 212                 K_MIN1_INV = 1.0f / ( K - 1.0f );
 213
 214                 SIGMA = sigma;
 215                 SIGMA_DIFF = sigma_diff;
 216                 KERNEL_DIFF = kernel_diff;
 217
 218                 l = new FloatArray2D[ 1 ];
 219                 l[ 0 ] = img;
 220                 d = null;
 221                 l1 = null;
 222         }
 223
 224         /**
 225          * build only the gaussian image with 2 * INITIAL_SIGMA
 226          *
 227          * Use this method for the partial creation of an octaved scale space
 228          * without creating each scale octave.  Like proposed by Lowe
 229          * \citep{Lowe04}, you can use this image to build the next scale octave.
 230          * Taking every second pixel of this image, you get a gaussian  image with
 231          * INITIAL_SIGMA of the half image size.
 232          */
 233         public void buildStub()
 234         {
 235                 FloatArray2D img = l[ 0 ];
 236                 l = new FloatArray2D[ 2 ];
 237                 l[ 0 ] = img;
 238                 l[ 1 ] = ImageFilter.convolveSeparable( l[ 0 ], KERNEL_DIFF[ STEPS ], KERNEL_DIFF[ STEPS ] );
 239
 240                 state = State.STUB;
 241         }
 242
 243
 244         /**
 245          * build the scale octave
 246          */
 247         public boolean build()
 248         {
 249                 FloatArray2D img = l[ 0 ];
 250                 FloatArray2D img2;
 251                 if ( state == State.STUB )
 252                 {
 253                         img2 = l[ 1 ];
 254                         l = new FloatArray2D[ STEPS + 3 ];
 255                         l[ STEPS ] = img2;
 256                 }
 257                 else l = new FloatArray2D[ STEPS + 3 ];
 258                 l[ 0 ] = img;
 259                 for ( int i = 1; i < SIGMA_DIFF.length; ++i )
 260                 {
 261                         if ( state == State.STUB && i == STEPS ) continue;
 262                         // use precomputed kernels
 263                         l[ i ] = ImageFilter.convolveSeparable( l[ 0 ], KERNEL_DIFF[ i ], KERNEL_DIFF[ i ] );
 264                         //l[ i ] = ImageFilter.computeGaussian( l[ 0 ], SIGMA_DIFF[ i ] );
 265                 }
 266                 d = new FloatArray2D[ STEPS + 2 ];
 267                 for ( int i = 0; i < d.length; ++i )
 268                 {
 269                         d[ i ] = new FloatArray2D( l[ i ].width, l[ i ].height );
 270                         int j = i + 1;
 271                         for ( int k = 0; k < l[ i ].data.length; ++k )
 272                         {
 273                                 d[ i ].data[ k ] = ( l[ j ].data[ k ] - l[ i ].data[ k ] ) * K_MIN1_INV;
 274                         }
 275                 }
 276                 l1 = new FloatArray2D[ STEPS + 3 ][];
 277                 for ( int i = 0; i < l1.length; ++i )
 278                 {
 279                         l1[ i ] = null;
 280                 }
 281
 282                 state = State.COMPLETE;
 283
 284                 return true;
 285         }
 286
 287         /**
 288          * clear the scale octave to save memory
 289          */
 290         public void clear()
 291         {
 292                 this.state = State.EMPTY;
 293                 this.d = null;
 294                 this.l = null;
 295                 this.l1 = null;
 296         }
 297
 298
 299         /**
 300          * downsample {@link src} by simply using every second pixel into
 301          * {@link dst}
 302          *
 303          * For efficiency reasons, the dimensions of {@link dst} are not checked,
 304          * that is, you have to take care, that
 305          * dst.width == src.width / 2 + src.width % 2 &&
 306          * dst.height == src.height / 2 + src.height % 2 .
 307          *
 308          * @param src the source image
 309          * @param dst destination image
 310          */
 311         public static void downsample( FloatArray2D src, FloatArray2D dst )
 312         {
 313                 int ws = 2 * src.width;
 314                 int rs = 0;
 315                 for ( int r = 0; r < dst.data.length; r += dst.width )
 316                 {
 317                         int xs = 0;
 318                         for ( int x = 0; x < dst.width; ++x )
 319                         {
 320                                 dst.data[ r + x ] = src.data[ rs + xs ];
 321                                 xs += 2;
 322                         }
 323                         rs += ws;
 324                 }
 325         }
 326
 327         /**
 328          * upsample {@link src} by linearly interpolating into {@link dst}
 329          *
 330          * For efficiency reasons, the dimensions of {@link dst} are not checked,
 331          * that is, you have to take care, that
 332          * src.width == dst.width / 2 + dst.width % 2 &&
 333          * src.height == dst.height / 2 + dst.height % 2 .
 334          *
 335          * @param src the source image
 336          * @param dst destination image
 337          */
 338         public static void upsample( FloatArray2D src, FloatArray2D dst )
 339         {
 340                 int rdw = 2 * dst.width;
 341                 int rd1 = rdw;
 342                 int rd2 = dst.width;
 343                 int xd1 = 2;
 344                 int xd2 = 1;
 345                 dst.data[ 0 ] = src.data[ 0 ];
 346                 for ( int xs1 = 1; xs1 < src.width; ++xs1 )
 347                 {
 348                         int xs2 = xs1 - 1;
 349                         dst.data[ xd1 ] = src.data[ xs1 ];
 350                         dst.data[ xd2 ] = ( src.data[ xs1 ] + src.data[ xs2 ] ) / 2.0f;
 351                         xd1 += 2;
 352                         xd2 += 2;
 353                 }
 354                 for ( int rs1 = src.width; rs1 < src.data.length; rs1 += src.width )
 355                 {
 356                         int rs2 = rs1 - src.width;
 357                         xd1 = 2;
 358                         xd2 = 1;
 359                         dst.data[ rd1 ] = src.data[ rs1 ];
 360                         dst.data[ rd2 ] = ( src.data[ rs1 ] + src.data[ rs2 ] ) / 2;
 361
 362                         for ( int xs1 = 1; xs1 < src.width; ++xs1 )
 363                         {
 364                                 int xs2 = xs1 - 1;
 365                                 dst.data[ rd1 + xd1 ] = src.data[ rs1 + xs1 ];
 366                                 dst.data[ rd1 + xd2 ] = ( src.data[ rs1 + xs1 ] + src.data[ rs1 + xs2 ] ) / 2.0f;
 367                                 dst.data[ rd2 + xd1 ] = ( src.data[ rs1 + xs1 ] + src.data[ rs2 + xs1 ] ) / 2.0f;
 368                                 dst.data[ rd2 + xd2 ] = ( src.data[ rs1 + xs1 ] + src.data[ rs2 + xs2 ] ) / 2.0f;
 369                                 xd1 += 2;
 370                                 xd2 += 2;
 371                         }
 372                         rd1 += rdw;
 373                         rd2 += rdw;
 374                 }
 375                 if ( dst.height % 2 == 0 )
 376                 {
 377                         rd1 = dst.data.length - dst.width;
 378                         rd2 = rd1 - dst.width;
 379                         for ( xd1 = 0; xd1 < dst.width; ++xd1 )
 380                         {
 381                                 dst.data[ rd1 + xd1 ] = dst.data[ rd2 + xd1 ];
 382                         }
 383                 }
 384                 if ( dst.height % 2 == 0 )
 385                 {
 386                         xd1 = dst.width - 1;
 387                         xd2 = dst.width - 2;
 388                         for ( rd1 = 0; rd1 < dst.data.length; rd1 += dst.width )
 389                         {
 390                                 dst.data[ rd1 + xd1 ] = dst.data[ rd1 + xd2 ];
 391                         }
 392                 }
 393         }
 394 }