?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Make CachingThread aware of potential concurrent access,
[trakem2/tony-azevedo.git] / lenscorrection / NonLinearTransform.java
bloba4af95d857f8a0d20fb6a29c3c0fa3fb4bdcb67f
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 }