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 import mpicbg.trakem2.transform.TranslationModel2D;
  55
  56
  57 public class NonLinearTransform implements mpicbg.trakem2.transform.CoordinateTransform{
  58
  59                 private double[][] beta = null;
  60                 private double[] normMean = null;
  61                 private double[] normVar = null;
  62                 private double[][][] transField = null;
  63                 private int dimension = 0;
  64                 private int length = 0;
  65                 private int width = 0;
  66                 private int height = 0;
  67                 private boolean precalculated = false;
  68
  69                 public NonLinearTransform(double[][] b, double[] nm, double[] nv, int d, int w, int h){
  70                                 beta = b;
  71                                 normMean = nm;
  72                                 normVar = nv;
  73                                 dimension = d;
  74                                 length = (dimension + 1)*(dimension + 2)/2;
  75                                 width = w;
  76                                 height = h;
  77                 }
  78
  79                 public NonLinearTransform(int d, int w, int h){
  80                                 dimension = d;
  81                                 length = (dimension + 1)*(dimension + 2)/2;
  82
  83                                 beta = new double[length][2];
  84                                 normMean = new double[length];
  85                                 normVar = new double[length];
  86
  87                                 for (int i=0; i < length; i++){
  88                                                 normMean[i] = 0;
  89                                                 normVar[i] = 1;
  90                                 }
  91
  92                                 width = w;
  93                                 height = h;
  94                 }
  95
  96                 public NonLinearTransform(){};
  97
  98                 public NonLinearTransform(String filename){
  99                                 this.load(filename);
 100                 }
 101
 102                 public NonLinearTransform(double[][] coeffMatrix, int w, int h){
 103                                 length = coeffMatrix.length;
 104                                 beta = new double[length][2];
 105                                 normMean = new double[length];
 106                                 normVar = new double[length];
 107                                 width = w;
 108                                 height = h;
 109                                 dimension = (int)(-1.5 + Math.sqrt(0.25 + 2*length));
 110
 111                                 for(int i=0; i<length; i++){
 112                                                 beta[i][0] = coeffMatrix[0][i];
 113                                                 beta[i][1] = coeffMatrix[1][i];
 114                                                 normMean[i] = coeffMatrix[2][i];
 115                                                 normVar[i] = coeffMatrix[3][i];
 116                                 }
 117                 }
 118
 119
 120                 //implements mpicbg.trakem2
 121                 public void init( String data ) throws NumberFormatException{
 122                                 final String[] fields = data.split( " " );
 123                                 int c = 0;
 124
 125                                 dimension = Integer.parseInt(fields[c]); c++;
 126                                 length = Integer.parseInt(fields[c]); c++;
 127
 128                                 beta = new double[length][2];
 129                                 normMean = new double[length];
 130                                 normVar = new double[length];
 131
 132                                 if ( fields.length == 4 + 4*length )
 133                                                 {
 134                                                                 for (int i=0; i < length; i++){
 135                                                                                 beta[i][0] = Double.parseDouble(fields[c]); c++;
 136                                                                                 beta[i][1] = Double.parseDouble(fields[c]); c++;
 137                                                                 }
 138
 139                                                                 System.out.println("c: " + c);
 140
 141                                                                 for (int i=0; i < length; i++){
 142                                                                                 normMean[i] = Double.parseDouble(fields[c]); c++;
 143                                                                 }
 144
 145                                                                 System.out.println("c: " + c);
 146
 147                                                                 for (int i=0; i < length; i++){
 148                                                                                 normVar[i] = Double.parseDouble(fields[c]); c++;
 149                                                                 }
 150
 151                                                                 width = Integer.parseInt(fields[c]); c++;
 152                                                                 height = Integer.parseInt(fields[c]); c++;
 153                                                                 System.out.println("c: " + c);
 154
 155                                                 }
 156                                 else throw new NumberFormatException( "Inappropriate parameters for " + this.getClass().getCanonicalName() );
 157                 }
 158
 159
 160
 161                 public String toXML(final String indent){
 162                                 return new StringBuffer(indent).append("<ict_transform class=\"").append(this.getClass().getCanonicalName()).append("\" data=\"").append(toDataString()).append("\"/>").toString();
 163                 }
 164
 165                 public String toDataString(){
 166                                 String data = "";
 167                                 data += Integer.toString(dimension) + " ";
 168                                 data += Integer.toString(length) + " ";
 169
 170                                 for (int i=0; i < length; i++){
 171                                                 data += Double.toString(beta[i][0]) + " ";
 172                                                 data += Double.toString(beta[i][1]) + " ";
 173                                 }
 174
 175                                 for (int i=0; i < length; i++){
 176                                                 data += Double.toString(normMean[i]) + " ";
 177                                 }
 178
 179                                 for (int i=0; i < length; i++){
 180                                                 data += Double.toString(normVar[i]) + " ";
 181                                 }
 182                                 data += Integer.toString(width) + " ";
 183                                 data += Integer.toString(height) + " ";
 184
 185                                 return data;
 186
 187                 }
 188
 189                 public float[] apply( float[] location ){
 190
 191                                 double[] position = {(double) location[0], (double) location[1]};
 192                                 double[] featureVector = kernelExpand(position);
 193                                 double[] newPosition = multiply(beta, featureVector);
 194
 195                                 float[] newLocation = new float[2];
 196                                 newLocation[0] = (float) newPosition[0];
 197                                 newLocation[1] = (float) newPosition[1];
 198
 199                                 return newLocation;
 200                 }
 201
 202                 public void applyInPlace( float[] location ){
 203                                 double[] position = {(double) location[0], (double) location[1]};
 204                                 double[] featureVector = kernelExpand(position);
 205                                 double[] newPosition = multiply(beta, featureVector);
 206
 207                                 location[0] = (float) newPosition[0];
 208                                 location[1] = (float) newPosition[1];
 209                 }
 210
 211                 void precalculateTransfom(){
 212                                 transField = new double[width][height][2];
 213                                 //double minX = width, minY = height, maxX = 0, maxY = 0;
 214
 215                                 for (int x=0; x<width; x++){
 216                                                 for (int y=0; y<height; y++){
 217                                                                 double[] position = {x,y};
 218                                                                 double[] featureVector = kernelExpand(position);
 219                                                                 double[] newPosition = multiply(beta, featureVector);
 220
 221                                                                 if ((newPosition[0] < 0) || (newPosition[0] >= width) ||
 222                                                                                 (newPosition[1] < 0) || (newPosition[1] >= height))
 223                                                                                 {
 224                                                                                                 transField[x][y][0] = -1;
 225                                                                                                 transField[x][y][1] = -1;
 226                                                                                                 continue;
 227                                                                                 }
 228
 229                                                                 transField[x][y][0] = newPosition[0];
 230                                                                 transField[x][y][1] = newPosition[1];
 231
 232                                                                 //minX = Math.min(minX, x);
 233                                                                 //minY = Math.min(minY, y);
 234                                                                 //maxX = Math.max(maxX, x);
 235                                                                 //maxY = Math.max(maxY, y);
 236
 237                                                 }
 238                                 }
 239
 240                                 precalculated = true;
 241                 }
 242
 243                 public double[][] getCoefficients(){
 244                                 double[][] coeffMatrix = new double[4][length];
 245
 246                                 for(int i=0; i<length; i++){
 247                                                 coeffMatrix[0][i] = beta[i][0];
 248                                                 coeffMatrix[1][i] = beta[i][1];
 249                                                 coeffMatrix[2][i] = normMean[i];
 250                                                 coeffMatrix[3][i] = normVar[i];
 251
 252                                 }
 253                                 return coeffMatrix;
 254                 }
 255
 256                 public void setBeta(double[][] b){
 257                                 beta = b;
 258                                 //FIXME: test if normMean and normVar are still valid for this beta
 259                 }
 260
 261                 public void print(){
 262                                 System.out.println("beta:");
 263                                 for (int i=0; i < beta.length; i++){
 264                                                 for (int j=0; j < beta[i].length; j++){
 265                                                                 System.out.print(beta[i][j]);
 266                                                                 System.out.print(" ");
 267                                                 }
 268                                                 System.out.println();
 269                                 }
 270
 271                                 System.out.println("normMean:");
 272                                 for (int i=0; i < normMean.length; i++){
 273                                                 System.out.print(normMean[i]);
 274                                                 System.out.print(" ");
 275                                 }
 276
 277                                 System.out.println("normVar:");
 278                                 for (int i=0; i < normVar.length; i++){
 279                                                 System.out.print(normVar[i]);
 280                                                 System.out.print(" ");
 281                                 }
 282
 283                                 System.out.println("Image size:");
 284                                 System.out.println("width: " + width + " height: " + height);
 285
 286                                 System.out.println();
 287
 288                 }
 289
 290                 public void save( final String filename )
 291                 {
 292                                 try{
 293                                                 BufferedWriter out = new BufferedWriter(
 294                                                         new OutputStreamWriter(
 295                                                         new FileOutputStream( filename) ) );
 296                                                 try{
 297                                                                 out.write("Kerneldimension");
 298                                                                 out.newLine();
 299                                                                 out.write(Integer.toString(dimension));
 300                                                                 out.newLine();
 301                                                                 out.newLine();
 302                                                                 out.write("number of rows");
 303                                                                 out.newLine();
 304                                                                 out.write(Integer.toString(length));
 305                                                                 out.newLine();
 306                                                                 out.newLine();
 307                                                                 out.write("Coefficients of the transform matrix:");
 308                                                                 out.newLine();
 309                                                                 for (int i=0; i < length; i++){
 310                                                                                 String s = Double.toString(beta[i][0]);
 311                                                                                 s += "    ";
 312                                                                                 s += Double.toString(beta[i][1]);
 313                                                                                 out.write(s);
 314                                                                                 out.newLine();
 315                                                                 }
 316                                                                 out.newLine();
 317                                                                 out.write("normMean:");
 318                                                                 out.newLine();
 319                                                                 for (int i=0; i < length; i++){
 320                                                                                 out.write(Double.toString(normMean[i]));
 321                                                                                 out.newLine();
 322                                                                 }
 323                                                                 out.newLine();
 324                                                                 out.write("normVar: ");
 325                                                                 out.newLine();
 326                                                                 for (int i=0; i < length; i++){
 327                                                                                 out.write(Double.toString(normVar[i]));
 328                                                                                 out.newLine();
 329                                                                 }
 330                                                                 out.newLine();
 331                                                                 out.write("image size: ");
 332                                                                 out.newLine();
 333                                                                 out.write(width + "    " + height);
 334                                                                 out.close();
 335                                                 }
 336                                                 catch(IOException e){System.out.println("IOException");}
 337                                 }
 338                                 catch(FileNotFoundException e){System.out.println("File not found!");}
 339                 }
 340
 341                 public void load(String filename){
 342                                 try{
 343                                                 BufferedReader in = new BufferedReader(new FileReader(filename));
 344                                                 try{
 345                                                                 String line = in.readLine(); //comment;
 346                                                                 dimension = Integer.parseInt(in.readLine());
 347                                                                 line = in.readLine(); //comment;
 348                                                                 line = in.readLine(); //comment;
 349                                                                 length = Integer.parseInt(in.readLine());
 350                                                                 line = in.readLine(); //comment;
 351                                                                 line = in.readLine(); //comment;
 352
 353                                                                 beta = new double[length][2];
 354
 355                                                                 for (int i=0; i < length; i++){
 356                                                                                 line = in.readLine();
 357                                                                                 int ind = line.indexOf(" ");
 358                                                                                 beta[i][0] = Double.parseDouble(line.substring(0, ind));
 359                                                                                 beta[i][1] = Double.parseDouble(line.substring(ind+4));
 360                                                                 }
 361
 362                                                                 line = in.readLine(); //comment;
 363                                                                 line = in.readLine(); //comment;
 364
 365                                                                 normMean = new double[length];
 366
 367                                                                 for (int i=0; i < length; i++){
 368                                                                                 normMean[i]=Double.parseDouble(in.readLine());
 369                                                                 }
 370
 371                                                                 line = in.readLine(); //comment;
 372                                                                 line = in.readLine(); //comment;
 373
 374                                                                 normVar = new double[length];
 375
 376                                                                 for (int i=0; i < length; i++){
 377                                                                                 normVar[i]=Double.parseDouble(in.readLine());
 378                                                                 }
 379                                                                 line = in.readLine(); //comment;
 380                                                                 line = in.readLine(); //comment;
 381                                                                 line = in.readLine();
 382                                                                 int ind = line.indexOf(" ");
 383                                                                 width = Integer.parseInt(line.substring(0, ind));
 384                                                                 height = Integer.parseInt(line.substring(ind+4));
 385                                                                 in.close();
 386
 387                                                                 print();
 388                                                 }
 389                                                 catch(IOException e){System.out.println("IOException");}
 390                                 }
 391                                 catch(FileNotFoundException e){System.out.println("File not found!");}
 392                 }
 393
 394                 public ImageProcessor[] transform(ImageProcessor ip){
 395                                 if (!precalculated)
 396                                                 this.precalculateTransfom();
 397
 398                                 ImageProcessor newIp = ip.createProcessor(ip.getWidth(), ip.getHeight());
 399                                 if (ip instanceof ColorProcessor) ip.max(0);
 400                                 ImageProcessor maskIp = new ByteProcessor(ip.getWidth(),ip.getHeight());
 401
 402                                 for (int x=0; x < width; x++){
 403                                                 for (int y=0; y < height; y++){
 404                                                                 if (transField[x][y][0] == -1){
 405                                                                                 continue;
 406                                                                 }
 407                                                                 newIp.set(x, y, (int) ip.getInterpolatedPixel((int)transField[x][y][0],(int)transField[x][y][1]));
 408                                                                 maskIp.set(x,y,255);
 409                                                 }
 410                                 }
 411                                 return new ImageProcessor[]{newIp, maskIp};
 412                 }
 413
 414                 private double[] multiply(double beta[][], double featureVector[]){
 415                                 double[] result = {0.0,0.0};
 416
 417                                 if (beta.length != featureVector.length){
 418                                                 IJ.showMessage("Dimension of TransformMatrix and featureVector do not match!");
 419                                                 return new double[2];
 420                                 }
 421
 422                                 for (int i=0; i<featureVector.length; i++){
 423                                                 result[0] = result[0] + featureVector[i] * beta[i][0];
 424                                                 result[1] = result[1] + featureVector[i] * beta[i][1];
 425                                 }
 426
 427                                 return result;
 428                 }
 429
 430                 public double[] kernelExpand(double position[]){
 431                                 double expanded[] = new double[length];
 432
 433                                 int counter = 0;
 434                                 for (int i=1; i<=dimension; i++){
 435                                                 for (double j=i; j>=0; j--){
 436                                                                 double val = Math.pow(position[0],j) * Math.pow(position[1],i-j);
 437                                                                 expanded[counter] = val;
 438                                                                 ++counter;
 439                                                 }
 440                                 }
 441
 442                                 for (int i=0; i<length-1; i++){
 443                                                 expanded[i] = expanded[i] - normMean[i];
 444                                                 expanded[i] = expanded[i] / normVar[i];
 445                                 }
 446
 447                                 expanded[length-1] = 100;
 448
 449                                 return expanded;
 450                 }
 451
 452
 453                 public double[][] kernelExpandMatrixNormalize(double positions[][]){
 454                                 normMean = new double[length];
 455                                 normVar = new double[length];
 456
 457                                 for (int i=0; i < length; i++){
 458                                                 normMean[i] = 0;
 459                                                 normVar[i] = 1;
 460                                 }
 461
 462                                 double expanded[][] = new double[positions.length][length];
 463
 464                                 for (int i=0; i < positions.length; i++){
 465                                                 expanded[i] = kernelExpand(positions[i]);
 466                                 }
 467
 468                                 for (int i=0; i < length; i++){
 469                                                 double mean = 0;
 470                                                 double var = 0;
 471                                                 for (int j=0; j < expanded.length; j++){
 472                                                                 mean += expanded[j][i];
 473                                                 }
 474
 475                                                 mean /= expanded.length;
 476
 477                                                 for (int j=0; j < expanded.length; j++){
 478                                                                 var += (expanded[j][i] - mean)*(expanded[j][i] - mean);
 479                                                 }
 480                                                 var /= (expanded.length -1);
 481                                                 var = Math.sqrt(var);
 482
 483                                                 normMean[i] = mean;
 484                                                 normVar[i] = var;
 485                                 }
 486
 487                                 return kernelExpandMatrix(positions);
 488
 489                 }
 490
 491                 //this function uses the parameters already stored
 492                 //in this object to normalize the positions given.
 493                 public double[][] kernelExpandMatrix(double positions[][]){
 494
 495
 496                                 double expanded[][] = new double[positions.length][length];
 497
 498                                 for (int i=0; i < positions.length; i++){
 499                                                 expanded[i] = kernelExpand(positions[i]);
 500                                 }
 501
 502                                 return expanded;
 503
 504                 }
 505
 506                 public void inverseTransform(double range[][]){
 507                                 Matrix expanded = new Matrix(kernelExpandMatrix(range));
 508                                 Matrix b = new Matrix(beta);
 509
 510                                 Matrix transformed = expanded.times(b);
 511                                 expanded = new Matrix(kernelExpandMatrixNormalize(transformed.getArray()));
 512
 513                                 Matrix r = new Matrix(range);
 514                                 Matrix invBeta = expanded.transpose().times(expanded).inverse().times(expanded.transpose()).times(r);
 515                                 setBeta(invBeta.getArray());
 516                 }
 517
 518     //FIXME this takes way too much memory
 519                 public void visualize(){
 520
 521                                 int density = Math.max(width,height)/32;
 522                                 int border = Math.max(width,height)/8;
 523
 524                                 double[][] orig = new double[width *  height][2];
 525                                 double[][] trans = new double[height * width][2];
 526                                 double[][] gridOrigVert = new double[width*height][2];
 527                                 double[][] gridTransVert = new double[width*height][2];
 528                                 double[][] gridOrigHor = new double[width*height][2];
 529                                 double[][] gridTransHor = new double[width*height][2];
 530
 531                                 FloatProcessor magnitude = new FloatProcessor(width, height);
 532                                 FloatProcessor angle = new FloatProcessor(width, height);
 533                                 ColorProcessor quiver = new ColorProcessor(width, height);
 534                                 ByteProcessor empty = new ByteProcessor(width+2*border, height+2*border);
 535                                 quiver.setLineWidth(1);
 536                                 quiver.setColor(Color.green);
 537
 538                                 GeneralPath quiverField = new GeneralPath();
 539
 540                                 float minM = 1000, maxM = 0;
 541                                 float minArc = 5, maxArc = -6;
 542                                 int countVert = 0, countHor = 0, countHorWhole = 0;
 543
 544                                 for (int i=0; i < width; i++){
 545                                                 countHor = 0;
 546                                                 for (int j=0; j < height; j++){
 547                                                                 double[] position = {(double) i,(double) j};
 548                                                                 double[] posExpanded = kernelExpand(position);
 549                                                                 double[] newPosition = multiply(beta, posExpanded);
 550
 551                                                                 orig[i*j][0] = position[0];
 552                                                                 orig[i*j][1] = position[1];
 553
 554                                                                 trans[i*j][0] = newPosition[0];
 555                                                                 trans[i*j][1] = newPosition[1];
 556
 557                                                                 double m = (position[0] - newPosition[0]) * (position[0] - newPosition[0]);
 558                                                                 m += (position[1] - newPosition[1]) * (position[1] - newPosition[1]);
 559                                                                 m = Math.sqrt(m);
 560                                                                 magnitude.setf(i,j, (float) m);
 561                                                                 minM = Math.min(minM, (float) m);
 562                                                                 maxM = Math.max(maxM, (float) m);
 563
 564                                                                 double a = Math.atan2(position[0] - newPosition[0], position[1] - newPosition[1]);
 565                                                                 minArc = Math.min(minArc, (float) a);
 566                                                                 maxArc = Math.max(maxArc, (float) a);
 567                                                                 angle.setf(i,j, (float) a);
 568
 569                                                                 if (i%density == 0 && j%density == 0)
 570                                                                                 drawQuiverField(quiverField, position[0], position[1], newPosition[0], newPosition[1]);
 571                                                                 if (i%density == 0){
 572                                                                                 gridOrigVert[countVert][0] = position[0] + border;
 573                                                                                 gridOrigVert[countVert][1] = position[1] + border;
 574                                                                                 gridTransVert[countVert][0] = newPosition[0] + border;
 575                                                                                 gridTransVert[countVert][1] = newPosition[1] + border;
 576                                                                                 countVert++;
 577                                                                 }
 578                                                                 if (j%density == 0){
 579                                                                                 gridOrigHor[countHor*width+i][0] = position[0] + border;
 580                                                                                 gridOrigHor[countHor*width+i][1] = position[1] + border;
 581                                                                                 gridTransHor[countHor*width+i][0] = newPosition[0] + border;
 582                                                                                 gridTransHor[countHor*width+i][1] = newPosition[1] + border;
 583                                                                                 countHor++;
 584                                                                                 countHorWhole++;
 585                                                                 }
 586                                                 }
 587                                 }
 588
 589                                 magnitude.setMinAndMax(minM, maxM);
 590                                 angle.setMinAndMax(minArc, maxArc);
 591                                 //System.out.println(" " + minArc + " " + maxArc);
 592
 593                                 ImagePlus magImg = new ImagePlus("Magnitude of Distortion Field", magnitude);
 594                                 magImg.show();
 595
 596                                 //              ImagePlus angleImg = new ImagePlus("Angle of Distortion Field Vectors", angle);
 597                                 //              angleImg.show();
 598
 599                                 ImagePlus quiverImg = new ImagePlus("Quiver Plot of Distortion Field", magnitude);
 600                                 quiverImg.show();
 601                                 quiverImg.getCanvas().setDisplayList(quiverField, Color.green, null );
 602                                 quiverImg.updateAndDraw();
 603
 604                                 //              GeneralPath gridOrig = new GeneralPath();
 605                                 //              drawGrid(gridOrig, gridOrigVert, countVert, height);
 606                                 //              drawGrid(gridOrig, gridOrigHor, countHorWhole, width);
 607                                 //              ImagePlus gridImgOrig = new ImagePlus("Distortion Grid", empty);
 608                                 //              gridImgOrig.show();
 609                                 //              gridImgOrig.getCanvas().setDisplayList(gridOrig, Color.green, null );
 610                                 //              gridImgOrig.updateAndDraw();
 611
 612                                 GeneralPath gridTrans = new GeneralPath();
 613                                 drawGrid(gridTrans, gridTransVert, countVert, height);
 614                                 drawGrid(gridTrans, gridTransHor, countHorWhole, width);
 615                                 ImagePlus gridImgTrans = new ImagePlus("Distortion Grid", empty);
 616                                 gridImgTrans.show();
 617                                 gridImgTrans.getCanvas().setDisplayList(gridTrans, Color.green, null );
 618                                 gridImgTrans.updateAndDraw();
 619
 620                                 //new FileSaver(quiverImg.getCanvas().imp).saveAsTiff("QuiverCanvas.tif");
 621                                 new FileSaver(quiverImg).saveAsTiff("QuiverImPs.tif");
 622
 623                                 System.out.println("FINISHED");
 624                 }
 625
 626
 627                 public void visualizeSmall(double lambda){
 628                                 int density = Math.max(width,height)/32;
 629
 630                                 double[][] orig = new double[width *  height][2];
 631                                 double[][] trans = new double[height * width][2];
 632
 633                                 FloatProcessor magnitude = new FloatProcessor(width, height);
 634                                 ColorProcessor quiver = new ColorProcessor(width, height);
 635                                 quiver.setLineWidth(1);
 636                                 quiver.setColor(Color.green);
 637
 638                                 GeneralPath quiverField = new GeneralPath();
 639
 640                                 float minM = 1000, maxM = 0;
 641                                 float minArc = 5, maxArc = -6;
 642                                 int countVert = 0, countHor = 0, countHorWhole = 0;
 643
 644                                 for (int i=0; i < width; i++){
 645                                                 countHor = 0;
 646                                                 for (int j=0; j < height; j++){
 647                                                                 double[] position = {(double) i,(double) j};
 648                                                                 double[] posExpanded = kernelExpand(position);
 649                                                                 double[] newPosition = multiply(beta, posExpanded);
 650
 651                                                                 orig[i*j][0] = position[0];
 652                                                                 orig[i*j][1] = position[1];
 653
 654                                                                 trans[i*j][0] = newPosition[0];
 655                                                                 trans[i*j][1] = newPosition[1];
 656
 657                                                                 double m = (position[0] - newPosition[0]) * (position[0] - newPosition[0]);
 658                                                                 m += (position[1] - newPosition[1]) * (position[1] - newPosition[1]);
 659                                                                 m = Math.sqrt(m);
 660                                                                 magnitude.setf(i,j, (float) m);
 661                                                                 minM = Math.min(minM, (float) m);
 662                                                                 maxM = Math.max(maxM, (float) m);
 663
 664                                                                 if (i%density == 0 && j%density == 0)
 665                                                                                 drawQuiverField(quiverField, position[0], position[1], newPosition[0], newPosition[1]);
 666                                                 }
 667                                 }
 668
 669                                 magnitude.setMinAndMax(minM, maxM);
 670                                 ImagePlus quiverImg = new ImagePlus("Quiver Plot for lambda = "+lambda, magnitude);
 671                                 quiverImg.show();
 672                                 quiverImg.getCanvas().setDisplayList(quiverField, Color.green, null );
 673                                 quiverImg.updateAndDraw();
 674
 675                                 System.out.println("FINISHED");
 676                 }
 677
 678
 679                 public static void drawGrid(GeneralPath g, double[][] points, int count, int s){
 680                                 for (int i=0; i < count - 1; i++){
 681                                                 if ((i+1)%s != 0){
 682                                                                 g.moveTo((float)points[i][0], (float)points[i][1]);
 683                                                                 g.lineTo((float)points[i+1][0], (float)points[i+1][1]);
 684                                                 }
 685                                 }
 686                 }
 687
 688                 public static void drawQuiverField(GeneralPath qf, double x1, double y1, double x2, double y2)
 689                 {
 690                                 qf.moveTo((float)x1, (float)y1);
 691                                 qf.lineTo((float)x2, (float)y2);
 692                 }
 693
 694                 public int getWidth(){
 695                                 return width;
 696                 }
 697
 698                 public int getHeight(){
 699                                 return height;
 700                 }
 701
 702                 @Override
 703                 /**
 704                  * TODO Make this more efficient
 705                  */
 706                 final public NonLinearTransform clone()
 707                 {
 708                         final NonLinearTransform t = new NonLinearTransform();
 709                         t.init( toDataString() );
 710                         return t;
 711                 }
 712 }