lenscorrection/NonLinearTransform.java

   1 /**
   2
   3 Copyright (C) 2008 Verena Kaynig.
   4
   5 This program is free software; you can redistribute it and/or
   6 modify it under the terms of the GNU General Public License
   7 as published by the Free Software Foundation (http://www.gnu.org/licenses/gpl.txt )
   8
   9 This program is distributed in the hope that it will be useful,
  10 but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12 GNU General Public License for more details.
  13
  14 You should have received a copy of the GNU General Public License
  15 along with this program; if not, write to the Free Software
  16 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  17  **/
  18
  19 /* ****************************************************************  *
  20  * Representation of a non linear transform by explicit polynomial
  21  * kernel expansion.
  22  *
  23  * TODO:
  24  *      - make different kernels available
  25  *      - inverse transform for visualization
  26  *  - improve image interpolation
  27  *  - apply and applyInPlace should use precalculated transform?
  28  *    (What about out of image range pixels?)
  29  *
  30  *  Author: Verena Kaynig
  31  *  Kontakt: verena.kaynig@inf.ethz.ch
  32  *
  33  * ****************************************************************  */
  34
  35 package lenscorrection;
  36
  37 import ij.process.ByteProcessor;
  38 import ij.process.ColorProcessor;
  39 import ij.process.FloatProcessor;
  40 import ij.process.ImageProcessor;
  41 import ij.IJ;
  42 import ij.ImagePlus;
  43 import ij.io.FileSaver;
  44
  45 import java.awt.Color;
  46 import java.awt.geom.GeneralPath;
  47 import java.io.BufferedWriter;
  48 import java.lang.Math;
  49
  50 import Jama.Matrix;
  51
  52 import java.io.*;
  53
  54
  55 public class NonLinearTransform implements mpicbg.trakem2.transform.CoordinateTransform{
  56
  57         private double[][] beta = null;
  58         private double[] normMean = null;
  59         private double[] normVar = null;
  60         private double[][][] transField = null;
  61         private int dimension = 0;
  62         private int length = 0;
  63         private int width = 0;
  64         private int height = 0;
  65
  66         public int getDimension(){ return dimension; }
  67         /** Deletes all dimension dependent properties */
  68         public void setDimension( final int dimension )
  69         {
  70                 this.dimension = dimension;
  71                 length = (dimension + 1)*(dimension + 2)/2;
  72
  73                 beta = new double[length][2];
  74                 normMean = new double[length];
  75                 normVar = new double[length];
  76
  77                 for (int i=0; i < length; i++){
  78                         normMean[i] = 0;
  79                         normVar[i] = 1;
  80                 }
  81                 transField = null;
  82                 precalculated = false;
  83         }
  84
  85         private boolean precalculated = false;
  86
  87         public int getMinNumMatches()
  88         {
  89                 return length;
  90         }
  91
  92
  93         public void fit( final double x[][], final double y[][], final double lambda )
  94         {
  95                 final double[][] expandedX = kernelExpandMatrixNormalize( x );
  96
  97                 final Matrix phiX = new Matrix( expandedX, expandedX.length, length );
  98                 final Matrix phiXTransp = phiX.transpose();
  99
 100                 final Matrix phiXProduct = phiXTransp.times( phiX );
 101
 102                 final int l = phiXProduct.getRowDimension();
 103                 final double lambda2 = 2 * lambda;
 104
 105                 for (int i = 0; i < l; ++i )
 106                         phiXProduct.set( i, i, phiXProduct.get( i, i ) + lambda2 );
 107
 108                 final Matrix phiXPseudoInverse = phiXProduct.inverse();
 109                 final Matrix phiXProduct2 = phiXPseudoInverse.times( phiXTransp );
 110                 final Matrix betaMatrix = phiXProduct2.times( new Matrix( y, y.length, 2 ) );
 111
 112                 setBeta( betaMatrix.getArray() );
 113         }
 114
 115         public void estimateDistortion( double hack1[][], double hack2[][], double transformParams[][], double lambda, int w, int h )
 116         {
 117                 beta = new double[ length ][ 2 ];
 118                 normMean = new double[ length ];
 119                 normVar = new double[ length ];
 120
 121                 for ( int i = 0; i < length; i++ )
 122                 {
 123                         normMean[ i ] = 0;
 124                         normVar[ i ] = 1;
 125                 }
 126
 127                 width = w;
 128                 height = h;
 129
 130                 final double expandedX[][] = kernelExpandMatrixNormalize( hack1 );
 131                 final double expandedY[][] = kernelExpandMatrix( hack2 );
 132
 133                 int s = expandedX[ 0 ].length;
 134                 Matrix S1 = new Matrix( 2 * s, 2 * s );
 135                 Matrix S2 = new Matrix( 2 * s, 1 );
 136
 137                 for ( int i = 0; i < expandedX.length; ++i )
 138                 {
 139                         Matrix xk_ij = new Matrix( expandedX[ i ], 1 );
 140                         Matrix xk_ji = new Matrix( expandedY[ i ], 1 );
 141
 142                         Matrix yk1a = xk_ij.minus( xk_ji.times( transformParams[ i ][ 0 ] ) );
 143                         Matrix yk1b = xk_ij.times( 0.0 ).minus( xk_ji.times( -transformParams[ i ][ 2 ] ) );
 144                         Matrix yk2a = xk_ij.times( 0.0 ).minus( xk_ji.times( -transformParams[ i ][ 1 ] ) );
 145                         Matrix yk2b = xk_ij.minus( xk_ji.times( transformParams[ i ][ 3 ] ) );
 146
 147                         Matrix y = new Matrix( 2, 2 * s );
 148                         y.setMatrix( 0, 0, 0, s - 1, yk1a );
 149                         y.setMatrix( 0, 0, s, 2 * s - 1, yk1b );
 150                         y.setMatrix( 1, 1, 0, s - 1, yk2a );
 151                         y.setMatrix( 1, 1, s, 2 * s - 1, yk2b );
 152
 153                         Matrix xk = new Matrix( 2, 2 * expandedX[ 0 ].length );
 154                         xk.setMatrix( 0, 0, 0, s - 1, xk_ij );
 155                         xk.setMatrix( 1, 1, s, 2 * s - 1, xk_ij );
 156
 157                         double[] vals = { hack1[ i ][ 0 ], hack1[ i ][ 1 ] };
 158                         Matrix c = new Matrix( vals, 2 );
 159
 160                         Matrix X = xk.transpose().times( xk ).times( lambda );
 161                         Matrix Y = y.transpose().times( y );
 162
 163                         S1 = S1.plus( Y.plus( X ) );
 164
 165                         double trans1 = ( transformParams[ i ][ 2 ] * transformParams[ i ][ 5 ] - transformParams[ i ][ 0 ] * transformParams[ i ][ 4 ] );
 166                         double trans2 = ( transformParams[ i ][ 1 ] * transformParams[ i ][ 4 ] - transformParams[ i ][ 3 ] * transformParams[ i ][ 5 ] );
 167                         double[] trans = { trans1, trans2 };
 168
 169                         Matrix translation = new Matrix( trans, 2 );
 170                         Matrix YT = y.transpose().times( translation );
 171                         Matrix XC = xk.transpose().times( c ).times( lambda );
 172
 173                         S2 = S2.plus( YT.plus( XC ) );
 174                 }
 175                 Matrix regularize = Matrix.identity( S1.getRowDimension(), S1.getColumnDimension() );
 176                 Matrix beta = new Matrix( S1.plus( regularize.times( 0.001 ) ).inverse().times( S2 ).getColumnPackedCopy(), s );
 177
 178                 setBeta( beta.getArray() );
 179         }
 180
 181         public NonLinearTransform(double[][] b, double[] nm, double[] nv, int d, int w, int h){
 182                 beta = b;
 183                 normMean = nm;
 184                 normVar = nv;
 185                 dimension = d;
 186                 length = (dimension + 1)*(dimension + 2)/2;
 187                 width = w;
 188                 height = h;
 189         }
 190
 191         public NonLinearTransform(int d, int w, int h){
 192                 dimension = d;
 193                 length = (dimension + 1)*(dimension + 2)/2;
 194
 195                 beta = new double[length][2];
 196                 normMean = new double[length];
 197                 normVar = new double[length];
 198
 199                 for (int i=0; i < length; i++){
 200                         normMean[i] = 0;
 201                         normVar[i] = 1;
 202                 }
 203
 204                 width = w;
 205                 height = h;
 206         }
 207
 208         public NonLinearTransform(){};
 209
 210         public NonLinearTransform(String filename){
 211                 this.load(filename);
 212         }
 213
 214         public NonLinearTransform(double[][] coeffMatrix, int w, int h){
 215                 length = coeffMatrix.length;
 216                 beta = new double[length][2];
 217                 normMean = new double[length];
 218                 normVar = new double[length];
 219                 width = w;
 220                 height = h;
 221                 dimension = (int)(-1.5 + Math.sqrt(0.25 + 2*length));
 222
 223                 for(int i=0; i<length; i++){
 224                         beta[i][0] = coeffMatrix[0][i];
 225                         beta[i][1] = coeffMatrix[1][i];
 226                         normMean[i] = coeffMatrix[2][i];
 227                         normVar[i] = coeffMatrix[3][i];
 228                 }
 229         }
 230
 231
 232         //implements mpicbg.trakem2
 233         public void init( String data ) throws NumberFormatException{
 234                 final String[] fields = data.split( " " );
 235                 int c = 0;
 236
 237                 dimension = Integer.parseInt(fields[c]); c++;
 238                 length = Integer.parseInt(fields[c]); c++;
 239
 240                 beta = new double[length][2];
 241                 normMean = new double[length];
 242                 normVar = new double[length];
 243
 244                 if ( fields.length == 4 + 4*length )
 245                 {
 246                         for (int i=0; i < length; i++){
 247                                 beta[i][0] = Double.parseDouble(fields[c]); c++;
 248                                 beta[i][1] = Double.parseDouble(fields[c]); c++;
 249                         }
 250
 251                         //System.out.println("c: " + c);
 252
 253                         for (int i=0; i < length; i++){
 254                                 normMean[i] = Double.parseDouble(fields[c]); c++;
 255                         }
 256
 257                         //System.out.println("c: " + c);
 258
 259                         for (int i=0; i < length; i++){
 260                                 normVar[i] = Double.parseDouble(fields[c]); c++;
 261                         }
 262
 263                         width = Integer.parseInt(fields[c]); c++;
 264                         height = Integer.parseInt(fields[c]); c++;
 265                         //System.out.println("c: " + c);
 266
 267                 }
 268                 else throw new NumberFormatException( "Inappropriate parameters for " + this.getClass().getCanonicalName() );
 269         }
 270
 271
 272
 273         public String toXML(final String indent){
 274                 return new StringBuilder(indent).append("<ict_transform class=\"").append(this.getClass().getCanonicalName()).append("\" data=\"").append(toDataString()).append("\"/>").toString();
 275         }
 276
 277         public String toDataString(){
 278                 String data = "";
 279                 data += Integer.toString(dimension) + " ";
 280                 data += Integer.toString(length) + " ";
 281
 282                 for (int i=0; i < length; i++){
 283                         data += Double.toString(beta[i][0]) + " ";
 284                         data += Double.toString(beta[i][1]) + " ";
 285                 }
 286
 287                 for (int i=0; i < length; i++){
 288                         data += Double.toString(normMean[i]) + " ";
 289                 }
 290
 291                 for (int i=0; i < length; i++){
 292                         data += Double.toString(normVar[i]) + " ";
 293                 }
 294                 data += Integer.toString(width) + " ";
 295                 data += Integer.toString(height) + " ";
 296
 297                 return data;
 298
 299         }
 300
 301         //@Override
 302         public String toString(){ return toDataString(); }
 303
 304         public float[] apply( float[] location ){
 305
 306                 double[] position = {(double) location[0], (double) location[1]};
 307                 double[] featureVector = kernelExpand(position);
 308                 double[] newPosition = multiply(beta, featureVector);
 309
 310                 float[] newLocation = new float[2];
 311                 newLocation[0] = (float) newPosition[0];
 312                 newLocation[1] = (float) newPosition[1];
 313
 314                 return newLocation;
 315         }
 316
 317         public void applyInPlace( float[] location ){
 318                 double[] position = {(double) location[0], (double) location[1]};
 319                 double[] featureVector = kernelExpand(position);
 320                 double[] newPosition = multiply(beta, featureVector);
 321
 322                 location[0] = (float) newPosition[0];
 323                 location[1] = (float) newPosition[1];
 324         }
 325
 326         void precalculateTransfom(){
 327                 transField = new double[width][height][2];
 328                 //double minX = width, minY = height, maxX = 0, maxY = 0;
 329
 330                 for (int x=0; x<width; x++){
 331                         for (int y=0; y<height; y++){
 332                                 double[] position = {x,y};
 333                                 double[] featureVector = kernelExpand(position);
 334                                 double[] newPosition = multiply(beta, featureVector);
 335
 336                                 if ((newPosition[0] < 0) || (newPosition[0] >= width) ||
 337                                                 (newPosition[1] < 0) || (newPosition[1] >= height))
 338                                 {
 339                                         transField[x][y][0] = -1;
 340                                         transField[x][y][1] = -1;
 341                                         continue;
 342                                 }
 343
 344                                 transField[x][y][0] = newPosition[0];
 345                                 transField[x][y][1] = newPosition[1];
 346
 347                                 //minX = Math.min(minX, x);
 348                                 //minY = Math.min(minY, y);
 349                                 //maxX = Math.max(maxX, x);
 350                                 //maxY = Math.max(maxY, y);
 351
 352                         }
 353                 }
 354
 355                 precalculated = true;
 356         }
 357
 358         public double[][] getCoefficients(){
 359                 double[][] coeffMatrix = new double[4][length];
 360
 361                 for(int i=0; i<length; i++){
 362                         coeffMatrix[0][i] = beta[i][0];
 363                         coeffMatrix[1][i] = beta[i][1];
 364                         coeffMatrix[2][i] = normMean[i];
 365                         coeffMatrix[3][i] = normVar[i];
 366
 367                 }
 368                 return coeffMatrix;
 369         }
 370
 371         public void setBeta(double[][] b){
 372                 beta = b;
 373                 //FIXME: test if normMean and normVar are still valid for this beta
 374         }
 375
 376         public void print(){
 377                 System.out.println("beta:");
 378                 for (int i=0; i < beta.length; i++){
 379                         for (int j=0; j < beta[i].length; j++){
 380                                 System.out.print(beta[i][j]);
 381                                 System.out.print(" ");
 382                         }
 383                         System.out.println();
 384                 }
 385
 386                 System.out.println("normMean:");
 387                 for (int i=0; i < normMean.length; i++){
 388                         System.out.print(normMean[i]);
 389                         System.out.print(" ");
 390                 }
 391
 392                 System.out.println("normVar:");
 393                 for (int i=0; i < normVar.length; i++){
 394                         System.out.print(normVar[i]);
 395                         System.out.print(" ");
 396                 }
 397
 398                 System.out.println("Image size:");
 399                 System.out.println("width: " + width + " height: " + height);
 400
 401                 System.out.println();
 402
 403         }
 404
 405         public void save( final String filename )
 406         {
 407                 try{
 408                         BufferedWriter out = new BufferedWriter(
 409                                         new OutputStreamWriter(
 410                                                         new FileOutputStream( filename) ) );
 411                         try{
 412                                 out.write("Kerneldimension");
 413                                 out.newLine();
 414                                 out.write(Integer.toString(dimension));
 415                                 out.newLine();
 416                                 out.newLine();
 417                                 out.write("number of rows");
 418                                 out.newLine();
 419                                 out.write(Integer.toString(length));
 420                                 out.newLine();
 421                                 out.newLine();
 422                                 out.write("Coefficients of the transform matrix:");
 423                                 out.newLine();
 424                                 for (int i=0; i < length; i++){
 425                                         String s = Double.toString(beta[i][0]);
 426                                         s += "    ";
 427                                         s += Double.toString(beta[i][1]);
 428                                         out.write(s);
 429                                         out.newLine();
 430                                 }
 431                                 out.newLine();
 432                                 out.write("normMean:");
 433                                 out.newLine();
 434                                 for (int i=0; i < length; i++){
 435                                         out.write(Double.toString(normMean[i]));
 436                                         out.newLine();
 437                                 }
 438                                 out.newLine();
 439                                 out.write("normVar: ");
 440                                 out.newLine();
 441                                 for (int i=0; i < length; i++){
 442                                         out.write(Double.toString(normVar[i]));
 443                                         out.newLine();
 444                                 }
 445                                 out.newLine();
 446                                 out.write("image size: ");
 447                                 out.newLine();
 448                                 out.write(width + "    " + height);
 449                                 out.close();
 450                         }
 451                         catch(IOException e){System.out.println("IOException");}
 452                 }
 453                 catch(FileNotFoundException e){System.out.println("File not found!");}
 454         }
 455
 456         public void load(String filename){
 457                 try{
 458                         BufferedReader in = new BufferedReader(new FileReader(filename));
 459                         try{
 460                                 String line = in.readLine(); //comment;
 461                                 dimension = Integer.parseInt(in.readLine());
 462                                 line = in.readLine(); //comment;
 463                                 line = in.readLine(); //comment;
 464                                 length = Integer.parseInt(in.readLine());
 465                                 line = in.readLine(); //comment;
 466                                 line = in.readLine(); //comment;
 467
 468                                 beta = new double[length][2];
 469
 470                                 for (int i=0; i < length; i++){
 471                                         line = in.readLine();
 472                                         int ind = line.indexOf(" ");
 473                                         beta[i][0] = Double.parseDouble(line.substring(0, ind));
 474                                         beta[i][1] = Double.parseDouble(line.substring(ind+4));
 475                                 }
 476
 477                                 line = in.readLine(); //comment;
 478                                 line = in.readLine(); //comment;
 479
 480                                 normMean = new double[length];
 481
 482                                 for (int i=0; i < length; i++){
 483                                         normMean[i]=Double.parseDouble(in.readLine());
 484                                 }
 485
 486                                 line = in.readLine(); //comment;
 487                                 line = in.readLine(); //comment;
 488
 489                                 normVar = new double[length];
 490
 491                                 for (int i=0; i < length; i++){
 492                                         normVar[i]=Double.parseDouble(in.readLine());
 493                                 }
 494                                 line = in.readLine(); //comment;
 495                                 line = in.readLine(); //comment;
 496                                 line = in.readLine();
 497                                 int ind = line.indexOf(" ");
 498                                 width = Integer.parseInt(line.substring(0, ind));
 499                                 height = Integer.parseInt(line.substring(ind+4));
 500                                 in.close();
 501
 502                                 print();
 503                         }
 504                         catch(IOException e){System.out.println("IOException");}
 505                 }
 506                 catch(FileNotFoundException e){System.out.println("File not found!");}
 507         }
 508
 509         public ImageProcessor[] transform(ImageProcessor ip){
 510                 if (!precalculated)
 511                         this.precalculateTransfom();
 512
 513                 ImageProcessor newIp = ip.createProcessor(ip.getWidth(), ip.getHeight());
 514                 if (ip instanceof ColorProcessor) ip.max(0);
 515                 ImageProcessor maskIp = new ByteProcessor(ip.getWidth(),ip.getHeight());
 516
 517                 for (int x=0; x < width; x++){
 518                         for (int y=0; y < height; y++){
 519                                 if (transField[x][y][0] == -1){
 520                                         continue;
 521                                 }
 522                                 newIp.set(x, y, (int) ip.getInterpolatedPixel((int)transField[x][y][0],(int)transField[x][y][1]));
 523                                 maskIp.set(x,y,255);
 524                         }
 525                 }
 526                 return new ImageProcessor[]{newIp, maskIp};
 527         }
 528
 529         private double[] multiply(double beta[][], double featureVector[]){
 530                 double[] result = {0.0,0.0};
 531
 532                 if (beta.length != featureVector.length){
 533                         IJ.log("Dimension of TransformMatrix and featureVector do not match!");
 534                         return new double[2];
 535                 }
 536
 537                 for (int i=0; i<featureVector.length; i++){
 538                         result[0] = result[0] + featureVector[i] * beta[i][0];
 539                         result[1] = result[1] + featureVector[i] * beta[i][1];
 540                 }
 541
 542                 return result;
 543         }
 544
 545         public double[] kernelExpand(double position[]){
 546                 double expanded[] = new double[length];
 547
 548                 int counter = 0;
 549                 for (int i=1; i<=dimension; i++){
 550                         for (double j=i; j>=0; j--){
 551                                 double val = Math.pow(position[0],j) * Math.pow(position[1],i-j);
 552                                 expanded[counter] = val;
 553                                 ++counter;
 554                         }
 555                 }
 556
 557                 for (int i=0; i<length-1; i++){
 558                         expanded[i] = expanded[i] - normMean[i];
 559                         expanded[i] = expanded[i] / normVar[i];
 560                 }
 561
 562                 expanded[length-1] = 100;
 563
 564                 return expanded;
 565         }
 566
 567
 568         public double[][] kernelExpandMatrixNormalize(double positions[][]){
 569                 normMean = new double[length];
 570                 normVar = new double[length];
 571
 572                 for (int i=0; i < length; i++){
 573                         normMean[i] = 0;
 574                         normVar[i] = 1;
 575                 }
 576
 577                 double expanded[][] = new double[positions.length][length];
 578
 579                 for (int i=0; i < positions.length; i++){
 580                         expanded[i] = kernelExpand(positions[i]);
 581                 }
 582
 583                 for (int i=0; i < length; i++){
 584                         double mean = 0;
 585                         double var = 0;
 586                         for (int j=0; j < expanded.length; j++){
 587                                 mean += expanded[j][i];
 588                         }
 589
 590                         mean /= expanded.length;
 591
 592                         for (int j=0; j < expanded.length; j++){
 593                                 var += (expanded[j][i] - mean)*(expanded[j][i] - mean);
 594                         }
 595                         var /= (expanded.length -1);
 596                         var = Math.sqrt(var);
 597
 598                         normMean[i] = mean;
 599                         normVar[i] = var;
 600                 }
 601
 602                 return kernelExpandMatrix(positions);
 603
 604         }
 605
 606         //this function uses the parameters already stored
 607         //in this object to normalize the positions given.
 608         public double[][] kernelExpandMatrix(double positions[][]){
 609
 610
 611                 double expanded[][] = new double[positions.length][length];
 612
 613                 for (int i=0; i < positions.length; i++){
 614                         expanded[i] = kernelExpand(positions[i]);
 615                 }
 616
 617                 return expanded;
 618
 619         }
 620
 621         public void inverseTransform(double range[][]){
 622                 Matrix expanded = new Matrix(kernelExpandMatrix(range));
 623                 Matrix b = new Matrix(beta);
 624
 625                 Matrix transformed = expanded.times(b);
 626                 expanded = new Matrix(kernelExpandMatrixNormalize(transformed.getArray()));
 627
 628                 Matrix r = new Matrix(range);
 629                 Matrix invBeta = expanded.transpose().times(expanded).inverse().times(expanded.transpose()).times(r);
 630                 setBeta(invBeta.getArray());
 631         }
 632
 633         //FIXME this takes way too much memory
 634         public void visualize(){
 635
 636                 int density = Math.max(width,height)/32;
 637                 int border = Math.max(width,height)/8;
 638
 639                 double[][] orig = new double[width *  height][2];
 640                 double[][] trans = new double[height * width][2];
 641                 double[][] gridOrigVert = new double[width*height][2];
 642                 double[][] gridTransVert = new double[width*height][2];
 643                 double[][] gridOrigHor = new double[width*height][2];
 644                 double[][] gridTransHor = new double[width*height][2];
 645
 646                 FloatProcessor magnitude = new FloatProcessor(width, height);
 647                 FloatProcessor angle = new FloatProcessor(width, height);
 648                 ColorProcessor quiver = new ColorProcessor(width, height);
 649                 ByteProcessor empty = new ByteProcessor(width+2*border, height+2*border);
 650                 quiver.setLineWidth(1);
 651                 quiver.setColor(Color.green);
 652
 653                 GeneralPath quiverField = new GeneralPath();
 654
 655                 float minM = 1000, maxM = 0;
 656                 float minArc = 5, maxArc = -6;
 657                 int countVert = 0, countHor = 0, countHorWhole = 0;
 658
 659                 for (int i=0; i < width; i++){
 660                         countHor = 0;
 661                         for (int j=0; j < height; j++){
 662                                 double[] position = {(double) i,(double) j};
 663                                 double[] posExpanded = kernelExpand(position);
 664                                 double[] newPosition = multiply(beta, posExpanded);
 665
 666                                 orig[i*j][0] = position[0];
 667                                 orig[i*j][1] = position[1];
 668
 669                                 trans[i*j][0] = newPosition[0];
 670                                 trans[i*j][1] = newPosition[1];
 671
 672                                 double m = (position[0] - newPosition[0]) * (position[0] - newPosition[0]);
 673                                 m += (position[1] - newPosition[1]) * (position[1] - newPosition[1]);
 674                                 m = Math.sqrt(m);
 675                                 magnitude.setf(i,j, (float) m);
 676                                 minM = Math.min(minM, (float) m);
 677                                 maxM = Math.max(maxM, (float) m);
 678
 679                                 double a = Math.atan2(position[0] - newPosition[0], position[1] - newPosition[1]);
 680                                 minArc = Math.min(minArc, (float) a);
 681                                 maxArc = Math.max(maxArc, (float) a);
 682                                 angle.setf(i,j, (float) a);
 683
 684                                 if (i%density == 0 && j%density == 0)
 685                                         drawQuiverField(quiverField, position[0], position[1], newPosition[0], newPosition[1]);
 686                                 if (i%density == 0){
 687                                         gridOrigVert[countVert][0] = position[0] + border;
 688                                         gridOrigVert[countVert][1] = position[1] + border;
 689                                         gridTransVert[countVert][0] = newPosition[0] + border;
 690                                         gridTransVert[countVert][1] = newPosition[1] + border;
 691                                         countVert++;
 692                                 }
 693                                 if (j%density == 0){
 694                                         gridOrigHor[countHor*width+i][0] = position[0] + border;
 695                                         gridOrigHor[countHor*width+i][1] = position[1] + border;
 696                                         gridTransHor[countHor*width+i][0] = newPosition[0] + border;
 697                                         gridTransHor[countHor*width+i][1] = newPosition[1] + border;
 698                                         countHor++;
 699                                         countHorWhole++;
 700                                 }
 701                         }
 702                 }
 703
 704                 magnitude.setMinAndMax(minM, maxM);
 705                 angle.setMinAndMax(minArc, maxArc);
 706                 //System.out.println(" " + minArc + " " + maxArc);
 707
 708                 ImagePlus magImg = new ImagePlus("Magnitude of Distortion Field", magnitude);
 709                 magImg.show();
 710
 711                 //              ImagePlus angleImg = new ImagePlus("Angle of Distortion Field Vectors", angle);
 712                 //              angleImg.show();
 713
 714                 ImagePlus quiverImg = new ImagePlus("Quiver Plot of Distortion Field", magnitude);
 715                 quiverImg.show();
 716                 quiverImg.getCanvas().setDisplayList(quiverField, Color.green, null );
 717                 quiverImg.updateAndDraw();
 718
 719                 //              GeneralPath gridOrig = new GeneralPath();
 720                 //              drawGrid(gridOrig, gridOrigVert, countVert, height);
 721                 //              drawGrid(gridOrig, gridOrigHor, countHorWhole, width);
 722                 //              ImagePlus gridImgOrig = new ImagePlus("Distortion Grid", empty);
 723                 //              gridImgOrig.show();
 724                 //              gridImgOrig.getCanvas().setDisplayList(gridOrig, Color.green, null );
 725                 //              gridImgOrig.updateAndDraw();
 726
 727                 GeneralPath gridTrans = new GeneralPath();
 728                 drawGrid(gridTrans, gridTransVert, countVert, height);
 729                 drawGrid(gridTrans, gridTransHor, countHorWhole, width);
 730                 ImagePlus gridImgTrans = new ImagePlus("Distortion Grid", empty);
 731                 gridImgTrans.show();
 732                 gridImgTrans.getCanvas().setDisplayList(gridTrans, Color.green, null );
 733                 gridImgTrans.updateAndDraw();
 734
 735                 //new FileSaver(quiverImg.getCanvas().imp).saveAsTiff("QuiverCanvas.tif");
 736                 new FileSaver(quiverImg).saveAsTiff("QuiverImPs.tif");
 737
 738                 System.out.println("FINISHED");
 739         }
 740
 741
 742         public void visualizeSmall(double lambda){
 743                 int density = Math.max(width,height)/32;
 744
 745                 double[][] orig = new double[2][width *  height];
 746                 double[][] trans = new double[2][height * width];
 747
 748                 FloatProcessor magnitude = new FloatProcessor(width, height);
 749
 750                 GeneralPath quiverField = new GeneralPath();
 751
 752                 float minM = 1000, maxM = 0;
 753                 float minArc = 5, maxArc = -6;
 754                 int countVert = 0, countHor = 0, countHorWhole = 0;
 755
 756                 for (int i=0; i < width; i++){
 757                         countHor = 0;
 758                         for (int j=0; j < height; j++){
 759                                 double[] position = {(double) i,(double) j};
 760                                 double[] posExpanded = kernelExpand(position);
 761                                 double[] newPosition = multiply(beta, posExpanded);
 762
 763                                 orig[0][i*j] = position[0];
 764                                 orig[1][i*j] = position[1];
 765
 766                                 trans[0][i*j] = newPosition[0];
 767                                 trans[1][i*j] = newPosition[1];
 768
 769                                 double m = (position[0] - newPosition[0]) * (position[0] - newPosition[0]);
 770                                 m += (position[1] - newPosition[1]) * (position[1] - newPosition[1]);
 771                                 m = Math.sqrt(m);
 772                                 magnitude.setf(i,j, (float) m);
 773                                 minM = Math.min(minM, (float) m);
 774                                 maxM = Math.max(maxM, (float) m);
 775
 776                                 if (i%density == 0 && j%density == 0)
 777                                         drawQuiverField(quiverField, position[0], position[1], newPosition[0], newPosition[1]);
 778                         }
 779                 }
 780
 781                 magnitude.setMinAndMax(minM, maxM);
 782                 ImagePlus quiverImg = new ImagePlus("Quiver Plot for lambda = "+lambda, magnitude);
 783                 quiverImg.show();
 784                 quiverImg.getCanvas().setDisplayList(quiverField, Color.green, null );
 785                 quiverImg.updateAndDraw();
 786
 787                 System.out.println("FINISHED");
 788         }
 789
 790
 791         public static void drawGrid(GeneralPath g, double[][] points, int count, int s){
 792                 for (int i=0; i < count - 1; i++){
 793                         if ((i+1)%s != 0){
 794                                 g.moveTo((float)points[i][0], (float)points[i][1]);
 795                                 g.lineTo((float)points[i+1][0], (float)points[i+1][1]);
 796                         }
 797                 }
 798         }
 799
 800         public static void drawQuiverField(GeneralPath qf, double x1, double y1, double x2, double y2)
 801         {
 802                 qf.moveTo((float)x1, (float)y1);
 803                 qf.lineTo((float)x2, (float)y2);
 804         }
 805
 806         public int getWidth(){
 807                 return width;
 808         }
 809
 810         public int getHeight(){
 811                 return height;
 812         }
 813
 814         @Override
 815         /**
 816          * TODO Make this more efficient
 817          */
 818         final public NonLinearTransform copy()
 819         {
 820                 final NonLinearTransform t = new NonLinearTransform();
 821                 t.init( toDataString() );
 822                 return t;
 823         }
 824
 825         public void set( final NonLinearTransform nlt )
 826         {
 827                 this.dimension = nlt.dimension;
 828                 this.height = nlt.height;
 829                 this.length = nlt.length;
 830                 this.precalculated = nlt.precalculated;
 831                 this.width = nlt.width;
 832
 833                 /* arrays by deep cloning */
 834                 this.beta = new double[ nlt.beta.length ][];
 835                 for ( int i = 0; i < nlt.beta.length; ++i )
 836                         this.beta[ i ] = nlt.beta[ i ].clone();
 837
 838                 this.normMean = nlt.normMean.clone();
 839                 this.normVar = nlt.normVar.clone();
 840                 this.transField = new double[ nlt.transField.length ][][];
 841
 842                 for ( int a = 0; a < nlt.transField.length; ++a )
 843                 {
 844                         this.transField[ a ] = new double[ nlt.transField[ a ].length ][];
 845                         for ( int b = 0; b < nlt.transField[ a ].length; ++b )
 846                                 this.transField[ a ][ b ] = nlt.transField[ a ][ b ].clone();
 847                 }
 848         }
 849 }