[ptx] Looking for expertise on optimizer internals
alexandre jenny
alexandre.jenny at le-geo.com
Tue Apr 27 06:45:56 BST 2004
Hi Rik,
I'm currently working on my own optimizer which will use quite the
Same algorithm as Panotools which is a optimized version of the
Levenberg marquardt, a non linear least square fitting algorithm.
For details about the panotools lmdif, I think it's a direct conversion
of the lmdif.f using f2c (a fortran to c conversion software), both are
coming from the netlib repository. (You can guess that because there's a
f2c.h
header).
Bye
Alexandre.
> -----Message d'origine-----
> De : ptx-bounces at email-lists.org
> [mailto:ptx-bounces at email-lists.org] De la part de
> Littlefields - Rik, Janis, Kyle & Peter
> Envoyé : mardi 27 avril 2004 06:42
> À : Pablo d'Angelo
> Cc : ptx at email-lists.org
> Objet : Re: [ptx] Looking for expertise on optimizer internals
>
>
> Pablo,
>
> The code is small enough that I will just include it here,
> after my signature.
> There are mods to adjust.c and optimise.c. In adjust.c, I
> replace fcnPano and introduce a function newDistSphere,
> logically replacing the old distSphere. In optimise.C, I
> allocate o->numData*2 function slots instead of o->numData.
>
> Replies to your other comments interspersed below...
>
> Pablo d'Angelo wrote:
>
>
> On Mon, 26 Apr 2004, Littlefields - Rik, Janis, Kyle &
> Peter wrote:
>
>
>
> Briefly, I have working prototype code that
> changes the way that the
> lmdif optimizer is called. The original code
> seeks to optimize distance
> squared, where distance is what lmdif sees. My
> mods also seek to
> optimize distance squared, but the formulation
> is to optimize dx^2 +
> dy^2, where dx and dy are what lmdif sees.
>
> The change in formulation gives lmdif better
> information about how to
> change the parameters, leading to faster and
> more accurate convergence.
>
> But of course there are wrinkles dealing with
> wraparound and near-polar
> points, and I would like to tap into whatever
> other people know about
> the optimizer, to be sure that I don't mess things up..
>
>
>
> I think panotools uses angular differences (for normal
> control points), and
> its probably a not good idea to split the spherical
> coordinates and use them
> like they where normal cartesian coordinates, at least
> near the poles.
>
>
> Yes, it appears to be using latitude/longitude style
> coordinates. These can be handled by scaling longitude by
> sin(latitude) [in 0..pi], then treating as normal cartesian.
> See the code for details.
>
>
> If you have experience with the guts of
> adjust.c and are willing to
> share it, please chime in.
>
>
>
> Sounds very interesting. I'm still thinking about
> adding a (optional)
> penalty for the HFOV optimisation, so that it becomes
> possible to optimise
> it for partial panos as well.
>
>
> I regularly use HFOV optimization for partial panos, but yes,
> there are lots of pits to fall into. We should split this
> topic off to a separate thread.
>
>
>
> Can you post you code, or the patches, so that we can
> see exactly what you
> modified.
>
> ciao
> Pablo
>
>
> Current code follows, in Textpad Compare Files format.
>
> Best regards,
> --Rik
>
> Compare:
> (<)C:\Rik\PanoramaTools2.6.ML12\PanoTools\Sources\adjust.c
> (59753 bytes)
> with:
> (>)C:\Rik\PanoramaTools2.6.ML12_beforeHacking\PanoTools\Source
> s\adjust.c (55785 bytes)
>
> 1246d1246
> < #if 0 // old version, replaced by Rik Littlefield
> 1323,1484d1322 < #endif <
> < /* ---- Revision history ----
> <
> < 04/26/2004, Rik Littlefield, reworked fcnPano and
> newDistSphere to expose dx and dy
> < separately to the optimizer (e.g.lmdif).
> < */
> <
> < // Angular Distance of Control point "num"
> <
> < void newDistSphere (int num, double cdiff[2]) {
> < double x, y ; // Coordinates of control
> point in panorama
> < double w2, h2;
> < int j;
> < Image sph;
> < int n[2];
> < struct MakeParams mp;
> < struct fDesc stack[15];
> < double lat[2], lon[2]; // latitude & longitude, in degrees
> < double dlon;
> <
> < // Get image position in imaginary spherical image
> <
> < SetImageDefaults( &sph );
> <
> < sph.width = 360;
> < sph.height = 180;
> < sph.format = _equirectangular;
> < sph.hfov = 360.0;
> <
> < n[0] = g->cpt[num].num[0];
> < n[1] = g->cpt[num].num[1];
> <
> < // Calculate coordinates x/y in panorama
> <
> < for(j=0; j<2; j++){
> < SetInvMakeParams( stack, &mp, &g->im[ n[j] ], &sph, 0 );
> <
> < h2 = (double)g->im[ n[j] ].height / 2.0 - 0.5;
> < w2 = (double)g->im[ n[j] ].width / 2.0 - 0.5;
> <
> <
> < execute_stack( (double)g->cpt[num].x[j] - w2,
> // cartesian x-coordinate src
> < (double)g->cpt[num].y[j] - h2,
> // cartesian y-coordinate src
> < &x, &y, stack);
> <
> < x = DEG_TO_RAD( x );
> < y = DEG_TO_RAD( y ) + PI/2.0;
> <
> < // x is now in the range -PI to +PI, and y is 0 to PI
> < lat[j] = y;
> < lon[j] = x;
> <
> < }
> <
> < dlon = lon[0]-lon[1];
> < if (dlon <= -PI) dlon += PI;
> < if (dlon >= PI) dlon -= PI;
> < cdiff[0] = (dlon*sin(0.5*(lat[0]+lat[1]))) *
> g->pano.width / ( 2.0 * PI );
> < cdiff[1] = (lat[0]-lat[1]) * g->pano.width / ( 2.0 * PI );
> <
> < return;
> <
> < }
> <
> <
> < // Levenberg-Marquardt function measuring the quality of
> the fit in fvec[] <
> < int fcnPano(m,n,x,fvec,iflag)
> < int m,n;
> < int *iflag;
> < double x[],fvec[];
> < {
> < #pragma unused(n)
> < int i;
> < static int numIt;
> < double result;
> < int iresult;
> < double cdiff[2];
> <
> < if( *iflag == -100 ){ // reset
> < numIt = 0;
> < infoDlg ( _initProgress, "Optimizing Variables" );
> < return 0;
> < }
> < if( *iflag == -99 ){ //
> < infoDlg ( _disposeProgress, "" );
> < return 0;
> < }
> <
> <
> < if( *iflag == 0 )
> < {
> < char message[256];
> < result = 0.0;
> < for( i=0; i < m; i++)
> < {
> < result += fvec[i]*fvec[i] ;
> < }
> < result = sqrt( result/ (double)m );
> <
> < sprintf( message, "Average Difference between
> Controlpoints \nafter %d iteration(s): %25.15g pixels",
> numIt,result);//average);
> < numIt += 10;
> < if( !infoDlg ( _setProgress,message ) )
> < *iflag = -1;
> < return 0;
> < }
> <
> < // Set Parameters
> <
> < SetAlignParams( x ) ;
> <
> < // Calculate distances
> <
> < iresult = 0;
> < for( i=0; i < g->numPts; i++){
> < int j;
> < switch(g->cpt[i].type){
> < // For ordinary control points, return the x-
> and y-errors
> < // separately.
> < case 0: newDistSphere(i,cdiff);
> < fvec[iresult] = cdiff[0];
> < fvec[iresult+1] = cdiff[1];
> < break;
> < // For horizontal, vertical, and line control points,
> < // include two copies of the single distance result.
> < // Note that we return distance, not distance squared.
> < // The optimizer will handle the squaring.
> < case 1:
> < case 2: fvec[iresult] = sqrt(distSquared(i));
> < fvec[iresult+1] = fvec[iresult];
> < break;
> < default:for(j=0; j<g->numPts; j++){
> < if(j!=i && g->cpt[i].type ==
> g->cpt[j].type){
> < fvec[iresult] = sqrt(distLine(i,j));
> < fvec[iresult+1] = fvec[iresult];
> < break;
> < }
> < }
> < break;
> < }
> < iresult += 2;
> < }
> <
> < // If not enough control points are provided, then fill out
> < // the function vector with copies of the average error
> < // for the actual control points. This is really just staving
> < // off disaster, since we have an underdetermined system,
> < // but if we've gotten this far down, it's really about all
> < // we can do.
> <
> < result = 0.0;
> < for (i=0; i < iresult; i++) {
> < result += fvec[i]*fvec[i];
> < }
> < result = sqrt(result/(double)iresult);
> < for (i=iresult; i < m; i++) {
> < fvec[i] = result;
> < }
> <
> < return 0;
> < }
>
> Compare:
> (<)C:\Rik\PanoramaTools2.6.ML12\PanoTools\Sources\optimize.c
> (9861 bytes)
> with:
> (>)C:\Rik\PanoramaTools2.6.ML12_beforeHacking\PanoTools\Source
> s\optimize.c (9414 bytes)
>
> 33,42c33,36
> < // PrintError("RunLMOptimizer");
> < LM.n = o->numVars;
> <
> < // We assume here that each control point contributes
> two functions to be optimized.
> < // For ordinary control points, these functions are the
> X- and Y-coordinate errors.
> < // Horizontal, vertical, and line control points each
> contribute only a single
> < // function logically, but fcnPano (in adjust.c)
> replicates these to consistently
> < // provide two functions per control point.
> <
> < if( o->numData*2 < LM.n )
> ---
> >
> > LM.n = o->numVars;
> >
> > if( o->numData < LM.n )
> 46d40
> < PrintError("You have too few control points or too
> many parameters. Expect strange values!");
> 50c43
> < LM.m = o->numData*2;
> ---
> > LM.m = o->numData;
> 134,137c127,129
> < // PrintError("RunBROptimizer");
> < LM.n = o->numVars;
> <
> < if( o->numData*2 < LM.n )
> ---
> > LM.n = o->numVars;
> >
> > if( o->numData < LM.n )
> 143c135
> < LM.m = o->numData*2;
> ---
> > LM.m = o->numData;
>
>
>
>
More information about the ptX
mailing list