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