1 /* Panorama_Tools - Generate, Edit and Convert Panoramic Images
2 Copyright (C) 1998,1999 - Helmut Dersch der@fh-furtwangen.de
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
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 software; see the file COPYING. If not, a copy
16 can be downloaded from http://www.gnu.org/licenses/gpl.html, or
17 obtained by writing to the Free Software Foundation, Inc.,
18 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */
19
20 /*------------------------------------------------------------*/
21
22 /* ---- Revision history ----
23
24 May 2004, Rik Littlefield, reworked fcnPano and related functions as follows:
25 1) For normal control points, allows exposing latitude and longitude
26 error components separately to the optimizer. This trades
27 faster convergence for slight loss of stability. This
28 behavior can be changed at runtime by calling setFcnPanoNperCP().
29 (new capability)
30 2) Optimize distance^2 instead of distance^4 for hor, vert, and
31 line control points (bug fix)
32 3) Scale errors by change in average fov. This stabilizes
33 fov optimization and allows its use in more cases with partial panos.
34 (new capability)
35 4) Improve accuracy of angular distance calculation by using asin
36 instead of acos (results improvement)
37 5) Consistently report errors in units of pixels scaled to current
38 panorama size (feature change)
39 6) Report rms error during optimization (bug fix)
40 */
41 #include <math.h>
42 #include "filter.h"
43 #include "f2c.h"
44 #include <float.h>
45 #include <assert.h>
46 #include "PTcommon.h"
47
48 #define C_FACTOR 100.0
49
50 static AlignInfo *optInfo; // This struct holds all informations for the optimization
51
52 static double initialAvgFov; // these three for fov stabilization
53 static double avgfovFromSAP;
54 static int needInitialAvgFov;
55
56 #define ADJUST_LOG_FILENAME "PToolsLog.txt" // file name for logging, if enabled
57 #define ADJUST_LOGGING_ENABLED 0
58
59 FILE* adjustLogFile = 0;
60
61
62 void ColCorrect( Image *im, double ColCoeff[3][2] );
63 void GetColCoeff( Image *src, Image *buf, double ColCoeff[3][2] );
64 void getControlPoints( Image *im, struct controlPoint *cp );
65 void writeControlPoints( struct controlPoint *cp,char* cdesc );
66 int CheckParams( AlignInfo *g );
67 static int CheckMakeParams( aPrefs *aP);
68 //static int GetOverlapRect( PTRect *OvRect, PTRect *r1, PTRect *r2 );
69 int AddEdgePoints( AlignInfo *gl );
70 int pt_average( uint8_t* pixel, int BytesPerLine, double rgb[3], int bytesPerChannel );
71 double distsqLine(int N0, int N1);
72
73
74
panoAdjustPrintMakeParams(char * msg,struct MakeParams * mp,Image * im)75 void panoAdjustPrintMakeParams(char *msg, struct MakeParams *mp, Image *im)
76 {
77 printf("-------------%s\n", msg);
78 if (mp != NULL) {
79 printf("distance %f\n", mp->distance);
80 printf("shear[0] %f\n", mp->shear[0]);
81 printf("shear[1] %f\n", mp->shear[1]);
82 printf("rot[0] %f\n", mp->rot[0]);
83 printf("rot[1] %f\n", mp->rot[1]);
84 printf("tilt[0] %f\n", mp->tilt[0]);
85 printf("tilt[1] %f\n", mp->tilt[1]);
86 printf("tilt[2] %f\n", mp->tilt[2]);
87 printf("tilt[3] %f\n", mp->tilt[3]);
88
89 printf("trans[0] %f\n", mp->trans[0]);
90 printf("trans[1] %f\n", mp->trans[1]);
91 printf("trans[2] %f\n", mp->trans[2]);
92 printf("trans[3] %f\n", mp->trans[3]);
93 printf("trans[4] %f\n", mp->trans[4]);
94
95 printf("test[0] %f\n", mp->test[0]);
96 printf("test[1] %f\n", mp->test[1]);
97 printf("test[2] %f\n", mp->test[2]);
98 printf("test[3] %f\n", mp->test[3]);
99
100
101 printf("mp->horizontal %f\n", mp->horizontal);
102 printf("mp->vertical %f\n", mp->vertical);
103 }
104 panoPrintImage(msg, im);
105 printf("\n\n");
106
107 }
108
109
adjust(TrformStr * TrPtr,aPrefs * prefs)110 void adjust(TrformStr *TrPtr, aPrefs *prefs)
111 {
112 int destwidth, destheight;
113 aPrefs aP, *aPtr=NULL;
114 #if 0
115 int nt = 0; // Morph parameters
116 PTTriangle *ts=NULL;
117 PTTriangle *td=NULL;
118 #endif
119 SetAdjustDefaults(&aP);
120
121 switch( prefs->mode & 7 )// Should we use prefs, or read from script?
122 {
123 case _insert:
124 case _extract:
125 if( prefs->mode & _useScript ){
126 aPtr = readAdjustLine( &(prefs->scriptFile) );
127 if(aPtr==NULL){
128 PrintError("Error processing script file" );
129 TrPtr->success = 0;
130 return;
131 }
132 memcpy(&aP, aPtr, sizeof(aPrefs));
133 free(aPtr); aPtr = &aP;
134
135 if( (TrPtr->mode & 7) == _usedata ){ // Report panorama format and stitching info back to calling app.
136 memcpy( &prefs->pano, &aP.pano, sizeof( Image ) );
137 memcpy( &prefs->sBuf, &aP.sBuf, sizeof( stBuf ) );
138 }
139
140 TrPtr->interpolator = aP.interpolator;
141 TrPtr->gamma = aP.gamma;
142 TrPtr->fastStep = aP.fastStep;
143
144 #if 0
145 int readmode = 1;
146 aPtr = &aP;
147 gsPrPtr->interpolator = TrPtr->interpolator;
148 gsPrPtr->gamma = TrPtr->gamma;
149 gsPrPtr->fastStep = TrPtr->fastStep;
150 if( TrPtr->mode & _destSupplied ){
151 PTRect* p = &TrPtr->dest->selection;
152 if( !(p->bottom == 0 && p->right == 0) &&
153 !(p->right == TrPtr->dest->width &&
154 p->bottom == TrPtr->dest->height) )
155 readmode = 0;
156 }
157 if( readAdjust( aPtr, &(prefs->scriptFile), readmode, gsPrPtr ) != 0 )
158 {
159 PrintError("Error processing script file" );
160 TrPtr->success = 0;
161 return;
162 }
163 if( (TrPtr->mode & 7) == _usedata ) // Report panorama format and stitching info back to calling app.
164 {
165 memcpy( &prefs->pano, &aP.pano, sizeof( Image ) );
166 memcpy( &prefs->sBuf, &aP.sBuf, sizeof( stBuf ) );
167 }
168 // Use modevalues read from script
169 TrPtr->interpolator = gsPrPtr->interpolator;
170 TrPtr->gamma = gsPrPtr->gamma;
171 TrPtr->fastStep = gsPrPtr->fastStep;
172
173 // Parse script again, now reading triangles if morphing requested
174 if( aPtr->im.cP.correction_mode & correction_mode_morph )
175 {
176 char* script;
177 AlignInfo ainf;
178 int nIm, nPts; // Number of image being processed
179 Image im[2];
180
181 script = LoadScript( &(prefs->scriptFile) );
182 if( script != NULL ) // We can read the scriptfile
183 {
184 nIm = numLines( script, '!' ) - 1;
185
186 if( nIm < 0)
187 nIm = numLines( script, 'o' ) - 1;
188
189 // Set ainf
190 ainf.nt = 0;
191 ainf.t = NULL;
192 ainf.numIm = 2;
193 ainf.im = im;
194 memcpy( &ainf.pano, &aP.pano, sizeof( Image ));
195 memcpy( &ainf.im[0], &aP.pano, sizeof( Image ));
196 memcpy( &ainf.im[1], &aP.pano, sizeof( Image ));
197
198 nPts = ReadMorphPoints( script, &ainf, nIm );
199 if(nPts > 0) // Found Points
200 {
201 AddEdgePoints( &ainf );
202 TriangulatePoints( &ainf, 1 );
203 nt = ainf.nt;
204 if(nt > 0)
205 {
206 SortControlPoints ( &ainf, 1 );
207 SetSourceTriangles ( &ainf, 1, &td );
208 SetDestTriangles ( &ainf, 1, &ts );
209 }
210 }
211 if(ainf.numPts > 0) free(ainf.cpt);
212 free( script );
213 }
214 }
215 #endif
216 }else{
217 aPtr = prefs;
218 }
219 break;
220 default:
221 break;
222 }
223 switch( prefs->mode & 7)
224 {
225 case _insert: // Create a panoramic image using src; merge with buffer if required
226 // Find brightest rectangle if this is a circular fisheye image
227 {
228 Image ImCrop, *theSrc=NULL;
229 // Initialise at least the data pointer since cutTheFrame may not do it
230 ImCrop.data = NULL;
231
232 if( aPtr->im.format ==_fisheye_circ && aPtr->im.cP.cutFrame )
233 {
234 int fwidth = TrPtr->src->width, fheight = TrPtr->src->height;
235
236 if( aPtr->im.cP.frame ) // subtract framewidth from width/height
237 {
238 fwidth = TrPtr->src->width - aPtr->im.cP.frame;
239 if( aPtr->im.cP.frame < fwidth ) fwidth -= aPtr->im.cP.frame;
240 if( aPtr->im.cP.frame < fheight) fheight-= aPtr->im.cP.frame;
241 }
242 else
243 {
244 if( aPtr->im.cP.fwidth > 0)
245 fwidth = aPtr->im.cP.fwidth;
246 if( aPtr->im.cP.fheight > 0)
247 fheight = aPtr->im.cP.fheight;
248 }
249
250 if( cutTheFrame( &ImCrop, TrPtr->src, fwidth, fheight, TrPtr->mode & _show_progress ) != 0 )
251 {
252 PrintError("Error Cropping Image");
253 TrPtr->success = 0;
254 return;
255 }
256 theSrc = TrPtr->src;
257 TrPtr->src = &ImCrop;
258
259 }
260 // Image params are set as src
261 aPtr->im.width = TrPtr->src->width;
262 aPtr->im.height = TrPtr->src->height;
263
264 // Pano is set to buffer, if merging requested; else as prefs
265 if( *aPtr->sBuf.srcName != 0 )
266 {
267 if (LoadBufImage( &(aPtr->pano), aPtr->sBuf.srcName, 0) != 0 )
268 {
269 PrintError( "Error loading Buffer; trying without" );
270 }
271 }
272
273 if( aPtr->pano.width == 0 && aPtr->im.hfov != 0.0)
274 {
275 aPtr->pano.width = (aPtr->im.width * aPtr->pano.hfov / aPtr->im.hfov);
276 aPtr->pano.width/=10; aPtr->pano.width*=10;
277 }
278 if( aPtr->pano.height == 0 )
279 aPtr->pano.height = aPtr->pano.width/2;
280
281 destheight = aPtr->pano.height;
282 destwidth = aPtr->pano.width;
283
284 if( destheight == 0 || destwidth == 0 )
285 {
286 PrintError("Please set Panorama width/height" );
287 TrPtr->success = 0;
288 goto _insert_exit;
289 }
290
291
292 if( SetDestImage( TrPtr, destwidth, destheight) != 0)
293 {
294 PrintError("Could not allocate %ld bytes",TrPtr->dest->dataSize );
295 TrPtr->success = 0;
296 goto _insert_exit;
297 }
298 TrPtr->mode |= _honor_valid;
299 CopyPosition( TrPtr->src, &(aPtr->im) );
300 CopyPosition( TrPtr->dest, &(aPtr->pano) );
301 // JMW 2008/01/07 Alpha is valid data don't override it with blank data
302 // addAlpha( TrPtr->src ); // Add alpha channel to indicate valid data
303
304 aPtr->mode = prefs->mode; // For checkparam
305 MakePano( TrPtr, aPtr );
306
307 if(aPtr->ts) free(aPtr->ts);
308 if(aPtr->td) free(aPtr->td);
309
310 // Stitch images; Proceed only if panoramic image valid
311
312 if( TrPtr->success )
313 {
314 if( *(aPtr->sBuf.srcName) != 0 ){ // We have to merge in one images
315 // Load the bufferimage
316 if( LoadBufImage( &aPtr->pano, aPtr->sBuf.srcName, 1 ) != 0 )
317 {
318 PrintError( "Could not load buffer %s; Keeping Source",aPtr->sBuf.srcName );
319 goto _insert_exit;
320 }
321
322 if( HaveEqualSize( &aPtr->pano, TrPtr->dest ))
323 {
324
325 // At this point we have two valid, equally sized images
326 // Do Colour Correction on one or both images
327 DoColorCorrection( TrPtr->dest, &aPtr->pano, aPtr->sBuf.colcorrect & 3);
328
329 if( merge( TrPtr->dest , &aPtr->pano, aPtr->sBuf.feather, TrPtr->mode & _show_progress, aPtr->sBuf.seam ) != 0 )
330 {
331 PrintError( "Error merging images. Keeping Source" );
332 }
333 }
334 myfree( (void**)aPtr->pano.data );
335 } // src != 0
336
337 if( *(aPtr->sBuf.destName) != 0 ) // save buffer image
338 {
339 if( SaveBufImage( TrPtr->dest, aPtr->sBuf.destName ) != 0 )
340 PrintError( "Could not save to Buffer. Most likely your disk is full");
341 }
342 } // Tr.success
343
344
345 if( TrPtr->success == 0 && ! (TrPtr->mode & _destSupplied) )
346 myfree( (void**)TrPtr->dest->data );
347
348 _insert_exit:
349 if( aPtr->im.format ==_fisheye_circ && aPtr->im.cP.cutFrame ) // There is a cropped source image;
350 {
351 if( ImCrop.data != NULL )
352 myfree( (void**) ImCrop.data );
353 TrPtr->src = theSrc;
354 }
355
356 }
357 break;
358
359 case _extract:
360
361 if( aPtr->im.width == 0 )
362 {
363 aPtr->im.width = 500 ;
364 }
365 if( aPtr->im.height == 0 )
366 {
367 aPtr->im.height = aPtr->im.width * 4 / 5;
368 }
369
370 // Set pano-params to src-image irrespective of prefs
371 aPtr->pano.width = TrPtr->src->width; // width of panorama
372 aPtr->pano.height = TrPtr->src->height; // height of panorama
373
374 CopyPosition( TrPtr->src, &(aPtr->pano) );
375 // addAlpha( TrPtr->src );
376
377 if( *(aPtr->sBuf.destName) != 0 ) // save buffer image
378 {
379 if( SaveBufImage( TrPtr->src, aPtr->sBuf.destName ) != 0 )
380 PrintError( "Could not save Buffer Image. Most likely your disk is full");
381 }
382
383 // Set up Image Structure in TrPtr struct
384
385
386 destheight = aPtr->im.height;
387 destwidth = aPtr->im.width;
388
389
390 if( SetDestImage( TrPtr, destwidth, destheight) != 0)
391 {
392 PrintError("Could not allocate %ld bytes",TrPtr->dest->dataSize );
393 TrPtr->success = 0;
394 return;
395 }
396
397 CopyPosition( TrPtr->dest, &(aPtr->im) );
398
399 TrPtr->mode |= _honor_valid;
400 if( aPtr->pano.hfov == 360.0 )
401 TrPtr->mode |= _wrapX;
402
403 aPtr->mode = prefs->mode; // For checkparam
404 ExtractStill( TrPtr, aPtr );
405
406
407 if( TrPtr->success == 0 && ! (TrPtr->mode & _destSupplied))
408 myfree( (void**)TrPtr->dest->data );
409 break;
410
411 case _readControlPoints:
412 {
413 char *script, *newscript, cdesc[1000];
414 controlPoint cp[NUMPTS]; // List of Control points
415
416 script = LoadScript( &(prefs->scriptFile) );
417 if( script != NULL ) // We can read the scriptfile
418 {
419 newscript = (char*) malloc( strlen(script) + NUMPTS * 60 ); // One line per pair of points
420 if( newscript != NULL )
421 {
422 readControlPoints( script, cp ); // If this is the second image: get coordinates in first
423 getControlPoints( TrPtr->src, cp ); // Scan image and find control points
424 writeControlPoints( cp, cdesc ); // format control point coordinates
425
426 sprintf( newscript, "%s\n%s", script, cdesc );
427
428 if( WriteScript( newscript,&( prefs->scriptFile), 0 ) != 0 )
429 PrintError( "Could not write Scriptfile" );
430 free( newscript );
431 }
432 free( script );
433 }
434
435 }
436 TrPtr->success = 0; // Don't destroy image!
437 break;
438
439
440 case _runOptimizer:
441 // Run Optimizer; Dummy image needed but not changed
442 {
443 char* script;
444 OptInfo opt;
445 AlignInfo ainf;
446
447 script = LoadScript( &(prefs->scriptFile) );
448 if( script != NULL ) // We can read the scriptfile
449 {
450 if (ParseScript( script, &ainf ) == 0)
451 {
452 if( CheckParams( &ainf ) == 0 ) // and it seems to make sense
453 {
454 ainf.fcn = fcnPano;
455
456 // optInfo is a static variable that is used in all optimizations
457 SetGlobalPtr( &ainf ); // equivalent to optInfo = &ainf;
458
459 opt.numVars = optInfo->numParam;
460 opt.numData = optInfo->numPts;
461 opt.SetVarsToX = SetLMParams;
462 opt.SetXToVars = SetAlignParams;
463 opt.fcn = optInfo->fcn;
464 *opt.message = 0;
465 RunLMOptimizer( &opt );
466 optInfo->data = opt.message;
467 WriteResults( script, &(prefs->scriptFile), optInfo, distSquared ,
468 ( TrPtr->mode & 7 ) != _usedata );
469 }
470 DisposeAlignInfo( &ainf ); // These were allocated by 'ParseScript()'
471 }
472 free( script );
473 }
474 }
475
476 TrPtr->success = 0; // Don't destroy Dummy image!
477 break;
478 default:
479 TrPtr->success = 0;
480 break;
481
482 }
483 }
484
485
486 // Make a pano in TrPtr->dest (must be allocated and all set!)
487 // using parameters in aPrefs (ignore image parameters in TrPtr !)
488
MakePano(TrformStr * TrPtr,aPrefs * aP)489 void MakePano( TrformStr *TrPtr, aPrefs *aP )
490 {
491 MyMakePano( TrPtr, aP, 1 );
492 }
493
494
495 /*This function was added by Kekus Digital on 18/9/2002.
496 This function takes the parameter 'imageNum' which repesents the index
497 of the image that has to be converted.*/
MyMakePano(TrformStr * TrPtr,aPrefs * aP,int imageNum)498 void MyMakePano( TrformStr *TrPtr, aPrefs *aP, int imageNum )
499 {
500 struct MakeParams mp,mpinv;
501 fDesc stack[15], fD; // Parameters for execute
502 fDesc invstack[15], finvD; // Invers Parameters for execute
503 void *morph[3];
504
505 int i,k, kstart, kend, color;
506
507 TrPtr->success = 1;
508
509 if( CheckMakeParams( aP) != 0)
510 {
511 TrPtr->success = 0;
512 return;
513 }
514
515
516 if( isColorSpecific( &(aP->im.cP) ) ) // Color dependent
517 {
518 kstart = 1; kend = 4;
519 }
520 else // Color independent
521 {
522 kstart = 0; kend = 1;
523 }
524
525 for( k = kstart; k < kend; k++ )
526 {
527 color = k-1; if( color < 0 ) color = 0;
528 SetMakeParams( stack, &mp, &(aP->im) , &(aP->pano), color );
529 SetInvMakeParamsCorrect( invstack, &mpinv, &(aP->im) , &(aP->pano), color );
530
531 if( aP->nt > 0 ) // Morphing requested
532 {
533 morph[0] = (void*)aP->td;
534 morph[1] = (void*)aP->ts;
535 morph[2] = (void*)&aP->nt;
536
537 i=0; while( stack[i].func != NULL && i<14 ) i++;
538 if( i!=14 )
539 {
540 for(i=14; i>0; i--)
541 {
542 memcpy( &stack[i], &stack[i-1], sizeof( fDesc ));
543 }
544 stack[0].func = tmorph;
545 stack[0].param = (void*)morph;
546 }
547 }
548
549 if( TrPtr->success != 0)
550 {
551 fD.func = execute_stack_new; fD.param = stack;
552 finvD.func = execute_stack_new; finvD.param = invstack;
553
554 transFormEx( TrPtr, &fD , &finvD , k, imageNum );
555 }
556 }
557 }
558
559 // Extract image from pano in TrPtr->src
560 // using parameters in prefs (ignore image parameters
561 // in TrPtr)
562
ExtractStill(TrformStr * TrPtr,aPrefs * aP)563 void ExtractStill( TrformStr *TrPtr , aPrefs *aP )
564 {
565 struct MakeParams mp,mpinv;
566 fDesc stack[15], fD; // Parameters for execute
567 fDesc stackinv[15], fDinv; // Invers Parameters for execute
568
569 int k, kstart, kend, color;
570
571 TrPtr->success = 1;
572
573 if( CheckMakeParams( aP) != 0)
574 {
575 TrPtr->success = 0;
576 return;
577 }
578
579
580 if( isColorSpecific( &(aP->im.cP) ) ) // Color dependent
581 {
582 kstart = 1; kend = 4;
583 }
584 else // Color independent
585 {
586 kstart = 0; kend = 1;
587 }
588
589 for( k = kstart; k < kend; k++ )
590 {
591 color = k-1; if( color < 0 ) color = 0;
592 SetInvMakeParamsCorrect( stack, &mp, &(aP->im), &(aP->pano), color );
593 SetMakeParams( stackinv, &mpinv, &(aP->im), &(aP->pano), color );
594
595 if( TrPtr->success != 0)
596 {
597 fD.func = execute_stack_new; fD.param = stack;
598 fDinv.func = execute_stack_new; fDinv.param = stackinv;
599 transFormEx( TrPtr, &fD, &fDinv, k, 1 );
600 }
601 }
602 }
603
604
605 // Set Makeparameters depending on adjustprefs, color and source image
606
SetMakeParams(struct fDesc * stack,struct MakeParams * mp,Image * im,Image * pn,int color)607 void SetMakeParams( struct fDesc *stack, struct MakeParams *mp, Image *im , Image *pn, int color )
608 {
609 int i;
610 double a,b; // field of view in rad
611 double tx,ty, tpara; // temporary variables
612 /* Joost Nieuwenhuijse, 3 feb 2005: Fix for cropping bug
613 If a script containing the 'C' crop parameter was stitched by PTStitcher,
614 it would fail if the cropping area is partially outside the source image.
615
616 For 'inside' cropping, PTStitcher apparently pre-crops the images, such that
617 *im contains the cropped area of the source image.
618 For 'outside' cropping, PTStitcher apparently does nothing. The cropping area
619 is stored in im->selection, and im->cp.cutFrame is set, but this information
620 was not used at all.
621
622 This is fixed here: All processing is now done based on the width&height of the
623 cropped area (instead of the width&height of the image). And an additional horizontal
624 and vertical offset are added to compensate for the shift of the center of the
625 crop area relative to the center of the image.
626 */
627 int image_selection_width=im->width;
628 int image_selection_height=im->height;
629 mp->im = im;
630 mp->pn = pn;
631 if(im->cP.horizontal)
632 {
633 mp->horizontal=im->cP.horizontal_params[color];
634 }
635 else
636 {
637 mp->horizontal=0;
638 }
639 if(im->cP.vertical)
640 {
641 mp->vertical=im->cP.vertical_params[color];
642 }
643 else
644 {
645 mp->vertical=0;
646 }
647 if( (im->selection.left != 0) || (im->selection.top != 0) || (im->selection.bottom != 0) || (im->selection.right != 0) )
648 {
649 if(im->cP.cutFrame)
650 {
651 image_selection_width = im->selection.right - im->selection.left;
652 image_selection_height = im->selection.bottom - im->selection.top;
653 mp->horizontal += (im->selection.right + im->selection.left - (int32_t)im->width)/2.0;
654 mp->vertical += (im->selection.bottom + im->selection.top - (int32_t)im->height)/2.0;
655 }
656 }
657
658 a = DEG_TO_RAD( im->hfov ); // field of view in rad
659 b = DEG_TO_RAD( pn->hfov );
660
661 SetMatrix( - DEG_TO_RAD( im->pitch ),
662 0.0,
663 - DEG_TO_RAD( im->roll ),
664 mp->mt,
665 0 );
666
667 #if 0
668 switch (pn->format)
669 {
670 case _rectilinear:
671 mp->distance = (double) pn->width / (2.0 * tan(b/2.0));
672 if(im->format == _rectilinear) // rectilinear image
673 {
674 mp->scale[0] = ((double)pn->hfov / im->hfov) *
675 (a /(2.0 * tan(a/2.0))) * ((double)image_selection_width/(double) pn->width)
676 * 2.0 * tan(b/2.0) / b;
677
678 }
679 else // pamoramic or fisheye image
680 {
681 mp->scale[0] = ((double)pn->hfov / im->hfov) * ((double)image_selection_width/ (double) pn->width)
682 * 2.0 * tan(b/2.0) / b;
683 }
684 break;
685 case _equirectangular:
686 case _fisheye_ff:
687 case _panorama:
688 case _mercator:
689 case _sinusoidal:
690 // horizontal pixels per degree
691 mp->distance = ((double) pn->width) / b;
692 if(im->format == _rectilinear) // rectilinear image
693 {
694 mp->scale[0] = ((double)pn->hfov / im->hfov) * (a /(2.0 * tan(a/2.0))) * ((double)image_selection_width)/ ((double) pn->width);
695 }
696 else // pamoramic or fisheye image
697 {
698 mp->scale[0] = ((double)pn->hfov / im->hfov) * ((double)image_selection_width)/ ((double) pn->width);
699 }
700 break;
701 case _stereographic:
702 case _trans_mercator:
703 default:
704 break;
705 }
706 mp->scale[1] = mp->scale[0];
707
708 // printf("\nOrig params: mp->distance: %lf, mp->scale: %lf\n\n", mp->distance, mp->scale[0]);
709 #endif
710
711 /* Pablo d'Angelo, April 2006.
712 * Added more output projection types. Broke mp->distance and mp->scale factor calculation
713 * into separate parts, making it easier to add new projection types
714 */
715 // calculate distance
716 switch (pn->format)
717 {
718 case _rectilinear:
719 mp->distance = (double) pn->width / (2.0 * tan(b/2.0));
720 break;
721 case _equirectangular:
722 case _fisheye_ff:
723 case _fisheye_circ:
724 case _panorama:
725 case _lambert:
726 case _mercator:
727 case _millercylindrical:
728 case _sinusoidal:
729 case _mirror:
730 // horizontal pixels per degree
731 mp->distance = ((double) pn->width) / b;
732 break;
733 case _panini:
734 tpara = 1;
735 panini_erect(b/2.0, 0.0, &tx, &ty, & tpara);
736 mp->distance = pn->width/(2.0*tx);
737 break;
738 case _equipanini:
739 tpara = 1;
740 equipanini_erect(b/2.0, 0.0, &tx, &ty, & tpara);
741 mp->distance = pn->width/(2.0*tx);
742 break;
743 case _panini_general:
744 // call setup_panini_general() to set distanceparam
745 pn->precomputedCount = 0; // clear old settings
746 setup_panini_general( mp );
747 // should abort now if it returns NULL
748 break;
749 case _architectural:
750 tpara = 1;
751 arch_erect(b/2.0, 0.0, &tx, &ty, & tpara);
752 mp->distance = pn->width/(2.0*tx);
753 break;
754 case _lambertazimuthal:
755 tpara = 1;
756 lambertazimuthal_erect(b/2.0, 0.0, &tx, &ty, & tpara);
757 mp->distance = pn->width/(2.0*tx);
758 break;
759 case _hammer:
760 tpara = 1;
761 hammer_erect(b/2.0, 0.0, &tx, &ty, & tpara);
762 mp->distance = pn->width/(2.0*tx);
763 break;
764 case _stereographic:
765 tpara = 1;
766 stereographic_erect(b/2.0, 0.0, &tx, &ty, & tpara);
767 mp->distance = pn->width/(2.0*tx);
768 break;
769 case _trans_mercator:
770 tpara = 1;
771 transmercator_erect(b/2.0, 0.0, &tx, &ty, &tpara);
772 mp->distance = pn->width/(2.0*tx);
773 break;
774 case _albersequalareaconic:
775 mp->distance = 1.0;
776 //albersequalareaconic_erect(1.924913116, -PI/2.0, &tx, &ty, mp); //b/2.0
777 albersequalareaconic_distance(&tx, mp);
778 mp->distance = pn->width/(2.0*tx);
779 break;
780 case _equisolid:
781 mp->distance = (double) pn->width / (4.0 * sin(b/4.0));
782 break;
783 case _orthographic:
784 mp->distance = (double) pn->width / (2.0 * sin(b/2.0));
785 break;
786 case _thoby:
787 mp->distance = (double) pn->width / (2.0 * THOBY_K1_PARM * sin(b * THOBY_K2_PARM/2.0));
788 break;
789 case _biplane:
790 biplane_distance(pn->width,b,mp);
791 break;
792 case _triplane:
793 triplane_distance(pn->width,b,mp);
794 break;
795 default:
796 // unknown
797 PrintError ("SetMakeParams: Unsupported panorama projection");
798 // no way to report an error back to the caller...
799 mp->distance = 1;
800 break;
801 }
802
803 // calculate final scaling factor, that reverses the mp->distance
804 // scaling and applies the required output scaling factor
805 // printf("im format %d\n", im->format);
806 switch (im->format)
807 {
808 case _rectilinear:
809 // calculate distance for this projection
810 mp->scale[0] = (double) image_selection_width / (2.0 * tan(a/2.0)) / mp->distance;
811 break;
812 case _equirectangular:
813 case _panorama:
814 case _fisheye_ff:
815 case _fisheye_circ:
816 case _mercator:
817 case _sinusoidal:
818 mp->scale[0] = ((double) image_selection_width) / a / mp->distance;
819 break;
820 case _equisolid:
821 case _mirror:
822 mp->scale[0] = (double) image_selection_width / (4.0 * sin(a/4.0)) / mp->distance;
823 break;
824 case _orthographic:
825 {
826 //generate monotonic scale function to help optimizer
827 int t=(int)ceil((a-PI)/(2.0*PI));
828 mp->scale[0] = (double) image_selection_width / (2.0 * (2 * t + pow(-1.0, t) * sin(a/2.0))) / mp->distance;
829 };
830 break;
831 case _thoby:
832 mp->scale[0] = (double) image_selection_width / (2.0 * THOBY_K1_PARM * sin(a * THOBY_K2_PARM /2.0)) / mp->distance;
833 break;
834 case _stereographic:
835 mp->scale[0] = (double) image_selection_width / (4.0 * tan(a/4.0)) / mp->distance;
836 break;
837 default:
838 PrintError ("SetMakeParams: Unsupported input image projection");
839 // no way to report an error back to the caller...
840 mp->scale[0] = 1;
841 break;
842 }
843 mp->scale[1] = mp->scale[0];
844
845 // printf("new params: mp->distance: %lf, mp->scale: %lf\n\n", mp->distance, mp->scale[0]);
846
847 mp->shear[0] = im->cP.shear_x / image_selection_height;
848 mp->shear[1] = im->cP.shear_y / image_selection_width;
849 mp->rot[0] = mp->distance * PI; // 180 in screenpoints
850 mp->rot[1] = -im->yaw * mp->distance * PI / 180.0; // rotation angle in screenpoints
851
852 mp->tilt[0] = DEG_TO_RAD(im->cP.tilt_x);
853 mp->tilt[1] = DEG_TO_RAD(im->cP.tilt_y);
854 mp->tilt[2] = DEG_TO_RAD(im->cP.tilt_z);
855 mp->tilt[3] = im->cP.tilt_scale;
856
857 mp->trans[0] = im->cP.trans_x;
858 mp->trans[1] = im->cP.trans_y;
859 mp->trans[2] = im->cP.trans_z;
860 mp->trans[3] = DEG_TO_RAD(im->cP.trans_yaw);
861 mp->trans[4] = DEG_TO_RAD(im->cP.trans_pitch);
862
863 mp->test[0] = im->cP.test_p0;
864 mp->test[1] = im->cP.test_p1;
865 mp->test[2] = im->cP.test_p2;
866 mp->test[3] = im->cP.test_p3;
867
868 // panoAdjustPrintMakeParams("SetmakeParms", mp, im);
869
870 mp->perspect[0] = (void*)(mp->mt);
871 mp->perspect[1] = (void*)&(mp->distance);
872
873 for(i=0; i<4; i++)
874 mp->rad[i] = im->cP.radial_params[color][i];
875 mp->rad[5] = im->cP.radial_params[color][4];
876
877 if( (im->cP.correction_mode & 3) == correction_mode_radial )
878 mp->rad[4] = ( (double)( image_selection_width < image_selection_height ? image_selection_width : image_selection_height) ) / 2.0;
879 else
880 mp->rad[4] = ((double) image_selection_height) / 2.0;
881
882 // Joost: removed, see above
883 // mp->horizontal = im->cP.horizontal_params[color];
884 // mp->vertical = im->cP.vertical_params[color];
885
886 i = 0;
887
888
889 // Building the stack
890 //
891 // - Convert from panorama projection to equirectangular
892 // - Rotate horizontally
893 // - Convert to spherical from equirectangular
894 // - Apply perspective correction (pitch and roll) in spherical coordinates
895 // - Convert to image format (rectilinear, pano, equirectangular)
896 // - Scale output image
897 // - Do radial correction
898 // - Do tilt
899 // - Do vertical shift
900 // - Do horizontal shift
901 // - Do shear
902
903
904 //////////////////////////////////////////////////////////////////////
905 // Convert from output projection to spherical coordinates
906 //
907 if(pn->format == _rectilinear) // rectilinear panorama
908 {
909 SetDesc(stack[i], erect_rect, &(mp->distance) ); i++; // Convert rectilinear to equirect
910 }
911 else if(pn->format == _panorama)
912 {
913 SetDesc(stack[i], erect_pano, &(mp->distance) ); i++; // Convert panoramic to equirect
914 }
915 else if(pn->format == _fisheye_circ || pn->format == _fisheye_ff)
916 {
917 // the sphere coordinates are actually equivalent to the equidistant fisheye projection
918 SetDesc(stack[i], erect_sphere_tp, &(mp->distance) ); i++; // Convert fisheye to equirect
919 }
920 else if(pn->format == _equisolid)
921 {
922 SetDesc(stack[i], sphere_tp_equisolid, &(mp->distance) ); i++; // Convert fisheye equisolid to spherical
923 SetDesc(stack[i], erect_sphere_tp, &(mp->distance) ); i++; // Convert spherical to equirect
924 }
925 else if(pn->format == _mirror)
926 {
927 SetDesc(stack[i], sphere_cp_mirror, &(mp->distance) ); i++; // Convert mirror to spherical
928 SetDesc(stack[i], erect_sphere_cp, &(mp->distance) ); i++; // Convert spherical to equirect
929 }
930 else if(pn->format == _orthographic)
931 {
932 SetDesc(stack[i], sphere_tp_orthographic, &(mp->distance) ); i++; // Convert fisheye orthographic to spherical
933 SetDesc(stack[i], erect_sphere_tp, &(mp->distance) ); i++; // Convert spherical to equirect
934 }
935 else if(pn->format == _thoby)
936 {
937 SetDesc(stack[i], sphere_tp_thoby, &(mp->distance) ); i++; // Convert thoby to spherical
938 SetDesc(stack[i], erect_sphere_tp, &(mp->distance) ); i++; // Convert spherical to equirect
939 }
940 else if(pn->format == _mercator)
941 {
942 SetDesc(stack[i], erect_mercator, &(mp->distance) ); i++; // Convert mercator to equirect
943 }
944 else if(pn->format == _millercylindrical)
945 {
946 SetDesc(stack[i], erect_millercylindrical, &(mp->distance) ); i++; // Convert miller to equirect
947 }
948 else if(pn->format == _panini)
949 {
950 SetDesc(stack[i], erect_panini, &(mp->distance) ); i++; // Convert panini to sphere
951 }
952 else if(pn->format == _equipanini)
953 {
954 SetDesc(stack[i], erect_equipanini, &(mp->distance) ); i++; // Convert equipanini to sphere
955 }
956 else if(pn->format == _panini_general)
957 {
958 SetDesc(stack[i], erect_panini_general, mp ); i++; // Convert general panini to sphere
959 }
960 else if(pn->format == _architectural)
961 {
962 SetDesc(stack[i], erect_arch, &(mp->distance) ); i++; // Convert arch to sphere
963 }
964 else if(pn->format == _lambert)
965 {
966 SetDesc(stack[i], erect_lambert, &(mp->distance) ); i++; // Convert lambert to equirect
967 }
968 else if(pn->format == _lambertazimuthal)
969 {
970 SetDesc(stack[i], erect_lambertazimuthal, &(mp->distance) ); i++; // Convert lambert to equirect
971 }
972 else if(pn->format == _hammer)
973 {
974 SetDesc(stack[i], erect_hammer, &(mp->distance) ); i++; // Convert hammer to equirect
975 }
976 else if(pn->format == _trans_mercator)
977 {
978 SetDesc(stack[i], erect_transmercator, &(mp->distance) ); i++; // Convert transverse mercator to equirect
979 }
980 else if(pn->format == _stereographic)
981 {
982 SetDesc(stack[i], erect_stereographic, &(mp->distance) ); i++; // Convert stereographic to equirect
983 }
984 else if(pn->format == _sinusoidal)
985 {
986 SetDesc(stack[i], erect_sinusoidal, &(mp->distance) ); i++; // Convert sinusoidal to equirect
987 }
988 else if(pn->format == _albersequalareaconic)
989 {
990 SetDesc(stack[i], erect_albersequalareaconic, mp ); i++; // Convert albersequalareaconic to equirect
991 }
992 else if(pn->format == _biplane)
993 {
994 SetDesc(stack[i], erect_biplane, mp ); i++; // Convert biplane to equirect
995 }
996 else if(pn->format == _triplane)
997 {
998 SetDesc(stack[i], erect_triplane, mp ); i++; // Convert triplane to equirect
999 }
1000 else if(pn->format == _equirectangular)
1001 {
1002 // no conversion needed
1003 }
1004 else
1005 {
1006 PrintError("Projection type %d not supported. Assuming equirectangular", pn->format);
1007 }
1008
1009 if (im->cP.trans) {
1010 SetDesc(stack[i], plane_transfer_to_camera, mp); i++;
1011 }
1012
1013 SetDesc( stack[i], rotate_erect, mp->rot ); i++; // Rotate equirect. image horizontally
1014 SetDesc( stack[i], sphere_tp_erect, &(mp->distance) ); i++; // Convert spherical image to equirect.
1015 SetDesc( stack[i], persp_sphere, mp->perspect ); i++; // Perspective Control spherical Image
1016
1017 //////////////////////////////////////////////////////////////////////
1018 // Convert from spherical coordinates to input projection
1019 //
1020 if(im->format == _rectilinear) // rectilinear image
1021 {
1022 SetDesc(stack[i], rect_sphere_tp, &(mp->distance) ); i++; // Convert spherical to rectilinear
1023 }
1024 else if(im->format == _panorama) // pamoramic image
1025 {
1026 SetDesc(stack[i], pano_sphere_tp, &(mp->distance) ); i++; // Convert spherical to pano
1027 }
1028 else if(im->format == _equirectangular) // equirectangular image
1029 {
1030 SetDesc(stack[i], erect_sphere_tp, &(mp->distance) ); i++; // Convert spherical to equirect
1031 }
1032 else if (im->format == _fisheye_circ || im->format == _fisheye_ff)
1033 {
1034 ; // no conversion needed. It is already in spherical coordinates
1035 }
1036 else if (im->format == _mirror)
1037 {
1038 SetDesc(stack[i], mirror_sphere_tp, &(mp->distance) ); i++; // Convert spherical to mirror
1039 }
1040 else if (im->format == _stereographic)
1041 {
1042 SetDesc(stack[i], erect_sphere_tp, &(mp->distance) ); i++; // Convert spherical to equirectangular
1043 SetDesc(stack[i], stereographic_erect, &(mp->distance) ); i++; // Convert equirectangular to stereographic
1044 }
1045 else if (im->format == _orthographic)
1046 {
1047 SetDesc(stack[i], orthographic_sphere_tp, &(mp->distance) ); i++; // Convert spherical to orthographic
1048 }
1049 else if (im->format == _thoby)
1050 {
1051 SetDesc(stack[i], thoby_sphere_tp, &(mp->distance) ); i++; // Convert spherical to thoby
1052 }
1053 else if (im->format == _equisolid)
1054 {
1055 SetDesc(stack[i], erect_sphere_tp, &(mp->distance) ); i++; // Convert spherical to equirectangular
1056 SetDesc(stack[i], lambertazimuthal_erect, &(mp->distance) ); i++; // Convert equirectangular to stereographic
1057 }
1058 else
1059 {
1060 PrintError("Invalid input projection %d. Assumed fisheye.", im->format);
1061 }
1062
1063
1064 SetDesc( stack[i], resize, mp->scale ); i++; // Scale image
1065
1066 //////////////////////////////////////////////////////////////////////
1067 // Apply lens corrections
1068 //
1069
1070 if( im->cP.radial )
1071 {
1072 switch( im->cP.correction_mode & 3 )
1073 {
1074 case correction_mode_radial: SetDesc(stack[i],radial,mp->rad); i++; break;
1075 case correction_mode_vertical: SetDesc(stack[i],vertical,mp->rad); i++; break;
1076 case correction_mode_deregister:SetDesc(stack[i],deregister,mp->rad); i++; break;
1077 }
1078 }
1079 if (im->cP.tilt) {
1080 SetDesc(stack[i], tiltInverse, mp); i++;
1081 }
1082
1083 if (mp->vertical != 0.0)
1084 {
1085 SetDesc(stack[i], vert, &(mp->vertical)); i++;
1086 }
1087 if (mp->horizontal != 0.0)
1088 {
1089 SetDesc(stack[i], horiz, &(mp->horizontal)); i++;
1090 }
1091 if( im->cP.shear )
1092 {
1093 SetDesc( stack[i], shear, mp->shear ); i++;
1094 }
1095
1096 stack[i].func = (trfn)NULL;
1097
1098 // print stack for debugging
1099 #if 0
1100 printf( "Rotate params: %lg %lg\n" , mp->rot[0], mp->rot[1]);
1101 printf( "Distance : %lg\n" , mp->distance);
1102 printf( "Perspect params: %lg %lg %lg\n",a, beta , gammar );
1103 if(aP->format == _rectilinear) // rectilinear image
1104 {
1105 printf( "Rectilinear\n" );
1106 }
1107 else if (aP->format == _panorama) // pamoramic image
1108 {
1109 printf( "Panorama\n" );
1110 }
1111 else
1112 {
1113 printf( "Fisheye\n" );
1114 }
1115
1116 printf( "Scaling : %lg\n" , mp->scale[0]);
1117
1118 if( aP->correct )
1119 {
1120 printf( "Correct:\n" );
1121 if( aP->c_prefs.tilt )
1122 {
1123 printf( "Tilt: %lg\n", mp->tilt );
1124 }
1125 if( aP->c_prefs.shear )
1126 {
1127 printf( "Shear: %lg\n", mp->shear );
1128 }
1129 if ( aP->c_prefs.horizontal )
1130 {
1131 printf( "horiz:%lg\n", mp->horizontal );
1132 }
1133 if ( aP->c_prefs.vertical)
1134 {
1135 printf( "vert:%lg\n", mp->vertical );
1136 }
1137 if( aP->c_prefs.radial )
1138 {
1139 printf( "Polynomial:\n" );
1140 if( aP->c_prefs.isScanningSlit )
1141 {
1142 printf( "Scanning Slit:\n" );
1143 }
1144 else
1145 {
1146 printf( "Radial:\n" );
1147 printf( "Params: %lg %lg %lg %lg %lg\n", mp->rad[0],mp->rad[1],mp->rad[2],mp->rad[3],mp->rad[4] );
1148 }
1149 }
1150 }
1151
1152 #endif
1153 }
1154
1155
1156
1157 // Set inverse Makeparameters depending on adjustprefs, color and source image
1158 // This code is executed when optimization
SetInvMakeParams(struct fDesc * stack,struct MakeParams * mp,Image * im,Image * pn,int color)1159 void SetInvMakeParams( struct fDesc *stack, struct MakeParams *mp, Image *im , Image *pn, int color )
1160 {
1161
1162 int i;
1163 double a,b; // field of view in rad
1164 double tx,ty,tpara;
1165
1166 a = DEG_TO_RAD( im->hfov ); // field of view in rad
1167 b = DEG_TO_RAD( pn->hfov );
1168
1169 mp->im = im;
1170 mp->pn = pn;
1171
1172 SetMatrix( DEG_TO_RAD( im->pitch ),
1173 0.0,
1174 DEG_TO_RAD( im->roll ),
1175 mp->mt,
1176 1 );
1177
1178 // dangelo: added mercator, sinusoidal and stereographic projection
1179 switch (pn->format)
1180 {
1181 case _rectilinear:
1182 mp->distance = (double) pn->width / (2.0 * tan(b/2.0));
1183 break;
1184 case _equirectangular:
1185 case _fisheye_ff:
1186 case _fisheye_circ:
1187 case _panorama:
1188 case _lambert:
1189 case _mercator:
1190 case _millercylindrical:
1191 case _sinusoidal:
1192 case _mirror:
1193 // horizontal pixels per rads
1194 mp->distance = ((double) pn->width) / b;
1195 break;
1196 case _lambertazimuthal:
1197 tpara = 1;
1198 lambertazimuthal_erect(b/2.0, 0.0, &tx, &ty, & tpara);
1199 mp->distance = pn->width/(2.0*tx);
1200 break;
1201 case _hammer:
1202 tpara = 1;
1203 hammer_erect(b/2.0, 0.0, &tx, &ty, & tpara);
1204 mp->distance = pn->width/(2.0*tx);
1205 break;
1206 case _panini:
1207 tpara = 1;
1208 panini_erect(b/2.0, 0.0, &tx, &ty, & tpara);
1209 mp->distance = pn->width/(2.0*tx);
1210 break;
1211 case _equipanini:
1212 tpara = 1;
1213 equipanini_erect(b/2.0, 0.0, &tx, &ty, & tpara);
1214 mp->distance = pn->width/(2.0*tx);
1215 break;
1216 case _panini_general:
1217 // call setup_panini_general() to set distanceparam
1218 setup_panini_general( mp );
1219 // should abort now if it returns NULL
1220 break;
1221 case _architectural:
1222 tpara = 1;
1223 arch_erect(b/2.0, 0.0, &tx, &ty, & tpara);
1224 mp->distance = pn->width/(2.0*tx);
1225 break;
1226 case _stereographic:
1227 tpara = 1;
1228 stereographic_erect(b/2.0, 0.0, &tx, &ty, & tpara);
1229 mp->distance = pn->width/(2.0*tx);
1230 break;
1231 case _trans_mercator:
1232 tpara = 1;
1233 transmercator_erect(b/2.0, 0.0, &tx, &ty, & tpara);
1234 mp->distance = pn->width/(2.0*tx);
1235 break;
1236 case _albersequalareaconic:
1237 mp->distance = 1.0;
1238 //albersequalareaconic_erect(1.924913116, -PI/2.0, &tx, &ty, mp); //b/2.0
1239 albersequalareaconic_distance(&tx, mp);
1240 mp->distance = pn->width/(2.0*tx);
1241 break;
1242 case _equisolid:
1243 mp->distance = (double) pn->width / (4.0 * sin(b/4.0));
1244 break;
1245 case _orthographic:
1246 mp->distance = (double) pn->width / (2.0 * sin(b/2.0));
1247 break;
1248 case _thoby:
1249 mp->distance = (double) pn->width / (2.0 * THOBY_K1_PARM * sin(b * THOBY_K2_PARM/2.0));
1250 break;
1251 case _biplane:
1252 biplane_distance(pn->width,b,mp);
1253 break;
1254 case _triplane:
1255 triplane_distance(pn->width,b,mp);
1256 break; default:
1257 // unknown
1258 PrintError ("SetInvMakeParams: Unsupported panorama projection");
1259 // no way to report an error back to the caller...
1260 mp->distance = 1;
1261 break;
1262 }
1263
1264 // calculate final scaling factor, that reverses the mp->distance
1265 // scaling and applies the required output scaling factor
1266 switch (im->format)
1267 {
1268 case _rectilinear:
1269 // calculate distance for this projection
1270 mp->scale[0] = (double) im->width / (2.0 * tan(a/2.0)) / mp->distance;
1271 break;
1272 case _equirectangular:
1273 case _panorama:
1274 case _fisheye_ff:
1275 case _fisheye_circ:
1276 case _mercator:
1277 case _sinusoidal:
1278 mp->scale[0] = ((double) im->width) / a / mp->distance;
1279 break;
1280 case _equisolid:
1281 case _mirror:
1282 mp->scale[0] = (double) im->width / (4.0 * sin(a/4.0)) / mp->distance;
1283 break;
1284 case _orthographic:
1285 {
1286 //generate monotonic scale function to help optimizer
1287 int t=(int)ceil((a-PI)/(2.0*PI));
1288 mp->scale[0] = (double) im->width / (2.0 * (2 * t + pow(-1.0, t) * sin(a/2.0))) / mp->distance;
1289 };
1290 break;
1291 case _thoby:
1292 mp->scale[0] = (double) im->width / (2.0 * THOBY_K1_PARM * sin(a * THOBY_K2_PARM/2.0)) / mp->distance;
1293 break;
1294 case _stereographic:
1295 mp->scale[0] = (double) im->width / (4.0 * tan(a/4.0)) / mp->distance;
1296 break;
1297 default:
1298 PrintError ("SetInvMakeParams: Unsupported input image projection");
1299 // no way to report an error back to the caller...
1300 mp->scale[0] = 1;
1301 break;
1302 }
1303 mp->scale[1] = mp->scale[0];
1304
1305 mp->shear[0] = im->cP.shear_x / im->height;
1306 mp->shear[1] = im->cP.shear_y / im->width;
1307
1308
1309 mp->tilt[0] = DEG_TO_RAD(im->cP.tilt_x);
1310 mp->tilt[1] = DEG_TO_RAD(im->cP.tilt_y);
1311 mp->tilt[2] = DEG_TO_RAD(im->cP.tilt_z);
1312 mp->tilt[3] = im->cP.tilt_scale;
1313
1314 mp->trans[0] = im->cP.trans_x;
1315 mp->trans[1] = im->cP.trans_y;
1316 mp->trans[2] = im->cP.trans_z;
1317 mp->trans[3] = DEG_TO_RAD(im->cP.trans_yaw);
1318 mp->trans[4] = DEG_TO_RAD(im->cP.trans_pitch);
1319
1320 mp->test[0] = im->cP.test_p0;
1321 mp->test[1] = im->cP.test_p1;
1322 mp->test[2] = im->cP.test_p2;
1323 mp->test[3] = im->cP.test_p3;
1324
1325 // panoAdjustPrintMakeParams("Inverse 20",mp,im);
1326
1327 mp->scale[0] = 1.0 / mp->scale[0];
1328 mp->scale[1] = mp->scale[0];
1329 mp->horizontal = -im->cP.horizontal_params[color];
1330 mp->vertical = -im->cP.vertical_params[color];
1331 for(i=0; i<4; i++)
1332 mp->rad[i] = im->cP.radial_params[color][i];
1333 mp->rad[5] = im->cP.radial_params[color][4];
1334
1335 switch( im->cP.correction_mode & 3 )
1336 {
1337 case correction_mode_radial: mp->rad[4] = ((double)(im->width < im->height ? im->width : im->height) ) / 2.0;break;
1338 case correction_mode_vertical:
1339 case correction_mode_deregister: mp->rad[4] = ((double) im->height) / 2.0;break;
1340 }
1341
1342 mp->rot[0] = mp->distance * PI; // 180 in screenpoints
1343 mp->rot[1] = im->yaw * mp->distance * PI / 180.0; // rotation angle in screenpoints
1344
1345 mp->perspect[0] = (void*)(mp->mt);
1346 mp->perspect[1] = (void*)&(mp->distance);
1347
1348 // panoAdjustPrintMakeParams("Invert 30",mp,im);
1349
1350 i = 0; // Stack counter
1351
1352
1353 if( im->cP.shear )
1354 {
1355 SetDesc( stack[i],shearInv, mp->shear ); i++;
1356 }
1357
1358 if ( im->cP.horizontal )
1359 {
1360 SetDesc(stack[i],horiz, &(mp->horizontal)); i++;
1361 }
1362
1363 if ( im->cP.vertical)
1364 {
1365 SetDesc(stack[i],vert, &(mp->vertical)); i++;
1366 }
1367
1368 if( im->cP.tilt )
1369 {
1370 SetDesc( stack[i],tiltForward, mp ); i++;
1371 }
1372
1373 // Perform radial correction
1374
1375 if( im->cP.radial )
1376 {
1377 switch( im->cP.correction_mode & 3)
1378 {
1379 case correction_mode_radial: SetDesc(stack[i],inv_radial,mp->rad); i++; break;
1380 case correction_mode_vertical: SetDesc(stack[i],inv_vertical,mp->rad); i++; break;
1381 case correction_mode_deregister: break;
1382 }
1383 }
1384
1385 SetDesc( stack[i], resize, mp->scale ); i++; // Scale image
1386
1387 // printf("values %d %d\n", i, im->format);
1388
1389
1390 if(im->format == _rectilinear) // rectilinear image
1391 {
1392 SetDesc(stack[i], sphere_tp_rect, &(mp->distance) ); i++; // Convert rectilinear to spherical
1393 }
1394 else if (im->format == _panorama) // pamoramic image
1395 {
1396 SetDesc(stack[i], sphere_tp_pano, &(mp->distance) ); i++; // Convert panoramic to spherical
1397 }
1398 else if (im->format == _equirectangular) // equirectangular image
1399 {
1400 SetDesc(stack[i], sphere_tp_erect, &(mp->distance) ); i++; // Convert equirectangular to spherical
1401 }
1402 else if (im->format == _mirror) // Mirror image
1403 {
1404 SetDesc(stack[i], sphere_tp_mirror, &(mp->distance) ); i++; // Convert mirror to spherical
1405 }
1406 else if (im->format == _equisolid) // Fisheye equisolid image
1407 {
1408
1409 SetDesc(stack[i], erect_lambertazimuthal, &(mp->distance) ); i++; // Convert lambert to equirectangular
1410 SetDesc(stack[i], sphere_tp_erect, &(mp->distance) ); i++; // Convert equirectangular to spherical
1411 //SetDesc(stack[i], sphere_tp_equisolid, &(mp->distance) ); i++; // Convert equisolid to spherical
1412 }
1413 else if (im->format == _orthographic) // Fisheye orthographic image
1414 {
1415 SetDesc(stack[i], sphere_tp_orthographic, &(mp->distance) ); i++; // Convert orthographic to spherical
1416 }
1417 else if (im->format == _thoby) // thoby projected image
1418 {
1419 SetDesc(stack[i], sphere_tp_thoby, &(mp->distance) ); i++; // Convert thoby to spherical
1420 }
1421 else if (im->format == _stereographic) // Fisheye stereographic image
1422 {
1423 SetDesc(stack[i], erect_stereographic, &(mp->distance) ); i++; // Convert stereographic to spherical
1424 SetDesc(stack[i], sphere_tp_erect, &(mp->distance) ); i++; // Convert equirectangular to spherical
1425 }
1426
1427 // printf("values %d %d\n", i, im->format);
1428 SetDesc( stack[i], persp_sphere, mp->perspect ); i++; // Perspective Control spherical Image
1429 SetDesc( stack[i], erect_sphere_tp, &(mp->distance) ); i++; // Convert spherical image to equirect.
1430 SetDesc( stack[i], rotate_erect, mp->rot ); i++; // Rotate equirect. image horizontally
1431
1432 if( im->cP.trans )
1433 {
1434 SetDesc( stack[i], plane_transfer_from_camera, mp ); i++;
1435 }
1436
1437 // THESE ARE ALL FORWARD transforms
1438 if(pn->format == _rectilinear) // rectilinear panorama
1439 {
1440 SetDesc(stack[i], rect_erect, &(mp->distance) ); i++; // Convert equirectangular to rectilinear
1441 }
1442 else if(pn->format == _panorama)
1443 {
1444 SetDesc(stack[i], pano_erect, &(mp->distance) ); i++; // Convert equirectangular to Cylindrical panorama
1445 }
1446 else if(pn->format == _fisheye_circ || pn->format == _fisheye_ff )
1447 {
1448 SetDesc(stack[i], sphere_tp_erect, &(mp->distance) ); i++; // Convert equirectangular to spherical
1449 }
1450 else if(pn->format == _mercator)
1451 {
1452 SetDesc(stack[i], mercator_erect, &(mp->distance) ); i++; // Convert equirectangular to mercator
1453 }
1454 else if(pn->format == _millercylindrical)
1455 {
1456 SetDesc(stack[i], millercylindrical_erect, &(mp->distance) ); i++; // Convert equirectangular to miller cylindrical
1457 }
1458 else if(pn->format == _panini)
1459 {
1460 SetDesc(stack[i], panini_erect, &(mp->distance) ); i++; // Convert panini to sphere
1461 }
1462 else if(pn->format == _equipanini)
1463 {
1464 SetDesc(stack[i], equipanini_erect, &(mp->distance) ); i++; // Convert equi panini to sphere
1465 }
1466 else if(pn->format == _panini_general)
1467 {
1468 SetDesc(stack[i], panini_general_erect, mp ); i++; // Convert general panini to sphere
1469 }
1470 else if(pn->format == _architectural)
1471 {
1472 SetDesc(stack[i], arch_erect, &(mp->distance) ); i++; // Convert arch to sphere
1473 }
1474 else if(pn->format == _lambert)
1475 {
1476 SetDesc(stack[i], lambert_erect, &(mp->distance) ); i++; // Convert equirectangular to lambert
1477 }
1478 else if(pn->format == _lambertazimuthal)
1479 {
1480 SetDesc(stack[i], lambertazimuthal_erect, &(mp->distance) ); i++; // Convert equirectangular to lambert azimuthal
1481 }
1482 else if(pn->format == _hammer)
1483 {
1484 SetDesc(stack[i], hammer_erect, &(mp->distance) ); i++; // Convert equirectangular to hammer
1485 }
1486 else if(pn->format == _trans_mercator)
1487 {
1488 SetDesc(stack[i], transmercator_erect, &(mp->distance) ); i++; // Convert equirectangular to transverse mercator
1489 }
1490 else if(pn->format == _mirror)
1491 {
1492 SetDesc(stack[i], mirror_erect, &(mp->distance) ); i++; // Convert equirectangular to mirror
1493 }
1494 else if(pn->format == _stereographic)
1495 {
1496 SetDesc(stack[i], stereographic_erect, &(mp->distance) ); i++; // Convert equirectangular to stereographic
1497 }
1498 else if(pn->format == _sinusoidal)
1499 {
1500 SetDesc(stack[i], sinusoidal_erect, &(mp->distance) ); i++; // Convert equirectangular to sinusoidal
1501 }
1502 else if(pn->format == _albersequalareaconic)
1503 {
1504 SetDesc(stack[i], albersequalareaconic_erect, mp ); i++; // Convert equirectangular to albersequalareaconic
1505 }
1506 else if(pn->format == _equisolid )
1507 {
1508 SetDesc(stack[i], sphere_tp_erect, &(mp->distance) ); i++; // Convert equirectangular to spherical
1509 SetDesc(stack[i], equisolid_sphere_tp, &(mp->distance) ); i++; // Convert spherical to equisolid
1510 }
1511 else if(pn->format == _orthographic )
1512 {
1513 SetDesc(stack[i], sphere_tp_erect, &(mp->distance) ); i++; // Convert equirectangular to spherical
1514 SetDesc(stack[i], orthographic_sphere_tp, &(mp->distance) ); i++; // Convert spherical to orthographic
1515 }
1516 else if(pn->format == _thoby )
1517 {
1518 SetDesc(stack[i], sphere_tp_erect, &(mp->distance) ); i++; // Convert equirectangular to spherical
1519 SetDesc(stack[i], thoby_sphere_tp, &(mp->distance) ); i++; // Convert spherical to thoby
1520 }
1521 else if(pn->format == _biplane)
1522 {
1523 SetDesc(stack[i], biplane_erect, mp ); i++; // Convert equirectangular to biplane
1524 }
1525 else if(pn->format == _triplane)
1526 {
1527 SetDesc(stack[i], triplane_erect, mp ); i++; // Convert equirectangular to biplane
1528 } else if(pn->format == _equirectangular)
1529 {
1530 // no conversion needed
1531 }
1532 else
1533 {
1534 PrintError("Projection type %d not supported, using equirectangular", pn->format);
1535 }
1536
1537 stack[i].func = (trfn)NULL;
1538 }
1539
SetInvMakeParamsCorrect(struct fDesc * stack,struct MakeParams * mp,Image * im,Image * pn,int color)1540 void SetInvMakeParamsCorrect( struct fDesc *stack, struct MakeParams *mp, Image *im , Image *pn, int color )
1541 {
1542 /* Thomas Rauscher, Sep 2005: Transfered the changes of Joost Nieuwenhuijse for MakeParams
1543 to the inverse function. This has broken the optimizer, now there are two functions.
1544 */
1545
1546 Image imSel; /* create a tempory copy of the image to manipulate */
1547 memcpy( &imSel, im, sizeof(Image));
1548
1549 if(im->cP.horizontal)
1550 {
1551 mp->horizontal = im->cP.horizontal_params[color];
1552 }
1553 else
1554 {
1555 mp->horizontal = 0;
1556 }
1557
1558 if(im->cP.vertical)
1559 {
1560 mp->vertical = im->cP.vertical_params[color];
1561 }
1562 else
1563 {
1564 mp->vertical = 0;
1565 }
1566
1567 if( (im->selection.left != 0) || (im->selection.top != 0) || (im->selection.bottom != 0) || (im->selection.right != 0) )
1568 {
1569 if(im->cP.cutFrame)
1570 {
1571 imSel.width = im->selection.right - im->selection.left;
1572 imSel.height = im->selection.bottom - im->selection.top;
1573
1574 mp->horizontal += (im->selection.right + im->selection.left - im->width)/2.0;
1575 mp->vertical += (im->selection.bottom + im->selection.top - im->height)/2.0;
1576
1577 imSel.cP.horizontal_params[color] = mp->horizontal;
1578 imSel.cP.vertical_params[color] = mp->vertical;
1579 }
1580 }
1581
1582 SetInvMakeParams( stack, mp, &imSel, pn, color );
1583 }
1584
1585
1586
1587 // Add an alpha channel to the image, assuming rectangular or circular shape
1588 // subtract frame
1589
addAlpha(Image * im)1590 void addAlpha( Image *im ){
1591 register int x,y,c1;
1592 int framex, framey;
1593 register unsigned char *src;
1594
1595 src = *(im->data);
1596 framex = 0; framey = 0;
1597
1598 if( im->cP.cutFrame ){
1599 if( im->cP.frame < 0 || im->cP.fwidth < 0 || im->cP.fheight < 0 ){ // Use supplied alpha channel
1600 return;
1601 }
1602
1603 if( im->cP.frame != 0 ){
1604 framex = (im->width/2 > im->cP.frame ? im->cP.frame : 0);
1605 framey = (im->height/2 > im->cP.frame ? im->cP.frame : 0);
1606 }
1607 else{
1608 if( im->width > im->cP.fwidth )
1609 framex = (im->width - im->cP.fwidth) / 2;
1610 if( im->height > im->cP.fheight )
1611 framey = (im->height - im->cP.fheight) / 2;
1612 }
1613 }
1614
1615
1616 if( im->bitsPerPixel == 32 || im->bitsPerPixel == 64 ) // leave 24/48 bit images unchanged
1617 {
1618 if( im->format != _fisheye_circ ) // Rectangle valid
1619 {
1620 int yend = im->height - framey;
1621 int xend = im->width - framex;
1622
1623 if( im->bitsPerPixel == 32 )
1624 {
1625 if( 0 != framey || 0 != framex )
1626 {
1627 for(y = 0; y < im->height; y++)
1628 {
1629 c1 = y * im->bytesPerLine;
1630 for(x = 0; x < im->width; x++)
1631 src[ c1 + 4 * x ] = 0;
1632 }
1633 }
1634 for(y = framey; y < yend; y++)
1635 {
1636 c1 = y * im->bytesPerLine;
1637 for(x = framex; x < xend; x++)
1638 src[ c1 + 4 * x ] = UCHAR_MAX;
1639 }
1640 }
1641 else // im->bitsPerPixel == 64
1642 {
1643 if( 0 != framey || 0 != framex )
1644 {
1645 for(y = 0; y < im->height; y++)
1646 {
1647 c1 = y * im->bytesPerLine;
1648 for(x = 0; x < im->width; x++)
1649 *((uint16_t*)(src + c1 + 8 * x )) = 0;
1650 }
1651 }
1652 for(y = framey; y < yend; y++)
1653 {
1654 c1 = y * im->bytesPerLine;
1655 for(x = framex; x < xend; x++)
1656 *((uint16_t*)(src + c1 + 8 * x )) = USHRT_MAX;
1657 }
1658 }
1659 }
1660 else if( im->format == _fisheye_circ ) // Circle valid
1661 {
1662 int topCircle = ( im->height - im->width ) / 2; // top of circle
1663 int botCircle = topCircle + im->width ; // bottom of circle
1664 int r = ( im->width / 2 ); // radius of circle
1665 int x1, x2, h;
1666
1667 if( im->bitsPerPixel == 32 )
1668 {
1669 for(y = 0; y < im->height ; y++)
1670 {
1671 if( (y < topCircle) || (y > botCircle) ) // Always invalid
1672 {
1673 for(x = 0; x < im->width; x++)
1674 src[y * im->bytesPerLine + 4 * x] = 0;
1675 }
1676 else
1677 {
1678 h = y - im->height/2;
1679 if( h*h > r*r ) h = r;
1680
1681 x1 = (int) (r - sqrt( r*r - h*h ));
1682 if( x1 < 0 ) x1 = 0;
1683 x2 = (int) (r + sqrt( r*r - h*h ));
1684 if( x2 > im->width ) x2 = im->width;
1685
1686 for(x = 0; x < x1; x++)
1687 src[y * im->bytesPerLine + 4 * x] = 0;
1688 for(x = x1; x < x2; x++)
1689 src[y * im->bytesPerLine + 4 * x] = UCHAR_MAX;
1690 for(x = x2; x < im->width; x++)
1691 src[y * im->bytesPerLine + 4 * x] = 0;
1692 }
1693 }
1694 }
1695 else // im->bitsPerPixel == 64
1696 {
1697 for(y = 0; y < im->height ; y++)
1698 {
1699 if( (y < topCircle) || (y > botCircle) ) // Always invalid
1700 {
1701 for(x = 0; x < im->width; x++)
1702 *((uint16_t*)(src + y * im->bytesPerLine + 8 * x)) = 0;
1703 }
1704 else
1705 {
1706 h = y - im->height/2;
1707 if( h*h > r*r ) h = r;
1708
1709 x1 = (int) (r - sqrt( r*r - h*h ));
1710 if( x1 < 0 ) x1 = 0;
1711 x2 = (int) (r + sqrt( r*r - h*h ));
1712 if( x2 > im->width ) x2 = im->width;
1713
1714 for(x = 0; x < x1; x++)
1715 *((uint16_t*)(src + y * im->bytesPerLine + 8 * x)) = 0;
1716 for(x = x1; x < x2; x++)
1717 *((uint16_t*)(src + y * im->bytesPerLine + 8 * x)) = USHRT_MAX;
1718 for(x = x2; x < im->width; x++)
1719 *((uint16_t*)(src + y * im->bytesPerLine + 8 * x)) = 0;
1720 }
1721 }
1722 }
1723 } // mode
1724 } // pixelsize
1725 }
1726
1727
1728 // Angular Distance for control point "num".
1729 // Function distSphere computes an exact angular distance and the
1730 // corresponding components in longitude/latitude directions.
1731 // These are returned in a slightly strange manner (distance as the
1732 // function result, components in a global array) to avoid changing the
1733 // calling sequence of distSphere, which might unnecessarily break
1734 // other code that we don't know about.
1735
1736 double distanceComponents[2];
1737
distSphere(int num)1738 double distSphere( int num ){
1739 double x, y ; // Coordinates of control point in panorama
1740 double w2, h2;
1741 int j;
1742 Image sph;
1743 int n[2];
1744 struct MakeParams mp;
1745 struct fDesc stack[15];
1746 CoordInfo b[2];
1747 CoordInfo cp;
1748 double lat[2], lon[2]; // latitude & longitude
1749 double dlon;
1750 double dangle;
1751 double dist;
1752 double radiansToPixelsFactor;
1753
1754 // Factor to convert angular error in radians to equivalent in pixels
1755
1756 radiansToPixelsFactor = optInfo->pano.width / (optInfo->pano.hfov * (PI/180.0));
1757
1758 // Get image position in imaginary spherical image
1759
1760 SetImageDefaults( &sph );
1761
1762 sph.width = 360;
1763 sph.height = 180;
1764 sph.format = _equirectangular;
1765 sph.hfov = 360.0;
1766
1767 n[0] = optInfo->cpt[num].num[0];
1768 n[1] = optInfo->cpt[num].num[1];
1769
1770 // Calculate coordinates using equirectangular mapping to get longitude/latitude
1771
1772 for(j=0; j<2; j++){
1773 SetInvMakeParams( stack, &mp, &optInfo->im[ n[j] ], &sph, 0 );
1774
1775 h2 = (double)optInfo->im[ n[j] ].height / 2.0 - 0.5;
1776 w2 = (double)optInfo->im[ n[j] ].width / 2.0 - 0.5;
1777
1778
1779 execute_stack_new( (double)optInfo->cpt[num].x[j] - w2, // cartesian x-coordinate src
1780 (double)optInfo->cpt[num].y[j] - h2, // cartesian y-coordinate src
1781 &x, &y, stack);
1782
1783 x = DEG_TO_RAD( x );
1784 y = DEG_TO_RAD( y ) + PI/2.0;
1785
1786 // x is now in the range -PI to +PI, and y is 0 to PI
1787 lat[j] = y;
1788 lon[j] = x;
1789
1790 b[j].x[0] = sin(x) * sin( y );
1791 b[j].x[1] = cos( y );
1792 b[j].x[2] = - cos(x) * sin(y);
1793 }
1794
1795 dlon = lon[0]-lon[1];
1796 if (dlon < -PI) dlon += 2.0*PI;
1797 if (dlon > PI) dlon -= 2.0*PI;
1798 distanceComponents[0] = (dlon*sin(0.5*(lat[0]+lat[1]))) * radiansToPixelsFactor;
1799 distanceComponents[1] = (lat[0]-lat[1]) * radiansToPixelsFactor;
1800
1801 // The original acos formulation (acos(SCALAR_PRODUCT(&b[0],&b[1]))
1802 // is inaccurate for angles near 0, because it essentially requires finding eps
1803 // based on the value of 1-eps^2/2. The asin formulation is much more
1804 // accurate under these conditions.
1805
1806 CROSS_PRODUCT(&b[0],&b[1],&cp);
1807 dangle = asin(ABS_VECTOR(&cp));
1808 if (SCALAR_PRODUCT(&b[0],&b[1]) < 0.0) dangle = PI - dangle;
1809 dist = dangle * radiansToPixelsFactor;
1810
1811 // Diagnostics to help debug various calculation errors.
1812 // Do not delete this code --- it has been needed surprisingly often.
1813 #if 0
1814 { double olddist;
1815 olddist = acos( SCALAR_PRODUCT( &b[0], &b[1] ) ) * radiansToPixelsFactor;
1816 // if (adjustLogFile != 0 && abs(dist-olddist) > 1.0) {
1817 if (adjustLogFile != 0 && num < 5) {
1818 fprintf(adjustLogFile,"***** DIST ***** dCoord = %g %g, lonlat0 = %g %g, lonlat1 = %g %g, dist=%g, olddist=%g, sumDcoordSq=%g, distSq=%g\n",
1819 distanceComponents[0],distanceComponents[1],lon[0],lat[0],lon[1],lat[1],dist,olddist,
1820 distanceComponents[0]*distanceComponents[0]+distanceComponents[1]*distanceComponents[1],dist*dist);
1821 }
1822 }
1823 #endif
1824
1825 return dist;
1826 }
1827
1828
1829
pt_getXY(int n,double x,double y,double * X,double * Y)1830 void pt_getXY(int n, double x, double y, double *X, double *Y){
1831 struct MakeParams mp;
1832 struct fDesc stack[15];
1833 double h2,w2;
1834
1835 SetInvMakeParams( stack, &mp, &optInfo->im[ n ], &optInfo->pano, 0 );
1836 h2 = (double)optInfo->im[ n ].height / 2.0 - 0.5;
1837 w2 = (double)optInfo->im[ n ].width / 2.0 - 0.5;
1838
1839
1840 execute_stack_new( x - w2, y - h2, X, Y, stack);
1841 }
1842
1843 // Return distance of points from a line
1844 // The line through the two farthest apart points is calculated
1845 // Returned is the sum distance squared of the other two points from the line
distsqLine(int N0,int N1)1846 double distsqLine(int N0, int N1){
1847 double x[4],y[4], del, delmax, A, B, C, mu, d0, d1;
1848 int n0, n1, n2=-1, n3=-1, i, k;
1849
1850 pt_getXY(optInfo->cpt[N0].num[0], (double)optInfo->cpt[N0].x[0], (double)optInfo->cpt[N0].y[0], &x[0], &y[0]);
1851 pt_getXY(optInfo->cpt[N0].num[1], (double)optInfo->cpt[N0].x[1], (double)optInfo->cpt[N0].y[1], &x[1], &y[1]);
1852 pt_getXY(optInfo->cpt[N1].num[0], (double)optInfo->cpt[N1].x[0], (double)optInfo->cpt[N1].y[0], &x[2], &y[2]);
1853 pt_getXY(optInfo->cpt[N1].num[1], (double)optInfo->cpt[N1].x[1], (double)optInfo->cpt[N1].y[1], &x[3], &y[3]);
1854
1855 delmax = 0.0;
1856 n0 = 0; n1 = 1;
1857
1858 for(i=0; i<4; i++){
1859 for(k=i+1; k<4; k++){
1860 del = (x[i]-x[k])*(x[i]-x[k])+(y[i]-y[k])*(y[i]-y[k]);
1861 if(del>delmax){
1862 n0=i; n1=k; delmax=del;
1863 }
1864 }
1865 }
1866 if(delmax==0.0) return 0.0;
1867
1868 for(i=0; i<4; i++){
1869 if(i!= n0 && i!= n1){
1870 n2 = i;
1871 break;
1872 }
1873 }
1874 for(i=0; i<4; i++){
1875 if(i!= n0 && i!= n1 && i!=n2){
1876 n3 = i;
1877 }
1878 }
1879
1880
1881 A=y[n1]-y[n0]; B=x[n0]-x[n1]; C=y[n0]*(x[n1]-x[n0])-x[n0]*(y[n1]-y[n0]);
1882
1883 mu=1.0/sqrt(A*A+B*B);
1884
1885 d0 = (A*x[n2]+B*y[n2]+C)*mu;
1886 d1 = (A*x[n3]+B*y[n3]+C)*mu;
1887 distanceComponents[0] = d0;
1888 distanceComponents[1] = d1;
1889
1890 return d0*d0 + d1*d1;
1891
1892 }
1893
1894
1895 // Calculate the distance of Control Point "num" between two images
1896 // in final pano. This is the old distSquared.
1897
rectDistSquared(int num)1898 double rectDistSquared( int num )
1899 {
1900 double x[2], y[2]; // Coordinates of control point in panorama
1901 double w2, h2;
1902 int j, n[2];
1903 double result;
1904
1905 struct MakeParams mp;
1906 struct fDesc stack[15];
1907
1908
1909
1910 n[0] = optInfo->cpt[num].num[0];
1911 n[1] = optInfo->cpt[num].num[1];
1912
1913 // Calculate coordinates x/y in panorama
1914
1915 for(j=0; j<2; j++)
1916 {
1917 SetInvMakeParams( stack, &mp, &optInfo->im[ n[j] ], &optInfo->pano, 0 );
1918
1919 h2 = (double)optInfo->im[ n[j] ].height / 2.0 - 0.5;
1920 w2 = (double)optInfo->im[ n[j] ].width / 2.0 - 0.5;
1921
1922
1923 execute_stack_new( (double)optInfo->cpt[num].x[j] - w2, // cartesian x-coordinate src
1924 (double)optInfo->cpt[num].y[j] - h2, // cartesian y-coordinate src
1925 &x[j], &y[j], stack);
1926 // test to check if inverse works
1927 #if 0
1928 {
1929 double xt, yt;
1930 struct MakeParams mtest;
1931 struct fDesc stacktest[15];
1932 SetMakeParams( stacktest, &mtest, &optInfo->im[ n[j] ], &optInfo->pano, 0 );
1933 execute_stack_new( x[j], // cartesian x-coordinate src
1934 y[j], // cartesian y-coordinate src
1935 &xt, &yt, stacktest);
1936
1937 printf("x= %lg, y= %lg, xb = %lg, yb = %lg \n", optInfo->cpt[num].x[j], optInfo->cpt[num].y[j], xt+w2, yt+h2);
1938
1939 }
1940 #endif
1941 }
1942
1943
1944 // printf("Coordinates 0: %lg:%lg 1: %lg:%lg\n",x[0] + g->pano->width/2,y[0]+ g->pano->height/2, x[1] + g->pano->width/2,y[1]+ g->pano->height/2);
1945
1946
1947 // take care of wrapping and points at edge of panorama
1948
1949 if( optInfo->pano.hfov == 360.0 )
1950 {
1951 double delta = abs( x[0] - x[1] );
1952
1953 if( delta > optInfo->pano.width / 2 )
1954 {
1955 if( x[0] < x[1] )
1956 x[0] += optInfo->pano.width;
1957 else
1958 x[1] += optInfo->pano.width;
1959 }
1960 }
1961
1962
1963 switch( optInfo->cpt[num].type ) // What do we want to optimize?
1964 {
1965 case 1: // x difference
1966 result = ( x[0] - x[1] ) * ( x[0] - x[1] );
1967 break;
1968 case 2: // y-difference
1969 result = ( y[0] - y[1] ) * ( y[0] - y[1] );
1970 break;
1971 default:
1972 result = ( y[0] - y[1] ) * ( y[0] - y[1] ) + ( x[0] - x[1] ) * ( x[0] - x[1] ); // square of distance
1973 distanceComponents[0] = y[0] - y[1];
1974 distanceComponents[1] = x[0] - x[1];
1975
1976 break;
1977 }
1978
1979
1980 return result;
1981 }
1982
1983 /// huber() is an M-Estimator function. Using an M-Estimator might
1984 /// work better if the control points contain outliers (eg. from autopano).
1985 /// this implementation accepts normal, non squared errors, and return non-squared errors,
1986 /// contrary to the definition in the literature (where the square is
1987 /// included in the function)
1988
1989 static double fcnPanoHuberSigma = 0; // sigma for Huber M-estimator. 0 disables M-estimator
1990
setFcnPanoHuberSigma(double sigma)1991 void setFcnPanoHuberSigma(double sigma)
1992 {
1993 fcnPanoHuberSigma = sigma;
1994 }
1995
huber(double x,double sigma)1996 double huber(double x, double sigma)
1997 {
1998 if (abs(x) < sigma)
1999 return x;
2000 else
2001 return sqrt(2.0*sigma*abs(x) - sigma*sigma);
2002 }
2003
2004
2005
2006 /// (function distSquared2 has been removed -- it was unused and redundant)
2007
2008 /// EvaluateControlPointErrorAndComponents is the central point-of-contact
2009 /// for determining the error for a specified control point pair.
2010 ///
2011 /// Arguments are
2012 /// num identifies the control point pair
2013 /// *errptr returns a single error measure (distance)
2014 /// errComponent[*] returns two components of that error, as nearly orthogonal
2015 /// as possible
2016 /// Return value is a success flag: 0 = successful, otherwise not.
2017
EvaluateControlPointErrorAndComponents(int num,double * errptr,double errComponent[2])2018 int EvaluateControlPointErrorAndComponents ( int num, double *errptr, double errComponent[2]) {
2019 int j;
2020 int result;
2021 switch(optInfo->cpt[num].type){
2022 case 0: // normal control points
2023 // Jim May 2004.
2024 // Optimizing cylindrical and rectilinear projection by calculating
2025 // distance error in pixel coordinates of the rendered image.
2026 // When using angular (spherical) distance for these projections,
2027 // larger errors are generated the further control points are from
2028 // the center.
2029 // In theory by optimizing in pixel coordinates all errors will be
2030 // distributed over the image. This is true.
2031 // In practice I have found that optimize large field of view
2032 // rectilinear projection images failed to resolve nicely if the
2033 // parameters were not very close to start with. I leave the
2034 // code here for others to play with and maybe get better results.
2035 /* if(optInfo->pano.format == _rectilinear || g->pano.format == _panorama)
2036 {
2037 *errptr = sqrt(rectDistSquared(num));
2038 errComponent[0] = distanceComponents[0];
2039 errComponent[1] = distanceComponents[1];
2040 result = 0;
2041 break;
2042 }
2043 else // _equirectangular, fisheye, spherical, mirror
2044 { */
2045 *errptr = distSphere(num);
2046 errComponent[0] = distanceComponents[0];
2047 errComponent[1] = distanceComponents[1];
2048 result = 0;
2049 break;
2050 //}
2051 case 1: // vertical
2052 case 2: // horizontal
2053 *errptr = sqrt(rectDistSquared(num));
2054 errComponent[0] = *errptr;
2055 errComponent[1] = 0.0;
2056 result = 0;
2057 break;
2058 default:// t+ controls = lines = sets of two control point pairs
2059 *errptr = 0.0; // in case there is no second pair
2060 errComponent[0] = 0.0;
2061 errComponent[1] = 0.0;
2062 result = 1;
2063 for(j=0; j<optInfo->numPts; j++){
2064 if(j!=num && optInfo->cpt[num].type == optInfo->cpt[j].type){
2065 *errptr = sqrt(distsqLine(num,j));
2066 // errComponent[0] = *errptr;
2067 // errComponent[1] = 0.0;
2068 errComponent[0] = distanceComponents[0];
2069 errComponent[1] = distanceComponents[1];
2070 result = 0;
2071 break;
2072 }
2073 }
2074 break;
2075 }
2076 return result;
2077 }
2078
2079 /// This distSquared is a convenience function, to be passed into
2080 /// WriteResults in the usual fashion, to avoid having to change
2081 /// other codes that call WriteResults. It replaces the old distSquared,
2082 /// which has been renamed rectDistSquared and is now used only
2083 /// internally by EvaluateControlPointErrorAndComponents.
2084
distSquared(int num)2085 double distSquared (int num ) {
2086 double result;
2087 double junk[2];
2088 EvaluateControlPointErrorAndComponents (num, &result, junk);
2089 return result*result;
2090 }
2091
2092 /// Function fcnPano calculates a vector of control points errors,
2093 /// for use during optimization. See lmdif.c and optimize.c for
2094 /// a description of its arguments. The helper functions that appear here
2095 /// allow to control the new capability, while preserving also the
2096 /// old calling sequences.
2097
2098 int fcnPanoNperCP = 1; // number of functions per control point, 1 or 2
2099
setFcnPanoNperCP(int i)2100 void setFcnPanoNperCP (int i) {
2101 fcnPanoNperCP = i;
2102 }
2103
getFcnPanoNperCP()2104 int getFcnPanoNperCP() {
2105 return fcnPanoNperCP;
2106 }
2107
setFcnPanoDoNotInitAvgFov()2108 void setFcnPanoDoNotInitAvgFov() { // must be called after iflag = -100 call to fcnPano
2109 needInitialAvgFov = 0;
2110 }
2111
forceFcnPanoReinitAvgFov()2112 void forceFcnPanoReinitAvgFov() { // applies to next call to fcnPano
2113 needInitialAvgFov = 1;
2114 }
2115
fcnPano(int m,int n,double x[],double fvec[],int * iflag)2116 int fcnPano(int m, int n, double x[], double fvec[], int *iflag)
2117 {
2118 int i;
2119 static int numIt;
2120 double result;
2121 int iresult;
2122 double junk;
2123 double junk2[2];
2124
2125 if( *iflag == -100 ){ // reset
2126 numIt = 0;
2127 needInitialAvgFov = 1;
2128 infoDlg ( _initProgress, "Optimizing Variables" );
2129 #if ADJUST_LOGGING_ENABLED
2130 if ((adjustLogFile = fopen(ADJUST_LOG_FILENAME,"a")) <= 0) {
2131 PrintError("Cannot Open Log File");
2132 adjustLogFile = 0;
2133 }
2134 #endif
2135 return 0;
2136 }
2137 if( *iflag == -99 ){ //
2138 infoDlg ( _disposeProgress, "" );
2139 if (adjustLogFile != 0) {
2140 result = 0.0;
2141 for( i=0; i < m; i++)
2142 {
2143 result += fvec[i]*fvec[i] ;
2144 }
2145 result = sqrt( result/ (double)m ) * sqrt((double)fcnPanoNperCP); // to approximate total distance vs dx, dy
2146 fprintf(adjustLogFile,"At iflag=-99 (dispose), NperCP=%d, m=%d, n=%d, err = %g, x= \n",
2147 fcnPanoNperCP,m,n,result);
2148 for (i=0; i<n; i++) {
2149 fprintf(adjustLogFile,"\t%20.10g",x[i]);
2150 }
2151 fprintf(adjustLogFile,"\n fvec = \n");
2152 for (i=0; i<m; i++) {
2153 fprintf(adjustLogFile,"\t%20.10g",fvec[i]);
2154 if (((i+1) % fcnPanoNperCP) == 0) fprintf(adjustLogFile,"\n");
2155 }
2156 fprintf(adjustLogFile,"\n");
2157 fclose(adjustLogFile);
2158 }
2159 return 0;
2160 }
2161
2162
2163 if( *iflag == 0 )
2164 {
2165 char message[256];
2166
2167 result = 0.0;
2168 for( i=0; i < m; i++)
2169 {
2170 result += fvec[i]*fvec[i] ;
2171 }
2172 result = sqrt( result/ (double)m ) * sqrt((double)fcnPanoNperCP); // to approximate total distance vs dx, dy
2173
2174 sprintf( message,"Strategy %d\nAverage (rms) distance between Controlpoints \nafter %d iteration(s): %25.15g units", getFcnPanoNperCP(), numIt,result);//average);
2175 numIt += 1; // 10;
2176 if( !infoDlg ( _setProgress,message ) )
2177 *iflag = -1;
2178
2179 if (adjustLogFile != 0) {
2180 fprintf(adjustLogFile,"At iteration %d, iflag=0 (print), NperCP=%d, m=%d, n=%d, err = %g, x= \n",
2181 numIt,fcnPanoNperCP,m,n,result);
2182 for (i=0; i<n; i++) {
2183 fprintf(adjustLogFile,"\t%20.10g",x[i]);
2184 }
2185 fprintf(adjustLogFile,"\n fvec = \n");
2186 for (i=0; i<m; i++) {
2187 fprintf(adjustLogFile,"\t%20.10g",fvec[i]);
2188 if (((i+1) % fcnPanoNperCP) == 0) fprintf(adjustLogFile,"\n");
2189 }
2190 fprintf(adjustLogFile,"\n");
2191 fflush(adjustLogFile);
2192 }
2193
2194 return 0;
2195 }
2196
2197 // Set Parameters
2198
2199 SetAlignParams( x ) ;
2200
2201 if (needInitialAvgFov) {
2202 initialAvgFov = avgfovFromSAP;
2203 needInitialAvgFov = 0;
2204 if (adjustLogFile != 0) {
2205 fprintf(adjustLogFile,"setting initialAvgFov = %g\n",initialAvgFov);
2206 fflush(adjustLogFile);
2207 }
2208 }
2209
2210 if (adjustLogFile != 0) {
2211 fprintf(adjustLogFile,"entering fcnPano, m=%d, n=%d, initialAvgFov=%g, avgfovFromSAP=%g, x = \n",
2212 m,n, initialAvgFov,avgfovFromSAP);
2213 for (i=0; i<n; i++) {
2214 fprintf(adjustLogFile,"\t%20.10g",x[i]);
2215 }
2216 fprintf(adjustLogFile,"\n");
2217 fflush(adjustLogFile);
2218 }
2219
2220 // Calculate distances
2221
2222 iresult = 0;
2223 for( i=0; i < optInfo->numPts; i++){
2224 if (fcnPanoNperCP == 1) {
2225 EvaluateControlPointErrorAndComponents ( i, &fvec[iresult], &junk2[0]);
2226 } else {
2227 EvaluateControlPointErrorAndComponents ( i, &junk, &fvec[iresult]);
2228 if (fcnPanoHuberSigma) {
2229 fvec[iresult] = huber(fvec[iresult], fcnPanoHuberSigma);
2230 fvec[iresult+1] = huber(fvec[iresult+1], fcnPanoHuberSigma);
2231 }
2232 }
2233
2234 // Field-of-view stabilization. Applying here means that the
2235 // errors seen by the optimizer may be different from those finally
2236 // reported, by the same factor for all errors. This introduces
2237 // the possibility of confusion for people who are paying really
2238 // close attention to the optimizer's periodic output versus the
2239 // final result. However, it seems like the right thing to do
2240 // because then the final reported errors will correspond to the
2241 // user's settings for pano size, total fov etc.
2242
2243 if ((initialAvgFov / avgfovFromSAP) > 1.0) {
2244 fvec[iresult] *= initialAvgFov / avgfovFromSAP;
2245 }
2246 iresult += 1;
2247 if (fcnPanoNperCP == 2) {
2248 if ((initialAvgFov / avgfovFromSAP) > 1.0) {
2249 fvec[iresult] *= initialAvgFov / avgfovFromSAP;
2250 }
2251 iresult += 1;
2252 }
2253 }
2254
2255 // If not enough control points are provided, then fill out
2256 // the function vector with copies of the average error
2257 // for the actual control points.
2258
2259 result = 0.0;
2260 for (i=0; i < iresult; i++) {
2261 result += fvec[i]*fvec[i];
2262 }
2263 result = sqrt(result/(double)iresult);
2264 for (i=iresult; i < m; i++) {
2265 fvec[i] = result;
2266 }
2267
2268 if (adjustLogFile != 0) {
2269 result *= sqrt((double)fcnPanoNperCP);
2270 fprintf(adjustLogFile,"leaving fcnPano, m=%d, n=%d, err=%g, fvec = \n",m,n,result);
2271 for (i=0; i<m; i++) {
2272 fprintf(adjustLogFile,"\t%20.10g",fvec[i]);
2273 if (((i+1) % fcnPanoNperCP) == 0) fprintf(adjustLogFile,"\n");
2274 }
2275 fprintf(adjustLogFile,"\n");
2276 fflush(adjustLogFile);
2277 }
2278
2279 return 0;
2280 }
2281
2282
2283
2284 // Find Colour correcting polynomial matching the overlap of src and buf
2285 // using least square fit.
2286 // Each RGB-Channel is fitted using the relation
2287 // buf = coeff[0] * src + coeff[1]
2288 #if 1
GetColCoeff(Image * src,Image * buf,double ColCoeff[3][2])2289 void GetColCoeff( Image *src, Image *buf, double ColCoeff[3][2] ){
2290 register int x,y,c1,c2,i, numPts;
2291 double xy[3], xi[3], xi2[3], yi[3], xav[3], yav[3];
2292 register unsigned char *source, *buff;
2293 int BitsPerChannel,bpp;
2294
2295
2296
2297 GetBitsPerChannel( src, BitsPerChannel );
2298 bpp = src->bitsPerPixel/8;
2299
2300
2301 source = *(src->data);
2302 buff = *(buf->data);
2303
2304 for(i=0;i<3;i++){
2305 xy[i] = xi[i] = xi2[i] = yi[i] = 0.0;
2306 }
2307 numPts = 0;
2308
2309 if( BitsPerChannel == 8 ){
2310 for( y=2; y<src->height-2; y++){
2311 c1 = y * src->bytesPerLine;
2312 for( x=2; x<src->width-2; x++){
2313 c2 = c1 + x*bpp;
2314 if( source[c2] != 0 && buff[c2] != 0 ){ // && // In overlap region?
2315 //(source[c2] != UCHAR_MAX || buff[c2] != UCHAR_MAX)){ // above seam?
2316 if( pt_average( source+c2, src->bytesPerLine, xav, 1 ) &&
2317 pt_average( buff+c2, src->bytesPerLine, yav, 1 ) ){
2318 numPts++;
2319 for( i=0; i<3; i++){
2320 xi[i] += xav[i];
2321 yi[i] += yav[i];
2322 xi2[i] += xav[i]*xav[i];
2323 xy[i] += xav[i]*yav[i];
2324 }
2325 }
2326 }
2327 }
2328 }
2329 }else{//16
2330 for( y=1; y<src->height-1; y++){
2331 c1 = y * src->bytesPerLine;
2332 for( x=1; x<src->width-1; x++){
2333 c2 = c1 + x*bpp;
2334 if( *((uint16_t*)(source + c2)) != 0 && *((uint16_t*)(buff + c2)) != 0 ) { //&& // In overlap region?
2335 //( *((uint16_t*)(source + c2)) != USHRT_MAX || *((uint16_t*)(buff + c2)) != USHRT_MAX ) ){ // above seam?
2336 if( pt_average( source + c2, src->bytesPerLine, xav, 2 ) &&
2337 pt_average( buff + c2, src->bytesPerLine, yav, 2 )){
2338 numPts++;
2339 for( i=0; i<3; i++){
2340 xi[i] += xav[i];
2341 yi[i] += yav[i];
2342 xi2[i] += xav[i]*xav[i];
2343 xy[i] += xav[i]*yav[i];
2344 }
2345 }
2346 }
2347 }
2348 }
2349 }
2350
2351
2352 if( numPts > 0 ){
2353 for( i=0; i<3; i++){
2354 ColCoeff[i][0] = ( numPts * xy[i] - xi[i] * yi[i] ) / ( numPts * xi2[i] - xi[i]*xi[i] );
2355 ColCoeff[i][1] = ( xi2[i] * yi[i] - xy[i] * xi[i] ) / ( numPts * xi2[i] - xi[i]*xi[i] );
2356 }
2357 }else{
2358 for( i=0; i<3; i++){
2359 ColCoeff[i][0] = 1.0;
2360 ColCoeff[i][1] = 0.0;
2361 }
2362 }
2363 }
2364 #endif
2365 // Average 9 pixels
pt_average(uint8_t * pixel,int BytesPerLine,double rgb[3],int bytesPerChannel)2366 int pt_average( uint8_t* pixel, int BytesPerLine, double rgb[3], int bytesPerChannel ){
2367 int x, y, i;
2368 uint8_t *px;
2369 double sum = 1.0 + 4 * 0.5 + 8 * 0.2 + 8 * 0.1 ;//2.6;
2370 #if 0
2371 double bl[3][3] = {{ 0.1, 0.3, 0.1}, // Blurr overlap using this matrix
2372 { 0.3, 1.0, 0.3},
2373 { 0.1, 0.3, 0.1}};
2374
2375 #endif
2376 double bl[5][5] = {{ 0.0, 0.1, 0.2, 0.1, 0.0},
2377 { 0.1, 0.2, 0.5, 0.2, 0.1},
2378 { 0.2, 0.5, 1.0, 0.5, 0.2},
2379 { 0.1, 0.2, 0.5, 0.2, 0.1},
2380 { 0.0, 0.1, 0.2, 0.1, 0.0}};
2381
2382
2383 rgb[0] = rgb[1] = rgb[2] = 0.0;
2384 if( bytesPerChannel != 1 ) return -1;
2385
2386 for(y=0; y<5; y++){
2387 for(x=0; x<5; x++){
2388 px = pixel + (y-2)*BytesPerLine + x-2;
2389 if( *px == 0 ) return 0;
2390 rgb[0] += *(++px) * bl[y][x];
2391 rgb[1] += *(++px) * bl[y][x];
2392 rgb[2] += *(++px) * bl[y][x];
2393 }
2394 }
2395 for( i=0; i<3; i++) rgb[i]/=sum;
2396 return 0;
2397
2398 }
2399
2400
2401 #if 0
2402
2403 // Backup
2404
2405 // Find Colour correcting polynomial matching the overlap of src and buf
2406 // using least square fit.
2407 // Each RGB-Channel is fitted using the relation
2408 // buf = coeff[0] * src + coeff[1]
2409
2410 void GetColCoeff( Image *src, Image *buf, double ColCoeff[3][2] )
2411 {
2412 register int x,y,c1,c2,i, numPts;
2413 double xy[3], xi[3], xi2[3], yi[3];
2414 register unsigned char *source, *buff;
2415 int BitsPerChannel,bpp;
2416
2417 GetBitsPerChannel( src, BitsPerChannel );
2418 bpp = src->bitsPerPixel/8;
2419
2420
2421 source = *(src->data);
2422 buff = *(buf->data);
2423 for(i=0;i<3;i++)
2424 {
2425 xy[i] = xi[i] = xi2[i] = yi[i] = 0.0;
2426 }
2427 numPts = 0;
2428
2429 if( BitsPerChannel == 8 )
2430 {
2431 for( y=0; y<src->height; y++)
2432 {
2433 c1 = y * src->bytesPerLine;
2434 for( x=0; x<src->width; x++)
2435 {
2436 c2 = c1 + x*bpp;
2437 if( source[c2] != 0 && buff[c2] != 0 ) // In overlap region?
2438 {
2439 numPts++;
2440 for( i=0; i<3; i++)
2441 {
2442 c2++;
2443 xi[i] += (double)source[c2];
2444 yi[i] += (double)buff[c2];
2445 xi2[i] += ((double)source[c2])*((double)source[c2]);
2446 xy[i] += ((double)source[c2])*((double)buff[c2]);
2447 }
2448 }
2449 }
2450 }
2451 }
2452 else // 16
2453 {
2454 for( y=0; y<src->height; y++)
2455 {
2456 c1 = y * src->bytesPerLine;
2457 for( x=0; x<src->width; x++)
2458 {
2459 c2 = c1 + x*bpp;
2460 if( *((uint16_t*)(source + c2)) != 0 && *((uint16_t*)(buff + c2)) != 0 ) // In overlap region?
2461 {
2462 numPts++;
2463 for( i=0; i<3; i++)
2464 {
2465 c2++;
2466 xi[i] += (double) *((uint16_t*)(source + c2));
2467 yi[i] += (double) *((uint16_t*)(buff + c2));
2468 xi2[i] += ((double) *((uint16_t*)(source + c2)))*((double) *((uint16_t*)(source + c2)));
2469 xy[i] += ((double) *((uint16_t*)(source + c2)))*((double) *((uint16_t*)(buff + c2)));
2470 }
2471 }
2472 }
2473 }
2474 }
2475
2476
2477 if( numPts > 0 )
2478 {
2479 for( i=0; i<3; i++)
2480 {
2481 ColCoeff[i][0] = ( numPts * xy[i] - xi[i] * yi[i] ) / ( numPts * xi2[i] - xi[i]*xi[i] );
2482 ColCoeff[i][1] = ( xi2[i] * yi[i] - xy[i] * xi[i] ) / ( numPts * xi2[i] - xi[i]*xi[i] );
2483 }
2484 }
2485 else
2486 {
2487 for( i=0; i<3; i++)
2488 {
2489 ColCoeff[i][0] = 1.0;
2490 ColCoeff[i][1] = 0.0;
2491 }
2492 }
2493 }
2494
2495 #endif
2496
2497
2498 // Colourcorrect the image im using polynomial coefficients ColCoeff
2499 // Each RGB-Channel is corrected using the relation
2500 // new = coeff[0] * old + coeff[1]
2501
ColCorrect(Image * im,double ColCoeff[3][2])2502 void ColCorrect( Image *im, double ColCoeff[3][2] )
2503 {
2504 register int x,y, c1, c2, i;
2505 register unsigned char* data;
2506 register double result;
2507 int bpp, BitsPerChannel;
2508
2509 GetBitsPerChannel( im, BitsPerChannel );
2510 bpp = im->bitsPerPixel/8;
2511
2512 data = *(im->data);
2513
2514 if( BitsPerChannel == 8 )
2515 {
2516 for( y=0; y<im->height; y++)
2517 {
2518 c1 = y * im->bytesPerLine;
2519 for( x=0; x<im->width; x++ )
2520 {
2521 c2 = c1 + x * bpp;
2522 if( data[ c2 ] != 0 ) // Alpha channel set
2523 {
2524 for( i=0; i<3; i++)
2525 {
2526 c2++;
2527 result = ColCoeff[i][0] * data[ c2 ] + ColCoeff[i][1];
2528 DBL_TO_UC( data[ c2 ], result );
2529 }
2530 }
2531 }
2532 }
2533 }
2534 else // 16
2535 {
2536 for( y=0; y<im->height; y++)
2537 {
2538 c1 = y * im->bytesPerLine;
2539 for( x=0; x<im->width; x++ )
2540 {
2541 c2 = c1 + x * bpp;
2542 if( *((uint16_t*)(data + c2 )) != 0 ) // Alpha channel set
2543 {
2544 for( i=0; i<3; i++)
2545 {
2546 c2++;
2547 result = ColCoeff[i][0] * *((uint16_t*)(data + c2 )) + ColCoeff[i][1];
2548 DBL_TO_US( *((uint16_t*)(data + c2 )) , result );
2549 }
2550 }
2551 }
2552 }
2553 }
2554 }
2555
2556
SetAdjustDefaults(aPrefs * prefs)2557 void SetAdjustDefaults( aPrefs *prefs )
2558 {
2559
2560 prefs->magic = 50; // File validity check, must be 50
2561 prefs->mode = _insert; //
2562
2563 SetImageDefaults( &(prefs->im) );
2564 SetImageDefaults( &(prefs->pano) );
2565
2566 SetStitchDefaults( &(prefs->sBuf) );
2567
2568 memset( &(prefs->scriptFile), 0, sizeof( fullPath ) );
2569
2570 prefs->nt = 0;
2571 prefs->ts = NULL;
2572 prefs->td = NULL;
2573
2574 prefs->interpolator = _spline36;
2575 prefs->gamma = 1.0;
2576 prefs->fastStep = FAST_TRANSFORM_STEP_NONE;
2577 }
2578
2579
2580
DisposeAlignInfo(struct AlignInfo * g)2581 void DisposeAlignInfo( struct AlignInfo *g )
2582 {
2583 if(g->im != NULL) free(g->im);
2584 if(g->opt!= NULL) free(g->opt);
2585 if(g->cpt!= NULL) free(g->cpt);
2586 if(g->t != NULL) free(g->t);
2587 if(g->cim != NULL) free(g->cim);
2588 }
2589
2590
2591
2592
2593 // Set global preferences structures using LM-params
2594
SetAlignParams(double * x)2595 int SetAlignParams( double *x )
2596 {
2597 // Set Parameters
2598 int i,j,k;
2599 double sumfov = 0.0;
2600
2601 j = 0;
2602 for( i=0; i<optInfo->numIm; i++ ){
2603
2604 if( (k = optInfo->opt[i].yaw) > 0 ){
2605 if( k == 1 ){ optInfo->im[i].yaw = x[j++]; NORM_ANGLE( optInfo->im[i].yaw );
2606 }else{ optInfo->im[i].yaw = optInfo->im[k-2].yaw; }
2607 }
2608 if( (k = optInfo->opt[i].pitch) > 0 ){
2609 if( k == 1 ){ optInfo->im[i].pitch = x[j++]; NORM_ANGLE( optInfo->im[i].pitch );
2610 }else{ optInfo->im[i].pitch = optInfo->im[k-2].pitch; }
2611 }
2612 if( (k = optInfo->opt[i].roll) > 0 ){
2613 if( k == 1 ){ optInfo->im[i].roll = x[j++]; NORM_ANGLE( optInfo->im[i].roll );
2614 }else{ optInfo->im[i].roll = optInfo->im[k-2].roll; }
2615 }
2616 if( (k = optInfo->opt[i].hfov) > 0 ){
2617 if( k == 1 ){
2618 optInfo->im[i].hfov = x[j++];
2619 if( optInfo->im[i].hfov < 0.0 )
2620 optInfo->im[i].hfov = - optInfo->im[i].hfov;
2621 }else{ optInfo->im[i].hfov = optInfo->im[k-2].hfov; }
2622 }
2623 sumfov += optInfo->im[i].hfov;
2624 if( (k = optInfo->opt[i].a) > 0 ){
2625 if( k == 1 ){ optInfo->im[i].cP.radial_params[0][3] = x[j++] / C_FACTOR;
2626 }else{ optInfo->im[i].cP.radial_params[0][3] = optInfo->im[k-2].cP.radial_params[0][3];}
2627 }
2628 if( (k = optInfo->opt[i].b) > 0 ){
2629 if( k == 1 ){
2630 optInfo->im[i].cP.radial_params[0][2] = x[j++] / C_FACTOR;
2631 }else{
2632 optInfo->im[i].cP.radial_params[0][2] = optInfo->im[k-2].cP.radial_params[0][2];
2633 }
2634 }
2635 if( (k = optInfo->opt[i].c) > 0 ){
2636 if( k == 1 ){ optInfo->im[i].cP.radial_params[0][1] = x[j++] / C_FACTOR;
2637 }else{ optInfo->im[i].cP.radial_params[0][1] = optInfo->im[k-2].cP.radial_params[0][1];}
2638 }
2639 if( (k = optInfo->opt[i].d) > 0 ){
2640 if( k == 1 ){ optInfo->im[i].cP.horizontal_params[0] = x[j++];
2641 }else{ optInfo->im[i].cP.horizontal_params[0] = optInfo->im[k-2].cP.horizontal_params[0];}
2642 }
2643 if( (k = optInfo->opt[i].e) > 0 ){
2644 if( k == 1 ){ optInfo->im[i].cP.vertical_params[0] = x[j++];
2645 }else{ optInfo->im[i].cP.vertical_params[0] = optInfo->im[k-2].cP.vertical_params[0];}
2646 }
2647 // tilt
2648 if( (k = optInfo->opt[i].tiltXopt) > 0 ){
2649 if( k == 1 ){ optInfo->im[i].cP.tilt_x = x[j++]; //NORM_ANGLE_RAD(optInfo->im[i].cP.tilt_x);
2650 }else{ optInfo->im[i].cP.tilt_x = optInfo->im[k-2].cP.tilt_x;}
2651 }
2652 if( (k = optInfo->opt[i].tiltYopt) > 0 ){
2653 if( k == 1 ){ optInfo->im[i].cP.tilt_y = x[j++]; //NORM_ANGLE_RAD(optInfo->im[i].cP.tilt_y);
2654 }else{ optInfo->im[i].cP.tilt_y = optInfo->im[k-2].cP.tilt_y;}
2655 }
2656 if( (k = optInfo->opt[i].tiltZopt) > 0 ){
2657 if( k == 1 ){ optInfo->im[i].cP.tilt_z =x[j++]; //NORM_ANGLE_RAD(optInfo->im[i].cP.tilt_z);
2658 }else{ optInfo->im[i].cP.tilt_z = optInfo->im[k-2].cP.tilt_z;}
2659 }
2660 if( (k = optInfo->opt[i].tiltScaleOpt) > 0 ){
2661 if( k == 1 ) {
2662 optInfo->im[i].cP.tilt_scale = x[j++];
2663 if (optInfo->im[i].cP.tilt_scale == 0) {
2664 optInfo->im[i].cP.tilt_scale = 0.001; //make sure it never becomes zero
2665 }
2666 optInfo->im[i].cP.tilt_scale = fabs(optInfo->im[i].cP.tilt_scale);
2667 /*
2668 if (optInfo->im[i].cP.tilt_scale > 10) {
2669 optInfo->im[i].cP.tilt_scale = 10; //make sure it never gets out of control
2670 }
2671 */
2672 } else{ optInfo->im[i].cP.tilt_scale = optInfo->im[k-2].cP.tilt_scale;}
2673 }
2674 // translate
2675 if( (k = optInfo->opt[i].transXopt) > 0 ){
2676 if( k == 1 ){ optInfo->im[i].cP.trans_x = x[j++];
2677 }else{ optInfo->im[i].cP.trans_x = optInfo->im[k-2].cP.trans_x;}
2678 }
2679 if( (k = optInfo->opt[i].transYopt) > 0 ){
2680 if( k == 1 ){ optInfo->im[i].cP.trans_y = x[j++];
2681 }else{ optInfo->im[i].cP.trans_y = optInfo->im[k-2].cP.trans_y;}
2682 }
2683 if( (k = optInfo->opt[i].transZopt) > 0 ){
2684 if( k == 1 ){ optInfo->im[i].cP.trans_z = x[j++];
2685 }else{ optInfo->im[i].cP.trans_z = optInfo->im[k-2].cP.trans_z;}
2686 }
2687 if( (k = optInfo->opt[i].transYawOpt) > 0 ){
2688 if( k == 1 ) { optInfo->im[i].cP.trans_yaw = x[j++]; //NORM_ANGLE(optInfo->im[i].cP.trans_yaw);
2689 while( optInfo->im[i].cP.trans_yaw > optInfo->im[i].yaw + 80) optInfo->im[i].cP.trans_yaw -= 180.0;
2690 while( optInfo->im[i].cP.trans_yaw < optInfo->im[i].yaw - 80) optInfo->im[i].cP.trans_yaw += 180.0;
2691 } else { optInfo->im[i].cP.trans_yaw = optInfo->im[k-2].cP.trans_yaw;}
2692 }
2693 if( (k = optInfo->opt[i].transPitchOpt) > 0 ){
2694 if( k == 1 ){ optInfo->im[i].cP.trans_pitch = x[j++]; //NORM_ANGLE(optInfo->im[i].cP.trans_pitch);
2695 while( optInfo->im[i].cP.trans_pitch > optInfo->im[i].pitch + 80) optInfo->im[i].cP.trans_pitch -= 180.0;
2696 while( optInfo->im[i].cP.trans_pitch < optInfo->im[i].pitch - 80) optInfo->im[i].cP.trans_pitch += 180.0;
2697 }else{ optInfo->im[i].cP.trans_pitch = optInfo->im[k-2].cP.trans_pitch;}
2698 }
2699 // test
2700 if( (k = optInfo->opt[i].testP0opt) > 0 ){
2701 if( k == 1 ){ optInfo->im[i].cP.test_p0 = x[j++];
2702 }else{ optInfo->im[i].cP.test_p0 = optInfo->im[k-2].cP.test_p0;}
2703 }
2704 if( (k = optInfo->opt[i].testP1opt) > 0 ){
2705 if( k == 1 ){ optInfo->im[i].cP.test_p1 = x[j++];
2706 }else{ optInfo->im[i].cP.test_p1 = optInfo->im[k-2].cP.test_p1;}
2707 }
2708 if( (k = optInfo->opt[i].testP2opt) > 0 ){
2709 if( k == 1 ){ optInfo->im[i].cP.test_p2 = x[j++];
2710 }else{ optInfo->im[i].cP.test_p2 = optInfo->im[k-2].cP.test_p2;}
2711 }
2712 if( (k = optInfo->opt[i].testP3opt) > 0 ){
2713 if( k == 1 ){ optInfo->im[i].cP.test_p3 = x[j++];
2714 }else{ optInfo->im[i].cP.test_p3 = optInfo->im[k-2].cP.test_p3;}
2715 }
2716
2717 //shear
2718 if( (k = optInfo->opt[i].shear_x) > 0 ){
2719 if( k == 1 ){ optInfo->im[i].cP.shear_x = x[j++];
2720 }else{ optInfo->im[i].cP.shear_x = optInfo->im[k-2].cP.shear_x;}
2721 }
2722
2723 if( (k = optInfo->opt[i].shear_y) > 0 ){
2724 if( k == 1 ){ optInfo->im[i].cP.shear_y = x[j++];
2725 }else{ optInfo->im[i].cP.shear_y = optInfo->im[k-2].cP.shear_y;}
2726 }
2727
2728
2729 optInfo->im[i].cP.radial_params[0][0] = 1.0 - ( optInfo->im[i].cP.radial_params[0][3]
2730 + optInfo->im[i].cP.radial_params[0][2]
2731 + optInfo->im[i].cP.radial_params[0][1] ) ;
2732
2733 }
2734 avgfovFromSAP = sumfov / optInfo->numIm;
2735 if( j != optInfo->numParam )
2736 return -1;
2737 else
2738 return 0;
2739
2740 }
2741
2742 // Set LM params using global preferences structure
2743 // Change to cover range 0....1 (roughly)
2744
SetLMParams(double * x)2745 int SetLMParams( double *x )
2746 {
2747 int i,j;
2748
2749 j=0; // Counter for optimization parameters
2750
2751
2752 for( i=0; i<optInfo->numIm; i++ ){
2753 if(optInfo->opt[i].yaw == 1) // optimize alpha? 0-no 1-yes
2754 x[j++] = optInfo->im[i].yaw;
2755
2756 if(optInfo->opt[i].pitch == 1) // optimize pitch? 0-no 1-yes
2757 x[j++] = optInfo->im[i].pitch;
2758
2759 if(optInfo->opt[i].roll == 1) // optimize gamma? 0-no 1-yes
2760 x[j++] = optInfo->im[i].roll ;
2761
2762 if(optInfo->opt[i].hfov == 1) // optimize hfov? 0-no 1-yes
2763 x[j++] = optInfo->im[i].hfov ;
2764
2765 if(optInfo->opt[i].a == 1) // optimize a? 0-no 1-yes
2766 x[j++] = optInfo->im[i].cP.radial_params[0][3] * C_FACTOR;
2767
2768 if(optInfo->opt[i].b == 1) // optimize b? 0-no 1-yes
2769 x[j++] = optInfo->im[i].cP.radial_params[0][2] * C_FACTOR;
2770
2771 if(optInfo->opt[i].c == 1) // optimize c? 0-no 1-yes
2772 x[j++] = optInfo->im[i].cP.radial_params[0][1] * C_FACTOR;
2773
2774 if(optInfo->opt[i].d == 1) // optimize d? 0-no 1-yes
2775 x[j++] = optInfo->im[i].cP.horizontal_params[0] ;
2776
2777 if(optInfo->opt[i].e == 1) // optimize e? 0-no 1-yes
2778 x[j++] = optInfo->im[i].cP.vertical_params[0] ;
2779
2780 // Tilt
2781 if(optInfo->opt[i].tiltXopt == 1) { // optimize tilt_x? 0-no 1-yes
2782 x[j++] = optInfo->im[i].cP.tilt_x ;
2783 }
2784 if(optInfo->opt[i].tiltYopt == 1) // optimize tilt_y? 0-no 1-yes
2785 x[j++] = optInfo->im[i].cP.tilt_y ;
2786
2787 if(optInfo->opt[i].tiltZopt == 1) { // optimize tilt_Z? 0-no 1-yes
2788 x[j++] = optInfo->im[i].cP.tilt_z ;
2789 }
2790 if(optInfo->opt[i].tiltScaleOpt == 1) { // optimize tilt_scale? 0-no 1-yes
2791 x[j++] = optInfo->im[i].cP.tilt_scale ;
2792 }
2793 // Trans
2794 if(optInfo->opt[i].transXopt == 1) { // optimize trans_x? 0-no 1-yes
2795 x[j++] = optInfo->im[i].cP.trans_x ;
2796 }
2797 if(optInfo->opt[i].transYopt == 1) { // optimize trans_y? 0-no 1-yes
2798 x[j++] = optInfo->im[i].cP.trans_y ;
2799 }
2800 if(optInfo->opt[i].transZopt == 1) { // optimize trans_Z? 0-no 1-yes
2801 x[j++] = optInfo->im[i].cP.trans_z ;
2802 }
2803 if(optInfo->opt[i].transYawOpt == 1) { // optimize trans_yaw? 0-no 1-yes
2804 x[j++] = optInfo->im[i].cP.trans_yaw ;
2805 }
2806 if(optInfo->opt[i].transPitchOpt == 1) { // optimize trans_pitch? 0-no 1-yes
2807 x[j++] = optInfo->im[i].cP.trans_pitch ;
2808 }
2809
2810 // Test
2811 if(optInfo->opt[i].testP0opt == 1) {
2812 x[j++] = optInfo->im[i].cP.test_p0; ;
2813 }
2814 if(optInfo->opt[i].testP1opt == 1) {
2815 x[j++] = optInfo->im[i].cP.test_p1; ;
2816 }
2817 if(optInfo->opt[i].testP2opt == 1) {
2818 x[j++] = optInfo->im[i].cP.test_p2; ;
2819 }
2820 if(optInfo->opt[i].testP3opt == 1) {
2821 x[j++] = optInfo->im[i].cP.test_p3; ;
2822 }
2823
2824 if(optInfo->opt[i].shear_x == 1) // optimize shear_x? 0-no 1-yes
2825 x[j++] = optInfo->im[i].cP.shear_x ;
2826
2827 if(optInfo->opt[i].shear_y == 1) // optimize shear_y? 0-no 1-yes
2828 x[j++] = optInfo->im[i].cP.shear_y ;
2829 }
2830
2831 if( j != optInfo->numParam )
2832 return -1;
2833 else
2834 return 0;
2835
2836 }
2837
2838
2839
2840
2841
2842
2843
2844 #define DX 3
2845 #define DY 14
2846
2847 // Read Control Point Position from flag pasted into image
2848
getControlPoints(Image * im,controlPoint * cp)2849 void getControlPoints( Image *im, controlPoint *cp )
2850 {
2851 int y, x, cy,cx, bpp, r,g,b,n, nim=0, k,i,np;
2852 register unsigned char *p,*ch;
2853
2854
2855 p = *(im->data);
2856 bpp = im->bitsPerPixel/8;
2857 if( bpp == 4 )
2858 {
2859 r = 1; g = 2; b = 3;
2860 }
2861 else if( bpp == 3 )
2862 {
2863 r = 0; g = 1; b = 2;
2864 }
2865 else
2866 {
2867 PrintError("Can't read ControlPoints from images with %d Bytes per Pixel", bpp);
2868 return;
2869 }
2870
2871 np = 0;
2872 for(y=0; y<im->height; y++)
2873 {
2874 cy = y * im->bytesPerLine;
2875 for(x=0; x<im->width; x++)
2876 {
2877 cx = cy + bpp * x;
2878 if( p[ cx + r ] == 0 && p[ cx + g ] == 255 && p[ cx + b ] == 0 &&
2879 p[ cx + bpp + r ] == 255 && p[ cx + bpp + g ] == 0 && p[ cx + bpp + b ] == 0 &&
2880 p[ cx + 2*bpp + r ] == 0 && p[ cx + 2*bpp + g ] == 0 && p[ cx + 2*bpp + b ] == 255 &&
2881 p[ cx - bpp + r ] == 0 && p[ cx - bpp + g ] == 0 && p[ cx - bpp + b ] == 0 )
2882 {
2883 if(p[cx + 3*bpp + r ] == 0 && p[ cx + 3*bpp + g ] == 255 && p[ cx + 3*bpp + b ] == 255)
2884 { // Control Point
2885 ch = &(p[cx + 4*bpp + r ]);
2886 n = 0;
2887 while( ch[0] == 255 && ch[1] == 0 && ch[2] == 0 )
2888 {
2889 n++;
2890 ch += bpp;
2891 }
2892 if( n >= 0 )
2893 {
2894 k = 0;
2895 if( cp[n].num[0] != -1 )
2896 k = 1;
2897 cp[n].x[k] = x + DX;
2898 cp[n].y[k] = y + DY;
2899 np++;
2900 }
2901 }
2902 else if(p[cx+3*bpp +r] == 255 && p[ cx + 3*bpp + g ] == 255 && p[ cx + 3*bpp + b ] == 0)
2903 { // Image number
2904 ch = &(p[cx + 4*bpp + r ]);
2905 n = 0;
2906 while( ch[0] == 255 && ch[1] == 0 && ch[2] == 0 )
2907 {
2908 n++;
2909 ch += bpp;
2910 }
2911 if( n >= 0 )
2912 {
2913 nim = n;
2914 }
2915 }
2916 }
2917 }
2918 }
2919 k = 0;
2920 if( cp[0].num[0] != -1 )
2921 k = 1;
2922 for(i=0; i<np; i++)
2923 cp[i].num[k] = nim;
2924
2925
2926 }
2927
2928
2929 // Write Control Point coordinates into script
2930
writeControlPoints(controlPoint * cp,char * cdesc)2931 void writeControlPoints( controlPoint *cp,char* cdesc )
2932 {
2933 int i;
2934 char line[80];
2935
2936 *cdesc = 0;
2937 for(i=0; i<NUMPTS && cp[i].num[0] != -1; i++)
2938 {
2939 //sprintf( line, "c n%d N%d x%d y%d X%d Y%d\n", cp[i].num[0], cp[i].num[1],
2940 sprintf( line, "c n%d N%d x%lf y%lf X%lf Y%lf\n", cp[i].num[0], cp[i].num[1],
2941 cp[i].x[0], cp[i].y[0],
2942 cp[i].x[1], cp[i].y[1]);
2943 strcat( cdesc, line );
2944 }
2945 }
2946
2947
SetStitchDefaults(struct stitchBuffer * sBuf)2948 void SetStitchDefaults( struct stitchBuffer *sBuf)
2949 {
2950 *sBuf->srcName = 0;
2951 *sBuf->destName = 0;
2952 sBuf->feather = 10;
2953 sBuf->colcorrect = 0;
2954 sBuf->seam = _middle;
2955 sBuf->psdOpacity = 255;
2956 sBuf->psdBlendingMode = PSD_NORMAL;
2957 }
2958
SetOptDefaults(optVars * opt)2959 void SetOptDefaults( optVars *opt )
2960 {
2961 opt->hfov = opt->yaw = opt->pitch = opt->roll = 0;
2962 opt->a = opt->b = opt->c = opt->d = opt->e = 0;
2963 opt->tiltXopt = opt->tiltYopt = opt->tiltZopt = opt->tiltScaleOpt = 0;
2964 opt->transXopt = opt->transYopt = opt->transZopt = opt->transYawOpt = opt->transPitchOpt = 0;
2965 opt->testP0opt = opt->testP1opt = opt->testP2opt = opt->testP3opt = 0;
2966 opt->shear_x = opt->shear_y = 0;
2967 }
2968
DoColorCorrection(Image * im1,Image * im2,int mode)2969 void DoColorCorrection( Image *im1, Image *im2, int mode )
2970 {
2971 double ColCoeff [3][2];
2972 int i;
2973
2974 switch( mode )
2975 {
2976 case 0:
2977 break; // no correction
2978 case 1: // Correct im1
2979 GetColCoeff( im1, im2, ColCoeff );
2980 ColCorrect( im1, ColCoeff );
2981 break;
2982 case 2: // Correct im2
2983 GetColCoeff( im1, im2, ColCoeff );
2984 // Invert coefficients
2985 for( i = 0; i<3; i++)
2986 {
2987 ColCoeff[i][1] = - ColCoeff[i][1] / ColCoeff[i][0];
2988 ColCoeff[i][0] = 1.0/ColCoeff[i][0];
2989 }
2990 ColCorrect( im2, ColCoeff );
2991 break;
2992 case 3: // Correct both halfs
2993 GetColCoeff( im1, im2, ColCoeff );
2994 for(i = 0; i<3; i++)
2995 {
2996 ColCoeff[i][1] = ColCoeff[i][1] / 2.0 ;
2997 ColCoeff[i][0] = (ColCoeff[i][0] + 1.0 ) / 2.0;
2998 }
2999 ColCorrect( im1, ColCoeff );
3000 for(i = 0; i<3; i++)
3001 {
3002 ColCoeff[i][1] = - ColCoeff[i][1] / ColCoeff[i][0];
3003 ColCoeff[i][0] = 1.0 / ColCoeff[i][0];
3004 }
3005 ColCorrect( im2, ColCoeff );
3006 break;
3007 default: break;
3008 } // switch
3009 }
3010
3011
3012 // Do some checks on Optinfo structure and reject if obviously nonsense
3013
CheckParams(AlignInfo * g)3014 int CheckParams( AlignInfo *g )
3015 {
3016 int i;
3017 int err = -1;
3018 char *errmsg[] = {
3019 "No Parameters to optimize",
3020 "No input images",
3021 "No Feature Points",
3022 "Image width must be positive",
3023 "Image height must be positive",
3024 "Field of View must be positive",
3025 "Field of View must be smaller than 180 degrees in rectilinear Images",
3026 "Unsupported Image Format (must be 0,1,2,3,4,7,8,10,14 or 19)",
3027 "Panorama Width must be positive",
3028 "Panorama Height must be positive",
3029 "Field of View must be smaller than 180 degrees in rectilinear Panos",
3030 "Unsupported Panorama Format",
3031 "Control Point Coordinates must be positive",
3032 "Invalid Image Number in Control Point Descriptions",
3033 };
3034
3035 if( g->numParam == 0 ) err = 0;
3036 if( g->numIm == 0 ) err = 1;
3037 if( g->numPts == 0 ) err = 2;
3038
3039 // Check images
3040
3041
3042
3043 for( i=0; i<g->numIm; i++)
3044 {
3045
3046 if (g->im[i].cP.tilt_scale == 0) {
3047 PrintError("Image [%d] has tilt_scale equal zero [%d]\n", i, g->im[i].cP.tilt_scale);
3048 return -1;
3049 }
3050
3051 if( g->im[i].width <= 0 ) err = 3;
3052 if( g->im[i].height <= 0 ) err = 4;
3053 if( g->im[i].hfov <= 0.0 ) err = 5;
3054 if( g->im[i].format == _rectilinear && g->im[i].hfov >= 180.0 ) err = 6;
3055 if( g->im[i].format != _rectilinear && g->im[i].format != _panorama &&
3056 g->im[i].format != _fisheye_circ && g->im[i].format != _fisheye_ff &&
3057 g->im[i].format != _equirectangular && g->im[i].format != _orthographic &&
3058 g->im[i].format != _mirror && g->im[i].format != _stereographic &&
3059 g->im[i].format != _equisolid && g->im[i].format != _thoby)
3060 err = 7;
3061 }
3062
3063
3064 // Check Panorama specs
3065
3066 if( g->pano.hfov <= 0.0 ) err = 5;
3067 if( g->pano.width <=0 ) err = 8;
3068 if( g->pano.height <=0 ) err = 9;
3069 if( g->pano.format == _rectilinear && g->pano.hfov >= 180.0 ) err = 10;
3070
3071
3072 if( g->pano.format != _rectilinear && g->pano.format != _panorama &&
3073 g->pano.format != _equirectangular && g->pano.format != _fisheye_ff &&
3074 g->pano.format != _stereographic && g->pano.format != _mercator &&
3075 g->pano.format != _trans_mercator && g->pano.format != _sinusoidal &&
3076 g->pano.format != _lambert && g->pano.format != _lambertazimuthal &&
3077 g->pano.format != _albersequalareaconic && g->pano.format != _millercylindrical &&
3078 g->pano.format != _panini && g->pano.format != _architectural &&
3079 g->pano.format != _equisolid && g->pano.format != _equipanini &&
3080 g->pano.format != _biplane && g->pano.format != _triplane &&
3081 g->pano.format != _panini_general && g->pano.format != _thoby &&
3082 g->pano.format != _orthographic && g->pano.format != _hammer
3083 ) err=11;
3084
3085 // Check Control Points
3086 for( i=0; i<g->numPts; i++) {
3087 // Joost: cp coordinates can be possible, no problem!
3088 // if( g->cpt[i].x[0] < 0 || g->cpt[i].y[0] < 0 || g->cpt[i].x[1] < 0 || g->cpt[i].y[1] < 0 )
3089 // err = 12;
3090 if( g->cpt[i].num[0] < 0 || g->cpt[i].num[0] >= g->numIm ||
3091 g->cpt[i].num[1] < 0 || g->cpt[i].num[1] >= g->numIm ) err = 13;
3092 }
3093
3094 if( err != -1 ) {
3095 PrintError( errmsg[ err ] );
3096 return -1;
3097 }
3098 else
3099 return 0;
3100 }
3101
3102
CheckMakeParams(aPrefs * aP)3103 static int CheckMakeParams( aPrefs *aP)
3104 {
3105 double im_vfov;
3106 im_vfov = aP->im.hfov / aP->im.width * aP->im.height;
3107
3108 if( (aP->pano.format == _rectilinear) && (aP->pano.hfov >= 180.0) )
3109 {
3110 PrintError("Rectilinear Panorama can not have 180 or more degrees field of view.");
3111 return -1;
3112 }
3113 if( (aP->im.format == _rectilinear) && (aP->im.hfov >= 180.0) )
3114 {
3115 PrintError("Rectilinear Image can not have 180 or more degrees field of view.");
3116 return -1;
3117 }
3118 if( (aP->mode & 7) == _insert ){
3119 if( (aP->im.format == _fisheye_circ || aP->im.format == _fisheye_ff) &&
3120 (aP->im.hfov > MAX_FISHEYE_FOV && im_vfov > MAX_FISHEYE_FOV) ){
3121 PrintError("Fisheye lens processing limited to fov <= %lg", MAX_FISHEYE_FOV);
3122 return -1;
3123 }
3124 }
3125
3126 return 0;
3127
3128 }
3129
3130
3131
3132
3133 // return 0, if overlap exists, else -1
3134 /*
3135 static int GetOverlapRect( PTRect *OvRect, PTRect *r1, PTRect *r2 )
3136 {
3137 OvRect->left = max( r1->left, r2->left );
3138 OvRect->right = min( r1->right, r2->right );
3139 OvRect->top = max( r1->top, r2->top );
3140 OvRect->bottom = min( r1->bottom, r2->bottom );
3141
3142 if( OvRect->right > OvRect->left && OvRect->bottom > OvRect->top )
3143 return 0;
3144 else
3145 return -1;
3146 }
3147 */
3148
3149
SetGlobalPtr(AlignInfo * p)3150 void SetGlobalPtr( AlignInfo *p )
3151 {
3152 optInfo = p;
3153 }
3154
3155
GetGlobalPtr(void)3156 AlignInfo* GetGlobalPtr(void)
3157 {
3158 return optInfo;
3159 }
3160
GetControlPointCoordinates(int i,double * x,double * y,AlignInfo * gl)3161 void GetControlPointCoordinates(int i, double *x, double *y, AlignInfo *gl )
3162 {
3163 double w2, h2;
3164 int j, n[2];
3165
3166 struct MakeParams mp;
3167 struct fDesc stack[15];
3168
3169
3170
3171 n[0] = gl->cpt[i].num[0];
3172 n[1] = gl->cpt[i].num[1];
3173
3174 // Calculate coordinates x/y in panorama
3175
3176 for(j=0; j<2; j++)
3177 {
3178 SetInvMakeParams( stack, &mp, &gl->im[ n[j] ], &gl->pano, 0 );
3179
3180 h2 = (double)gl->im[ n[j] ].height / 2.0 - 0.5;
3181 w2 = (double)gl->im[ n[j] ].width / 2.0 - 0.5;
3182
3183
3184 execute_stack_new( (double)gl->cpt[i].x[j] - w2, // cartesian x-coordinate src
3185 (double)gl->cpt[i].y[j] - h2, // cartesian y-coordinate src
3186 &x[j], &y[j], stack);
3187
3188 h2 = (double)gl->pano.height / 2.0 - 0.5;
3189 w2 = (double)gl->pano.width / 2.0 - 0.5;
3190 x[j] += w2;
3191 y[j] += h2;
3192 }
3193 }
3194
3195
AddEdgePoints(AlignInfo * gl)3196 int AddEdgePoints( AlignInfo *gl )
3197 {
3198 void *tmp;
3199
3200 tmp = realloc( gl->cpt, (gl->numPts+4) * sizeof( controlPoint ) );
3201 if( tmp == NULL ) return -1;
3202 gl->numPts+=4; gl->cpt = (controlPoint*)tmp;
3203
3204 gl->cpt[gl->numPts-4].num[0] = 0;
3205 gl->cpt[gl->numPts-4].num[1] = 1;
3206 gl->cpt[gl->numPts-4].x[0] = gl->cpt[gl->numPts-4].x[1] = -9.0 * (double)gl->pano.width;
3207 gl->cpt[gl->numPts-4].y[0] = gl->cpt[gl->numPts-4].y[1] = -9.0 * (double)gl->pano.height;
3208
3209 gl->cpt[gl->numPts-3].num[0] = 0;
3210 gl->cpt[gl->numPts-3].num[1] = 1;
3211 gl->cpt[gl->numPts-3].x[0] = gl->cpt[gl->numPts-3].x[1] = 10.0 * (double)gl->pano.width;
3212 gl->cpt[gl->numPts-3].y[0] = gl->cpt[gl->numPts-3].y[1] = -9.0 * (double)gl->pano.height;
3213
3214 gl->cpt[gl->numPts-2].num[0] = 0;
3215 gl->cpt[gl->numPts-2].num[1] = 1;
3216 gl->cpt[gl->numPts-2].x[0] = gl->cpt[gl->numPts-2].x[1] = -9.0 * (double)gl->pano.width;
3217 gl->cpt[gl->numPts-2].y[0] = gl->cpt[gl->numPts-2].y[1] = 10.0 * (double)gl->pano.height;
3218
3219 gl->cpt[gl->numPts-1].num[0] = 0;
3220 gl->cpt[gl->numPts-1].num[1] = 1;
3221 gl->cpt[gl->numPts-1].x[0] = gl->cpt[gl->numPts-1].x[1] = 10.0 * (double)gl->pano.width;
3222 gl->cpt[gl->numPts-1].y[0] = gl->cpt[gl->numPts-1].y[1] = 10.0 * (double)gl->pano.height;
3223
3224 return 0;
3225 }
3226
3227
3228