mpi/fruitfly/registration/TRModel2D.java

   1 package mpi.fruitfly.registration;
   2
   3 import ini.trakem2.utils.Utils;
   4
   5 import java.util.Iterator;
   6 import java.util.ArrayList;
   7 import java.util.Collection;
   8 import java.util.List;
   9 import java.util.Random;
  10 import java.awt.geom.AffineTransform;
  11
  12 public class TRModel2D extends Model {
  13
  14         static final public int MIN_SET_SIZE = 2;
  15
  16         final private AffineTransform affine = new AffineTransform();
  17         public AffineTransform getAffine() { return affine; }
  18
  19         @Override
  20         public float[] apply( float[] point )
  21         {
  22                 float[] transformed = new float[ 2 ];
  23                 affine.transform( point, 0, transformed, 0, 1 );
  24                 return transformed;
  25         }
  26
  27         @Override
  28         public void applyInPlace( float[] point )
  29         {
  30                 affine.transform( point, 0, point, 0, 1 );
  31         }
  32
  33         @Override
  34         public float[] applyInverse( float[] point )
  35         {
  36                 // the brilliant java.awt.geom.AffineTransform implements transform for float[] but inverseTransform for double[] only...
  37                 double[] double_point = new double[]{ point[ 0 ], point[ 1 ] };
  38                 double[] transformed = new double[ 2 ];
  39                 try
  40                 {
  41                         affine.inverseTransform( double_point, 0, transformed, 0, 1 );
  42                 }
  43                 catch ( Exception e )
  44                 {
  45                         System.err.println( "Noninvertible transformation." );
  46                 }
  47                 return new float[]{ ( float )transformed[ 0 ], ( float )transformed[ 1 ] };
  48         }
  49
  50         @Override
  51         public void applyInverseInPlace( float[] point )
  52         {
  53                 float[] temp_point = applyInverse( point );
  54                 point[ 0 ] = temp_point[ 0 ];
  55                 point[ 1 ] = temp_point[ 1 ];
  56         }
  57
  58         @Override
  59         public boolean fit( PointMatch[] min_matches )
  60         {
  61                 PointMatch m1 = min_matches[ 0 ];
  62                 PointMatch m2 = min_matches[ 1 ];
  63
  64                 float[] m1_p1 = m1.getP1().getL();
  65                 float[] m2_p1 = m2.getP1().getL();
  66                 float[] m1_p2 = m1.getP2().getW();
  67                 float[] m2_p2 = m2.getP2().getW();
  68
  69                 float x1 = m2_p1[ 0 ] - m1_p1[ 0 ];
  70                 float y1 = m2_p1[ 1 ] - m1_p1[ 1 ];
  71                 float x2 = m2_p2[ 0 ] - m1_p2[ 0 ];
  72                 float y2 = m2_p2[ 1 ] - m1_p2[ 1 ];
  73                 float l1 = ( float )Math.sqrt( x1 * x1 + y1 * y1 );
  74                 float l2 = ( float )Math.sqrt( x2 * x2 + y2 * y2 );
  75
  76                 x1 /= l1;
  77                 x2 /= l2;
  78                 y1 /= l1;
  79                 y2 /= l2;
  80
  81                 //! unrotate (x2,y2)^T to (x1,y1)^T = (1,0)^T getting the sinus and cosinus of the rotation angle
  82                 float cos = x1 * x2 + y1 * y2;
  83                 float sin = x1 * y2 - y1 * x2;
  84
  85                 //m.alpha = atan2( y, x );
  86
  87                 float tx = m1_p2[ 0 ] - cos * m1_p1[ 0 ] + sin * m1_p1[ 1 ];
  88                 float ty = m1_p2[ 1 ] - sin * m1_p1[ 0 ] - cos * m1_p1[ 1 ];
  89                 affine.setTransform( cos, sin, -sin, cos, tx, ty );
  90
  91 //              System.out.println( this );
  92
  93                 return true;
  94         }
  95
  96         @Override
  97         public String toString()
  98         {
  99                 return ( "[3,3](" + affine + ") " + error );
 100         }
 101
 102         public void minimize( Collection< PointMatch > matches )
 103         {
 104                 // center of mass:
 105                 float xo1 = 0, yo1 = 0;
 106                 float xo2 = 0, yo2 = 0;
 107                 // Implementing Johannes Schindelin's squared error minimization formula
 108                 // tan(angle) = Sum(x1*y1 + x2y2) / Sum(x1*y2 - x2*y1)
 109                 int length = matches.size();
 110                 // 1 - compute centers of mass, for displacement and origin of rotation
 111
 112                 if (0 == length) return;
 113
 114                 for ( PointMatch m : matches )
 115                 {
 116                         float[] m_p1 = m.getP1().getL();
 117                         float[] m_p2 = m.getP2().getW();
 118
 119                         xo1 += m_p1[ 0 ];
 120                         yo1 += m_p1[ 1 ];
 121                         xo2 += m_p2[ 0 ];
 122                         yo2 += m_p2[ 1 ];
 123                 }
 124                 xo1 /= length;
 125                 yo1 /= length;
 126                 xo2 /= length;
 127                 yo2 /= length;
 128
 129                 float dx = xo1 - xo2; // reversed, because the second will be moved relative to the first
 130                 float dy = yo1 - yo2;
 131                 float sum1 = 0, sum2 = 0;
 132                 float x1, y1, x2, y2;
 133                 for ( PointMatch m : matches )
 134                 {
 135                         float[] m_p1 = m.getP1().getL();
 136                         float[] m_p2 = m.getP2().getW();
 137
 138                         // make points local to the center of mass of the first landmark set
 139                         x1 = m_p1[ 0 ] - xo1; // x1
 140                         y1 = m_p1[ 1 ] - yo1; // x2
 141                         x2 = m_p2[ 0 ] - xo2 + dx; // y1
 142                         y2 = m_p2[ 1 ] - yo2 + dy; // y2
 143                         sum1 += x1 * y2 - y1 * x2; //   x1 * y2 - x2 * y1 // assuming p1 is x1,x2, and p2 is y1,y2
 144                         sum2 += x1 * x2 + y1 * y2; //   x1 * y1 + x2 * y2
 145                 }
 146                 float angle = ( float )Math.atan2( -sum1, sum2 );
 147
 148                 affine.setToIdentity();
 149                 affine.rotate( -angle, xo2, yo2 );
 150                 affine.translate( -dx, -dy );
 151         }
 152
 153         /**
 154          * change the model a bit
 155          *
 156          * estimates the necessary amount of shaking for each single dimensional
 157          * distance in the set of matches
 158          *
 159          * @param matches point matches
 160          * @param scale gives a multiplicative factor to each dimensional distance (increases the amount of shaking)
 161          * @param center local pivot point
 162          */
 163         final public void shake(
 164                         Collection< PointMatch > matches,
 165                         float scale,
 166                         float[] center )
 167         {
 168                 double xd = 0.0;
 169                 double yd = 0.0;
 170                 double rd = 0.0;
 171
 172                 int num_matches = matches.size();
 173                 if ( num_matches > 0 )
 174                 {
 175                         for ( PointMatch m : matches )
 176                         {
 177                                 float[] m_p1 = m.getP1().getW();
 178                                 float[] m_p2 = m.getP2().getW();
 179
 180                                 xd += Math.abs( m_p1[ 0 ] - m_p2[ 0 ] );;
 181                                 yd += Math.abs( m_p1[ 1 ] - m_p2[ 1 ] );;
 182
 183                                 // shift relative to the center
 184                                 float x1 = m_p1[ 0 ] - center[ 0 ];
 185                                 float y1 = m_p1[ 1 ] - center[ 1 ];
 186                                 float x2 = m_p2[ 0 ] - center[ 0 ];
 187                                 float y2 = m_p2[ 1 ] - center[ 1 ];
 188
 189                                 float l1 = ( float )Math.sqrt( x1 * x1 + y1 * y1 );
 190                                 float l2 = ( float )Math.sqrt( x2 * x2 + y2 * y2 );
 191
 192                                 x1 /= l1;
 193                                 x2 /= l2;
 194                                 y1 /= l1;
 195                                 y2 /= l2;
 196
 197                                 //! unrotate (x1,y1)^T to (x2,y2)^T = (1,0)^T getting the sinus and cosinus of the rotation angle
 198                                 float cos = x1 * x2 + y1 * y2;
 199                                 float sin = y1 * x2 - x1 * y2;
 200
 201                                 rd += Math.abs( Math.atan2( sin, cos ) );
 202                         }
 203                         xd /= matches.size();
 204                         yd /= matches.size();
 205                         rd /= matches.size();
 206
 207                         //System.out.println( rd );
 208                 }
 209
 210                 affine.rotate( rnd.nextGaussian() * ( float )rd * scale, center[ 0 ], center[ 1 ] );
 211         }
 212
 213         /**
 214          * estimate the transformation model for a set of feature correspondences
 215          * containing a high number of outliers using RANSAC
 216          *
 217          * @param candidates set of correspondence candidates
 218          * @param inliers set ot correspondences that fit the finally estimated model if any
 219          * @param iterations number of iterations
 220          * @param epsilon maximally allowed displacement
 221          * @param min_inlier_ratio minimal amount of inliers
 222          *
 223          * @return TRModel2D or null
 224          *
 225          *
 226          * Bibtex reference:
 227          * <pre>
 228          * @article{FischlerB81,
 229          *      author                  = {Martin A. Fischler and Robert C. Bolles},
 230          *      title                   = {Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography},
 231          *      journal                 = {Communications of the ACM},
 232          *      volume                  = {24},
 233          *      number                  = {6},
 234          *      year                    = {1981},
 235          *      pages                   = {381--395},
 236          *      publisher               = {ACM Press},
 237          *      address                 = {New York, NY, USA},
 238          *      issn                    = {0001-0782},
 239          *      doi                             = {http://doi.acm.org/10.1145/358669.358692},
 240          * }
 241          * </pre>
 242          *
 243          */
 244         static public TRModel2D estimateModel(
 245                         List< PointMatch > candidates,
 246                         Collection< PointMatch > inliers,
 247                         int iterations,
 248                         float epsilon,
 249                         float min_inlier_ratio )
 250         {
 251                 inliers.clear();
 252
 253                 if ( candidates.size() < MIN_SET_SIZE )
 254                 {
 255                         System.err.println( candidates.size() + " correspondence candidates are not enough to estimate a model, at least " + TRModel2D.MIN_SET_SIZE + " required." );
 256                         return null;
 257                 }
 258
 259                 TRModel2D model = new TRModel2D();              //!< the final model to be estimated
 260
 261                 int i = 0;
 262                 while ( i < iterations )
 263                 {
 264                         // choose T::MIN_SET_SIZE disjunctive matches randomly
 265                         PointMatch[] min_matches = new PointMatch[ MIN_SET_SIZE ];
 266                         int[] keys = new int[ MIN_SET_SIZE ];
 267
 268                         for ( int j = 0; j < MIN_SET_SIZE; ++j )
 269                         {
 270                                 int key;
 271                                 boolean in_set = false;
 272                                 do
 273                                 {
 274                                         key = ( int )( rnd.nextDouble() * candidates.size() );
 275                                         in_set = false;
 276
 277                                         // check if this key exists already
 278                                         for ( int k = 0; k < j; ++k )
 279                                         {
 280                                                 if ( key == keys[ k ] )
 281                                                 {
 282                                                         in_set = true;
 283                                                         break;
 284                                                 }
 285                                         }
 286                                 }
 287                                 while ( in_set );
 288                                 keys[ j ] = key;
 289                                 min_matches[ j ] = candidates.get( key );
 290                         }
 291
 292                         TRModel2D m = new TRModel2D();
 293                         final ArrayList< PointMatch > temp_inliers = new ArrayList< PointMatch >();
 294                         m.fit( min_matches );
 295                         int num_inliers = 0;
 296                         boolean is_good = m.test( candidates, temp_inliers, epsilon, min_inlier_ratio );
 297                         while ( is_good && num_inliers < temp_inliers.size() )
 298                         {
 299                                 num_inliers = temp_inliers.size();
 300                                 m.minimize( temp_inliers );
 301                                 is_good = m.test( candidates, temp_inliers, epsilon, min_inlier_ratio );
 302                         }
 303                         if (
 304                                         is_good &&
 305                                         m.betterThan( model ) &&
 306                                         temp_inliers.size() >= 3 * MIN_SET_SIZE ) // now at least 6 matches required
 307                         {
 308                                 model = m.clone();
 309                                 inliers.clear();
 310                                 inliers.addAll( temp_inliers );
 311                         }
 312                         ++i;
 313                 }
 314                 if ( inliers.size() == 0 )
 315                         return null;
 316
 317                 return model;
 318         }
 319
 320         /**
 321          * estimate the transformation model for a set of feature correspondences
 322          * containing a high number of outliers using RANSAC
 323          *
 324          * increase the error as long as not more inliers occur
 325          */
 326         static public TRModel2D estimateBestModel(
 327                         List< PointMatch > candidates,
 328                         Collection< PointMatch > inliers,
 329                         float min_epsilon,
 330                         float max_epsilon,
 331                         float min_inlier_ratio )
 332         {
 333                 inliers.clear();
 334                 TRModel2D model = null;
 335                 float epsilon = 0.0f;
 336                 if ( candidates.size() > MIN_SET_SIZE )
 337                 {
 338                         int highest_num_inliers = 0;
 339                         int convergence_count = 0;
 340                         TRModel2D m = null;
 341                         do
 342                         {
 343                                 final ArrayList< PointMatch > temp_inliers = new ArrayList< PointMatch >();
 344                                 epsilon += min_epsilon;
 345                                 // 1000 iterations lead to a probability of < 0.01% that only bad data values were found
 346                                 m = estimateModel(
 347                                                 candidates,                                     //!< point correspondence candidates
 348                                                 temp_inliers,
 349                                                 1000,                                           //!< iterations
 350                                                 epsilon,                                        //!< maximal alignment error for a good point pair when fitting the model
 351                                                 min_inlier_ratio );                             //!< minimal partition (of 1.0) of inliers
 352
 353                                 if ( m != null )
 354                                 {
 355                                         int num_inliers = temp_inliers.size();
 356                                         if ( num_inliers <= highest_num_inliers )
 357                                         {
 358                                                 ++convergence_count;
 359                                         }
 360                                         else
 361                                         {
 362                                                 model = m.clone();
 363                                                 inliers.clear();
 364                                                 inliers.addAll( temp_inliers );
 365                                                 convergence_count = 0;
 366                                                 highest_num_inliers = num_inliers;
 367                                         }
 368                                 }
 369                         }
 370                         while ( ( m == null || convergence_count < 4 ) && epsilon < max_epsilon );
 371                 }
 372                 if ( model == null )
 373                 {
 374                         Utils.log2( "No model found." );
 375                 }
 376                 else
 377                 {
 378                         Utils.log2( "Model with epsilon <= " + epsilon + " for " + inliers.size() + " inliers found." );
 379                         Utils.log2( "  Affine transform: " + model.getAffine().toString() );
 380                 }
 381
 382                 return model;
 383         }
 384
 385
 386         public TRModel2D clone()
 387         {
 388                 TRModel2D trm = new TRModel2D();
 389                 trm.affine.setTransform( affine );
 390                 trm.error = error;
 391                 return trm;
 392         }
 393
 394         public void preConcatenate(TRModel2D model) {
 395                 this.affine.preConcatenate(model.affine);
 396         }
 397
 398         public void concatenate(TRModel2D model) {
 399                 this.affine.concatenate(model.affine);
 400         }
 401 }