mpi/fruitfly/registration/Match.java

   1 package mpi.fruitfly.registration;
   2
   3 import static mpi.fruitfly.math.General.sq;
   4
   5 import Jama.EigenvalueDecomposition;
   6 import Jama.Matrix;
   7
   8 import java.util.Collection;
   9
  10 public class Match {
  11
  12         // coordinates
  13         final public float[] p1;
  14         final public float[] p2;
  15
  16         // weight
  17         final public float w1;
  18         final public float w2;
  19
  20         // orientation TODO not needed for visualisation, remove it
  21 //      final public float o1;
  22 //      final public float o2;
  23
  24         public Match( float[] p1, float[] p2 )
  25         {
  26                 this.p1 = p1;
  27                 this.p2 = p2;
  28                 w1 = w2 = 1.0f;
  29 //              o1 = o2 = 0.0f;
  30         }
  31
  32         public Match( float[] p1, float[] p2, float w1, float w2 )
  33         {
  34                 this.p1 = p1;
  35                 this.p2 = p2;
  36                 this.w1 = w1;
  37                 this.w2 = w2;
  38 //              o1 = o2 = 0.0f;
  39         }
  40
  41 //      public Match( float[] p1, float[] p2, float w1, float w2, float o1, float o2 )
  42 //      {
  43 //              this.p1 = p1;
  44 //              this.p2 = p2;
  45 //              this.w1 = w1;
  46 //              this.w2 = w2;
  47 //              this.o1 = o1;
  48 //              this.o2 = o2;
  49 //
  50 //      }
  51
  52         /**
  53          * estimate the covariance, main components (eigenvalues) and mean values
  54          * of the spatial distribution of the matches
  55          */
  56         static public void covariance(
  57                         Collection< Match > matches,
  58                         double[] covariance_P1,         //!< xx_1, xy_1, yy_1
  59                         double[] covariance_P2,         //!< xx_2, xy_2, yy_2
  60                         double[] o_P1,                          //!< x0_1, y0_1
  61                         double[] o_P2,                          //!< x0_2, y0_2
  62                         double[] eigenvalues_P1,        //!< E_1, e_1
  63                         double[] eigenvalues_P2 ,       //!< E_1, e_1
  64                         double[] eigenvector_P1,
  65                         double[] eigenvector_P2
  66                         )
  67         {
  68                 // estimate the covariance matrix of the matching locations in both
  69                 // images
  70                 //---------------------------------------------------------------------
  71
  72                 // first the mean values
  73                 o_P1[ 0 ] = 0;
  74                 o_P1[ 1 ] = 0;
  75                 o_P2[ 0 ] = 0;
  76                 o_P2[ 1 ] = 0;
  77                 for ( Match m : matches )
  78                 {
  79                         o_P1[ 0 ] += ( double )m.p1[ 0 ];
  80                         o_P1[ 1 ] += ( double )m.p1[ 1 ];
  81                         o_P2[ 0 ] += ( double )m.p2[ 0 ];
  82                         o_P2[ 1 ] += ( double )m.p2[ 1 ];
  83                 }
  84                 o_P1[ 0 ] /= matches.size();
  85                 o_P1[ 1 ] /= matches.size();
  86                 o_P2[ 0 ] /= matches.size();
  87                 o_P2[ 1 ] /= matches.size();
  88
  89                 // now the covariances cxx, cyy and cxy
  90                 covariance_P1[ 0 ] = 0;
  91                 covariance_P1[ 1 ] = 0;
  92                 covariance_P1[ 2 ] = 0;
  93                 covariance_P2[ 0 ] = 0;
  94                 covariance_P2[ 1 ] = 0;
  95                 covariance_P2[ 2 ] = 0;
  96
  97                 for ( Match m : matches )
  98                 {
  99                         double dx1 = m.p1[ 0 ] - o_P1[ 0 ];
 100                         double dy1 = m.p1[ 1 ] - o_P1[ 1 ];
 101                         double dx2 = m.p2[ 0 ] - o_P2[ 0 ];
 102                         double dy2 = m.p2[ 1 ] - o_P2[ 1 ];
 103                         covariance_P1[ 0 ] += sq( dx1 );
 104                         covariance_P1[ 1 ] += dx1 * dy1;
 105                         covariance_P1[ 2 ] += sq( dy1 );
 106                         covariance_P2[ 0 ] += sq( dx2 );
 107                         covariance_P2[ 1 ] += dx2 * dy2;
 108                         covariance_P2[ 2 ] += sq( dy2 );
 109                 }
 110                 covariance_P1[ 0 ] /= matches.size();
 111                 covariance_P1[ 1 ] /= matches.size();
 112                 covariance_P1[ 2 ] /= matches.size();
 113                 covariance_P2[ 0 ] /= matches.size();
 114                 covariance_P2[ 1 ] /= matches.size();
 115                 covariance_P2[ 2 ] /= matches.size();
 116
 117                 // estimate the main components being the eigenvalues of the
 118                 // covariance matrix
 119                 /*
 120                 double sqrt1 = Math.sqrt( sq( covariance_P1[ 0 ] + covariance_P1[ 2 ] ) / 4.0 + covariance_P1[ 1 ] * covariance_P1[ 1 ] + covariance_P1[ 0 ] * covariance_P1[ 2 ] );
 121                 double sqrt2 = Math.sqrt( sq( covariance_P2[ 0 ] + covariance_P2[ 2 ] ) / 4.0 + covariance_P2[ 1 ] * covariance_P2[ 1 ] + covariance_P2[ 0 ] * covariance_P2[ 2 ] );
 122
 123                 double p1 = ( covariance_P1[ 0 ] + covariance_P1[ 2 ] ) / 2.0;
 124                 double p2 = ( covariance_P2[ 0 ] + covariance_P2[ 2 ] ) / 2.0;
 125
 126                 eigenvalues_P1[ 0 ] = p1 + sqrt1;
 127                 eigenvalues_P1[ 1 ] = p1 - sqrt1;
 128                 eigenvalues_P2[ 0 ] = p2 + sqrt2;
 129                 eigenvalues_P2[ 1 ] = p2 - sqrt2;
 130
 131                 if ( eigenvalues_P1[ 0 ] < eigenvalues_P1[ 1 ] )
 132                 {
 133                         double t = eigenvalues_P1[ 0 ];
 134                         eigenvalues_P1[ 0 ] = eigenvalues_P1[ 1 ];
 135                         eigenvalues_P1[ 1 ] = t;
 136                 }
 137                 if ( eigenvalues_P2[ 0 ] < eigenvalues_P2[ 1 ] )
 138                 {
 139                         double t = eigenvalues_P2[ 0 ];
 140                         eigenvalues_P2[ 0 ] = eigenvalues_P2[ 1 ];
 141                         eigenvalues_P2[ 1 ] = t;
 142                 }
 143                 */
 144                 EigenvalueDecomposition evd1 =
 145                         new EigenvalueDecomposition(
 146                                         new Matrix(
 147                                                         new double[][]{
 148                                                                         { covariance_P1[ 0 ], covariance_P1[ 1 ] },
 149                                                                         { covariance_P1[ 1 ], covariance_P1[ 2 ] } } ) );
 150                 EigenvalueDecomposition evd2 =
 151                         new EigenvalueDecomposition(
 152                                         new Matrix(
 153                                                         new double[][]{
 154                                                                         { covariance_P2[ 0 ], covariance_P2[ 1 ] },
 155                                                                         { covariance_P2[ 1 ], covariance_P2[ 2 ] } } ) );
 156                 double[][] e = evd1.getD().getArray();
 157                 eigenvalues_P1[ 0 ] = e[ 0 ][ 0 ];
 158                 eigenvalues_P1[ 1 ] = e[ 1 ][ 1 ];
 159                 e = evd2.getD().getArray();
 160                 eigenvalues_P2[ 0 ] = e[ 0 ][ 0 ];
 161                 eigenvalues_P2[ 1 ] = e[ 1 ][ 1 ];
 162                 double[][] v = evd1.getV().getArray();
 163                 eigenvector_P1[ 0 ] = v[ 0 ][ 0 ];
 164                 eigenvector_P1[ 1 ] = v[ 1 ][ 0 ];
 165                 eigenvector_P1[ 2 ] = v[ 0 ][ 1 ];
 166                 eigenvector_P1[ 3 ] = v[ 1 ][ 1 ];
 167                 v = evd2.getV().getArray();
 168                 eigenvector_P2[ 0 ] = v[ 0 ][ 0 ];
 169                 eigenvector_P2[ 1 ] = v[ 1 ][ 0 ];
 170                 eigenvector_P2[ 2 ] = v[ 0 ][ 1 ];
 171                 eigenvector_P2[ 3 ] = v[ 1 ][ 1 ];
 172
 173         }
 174 }