algebraic lens distortion

- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

All Classes Namespaces Files Functions Variables Defines
lens_distortion_estimation.cpp
Go to the documentation of this file.

    00001 /*
    00002  * Copyright 2009, 2010 IPOL Image Processing On Line http://www.ipol.im/
    00003  *
    00004  * This program is free software: you can redistribute it and/or modify
    00005  * it under the terms of the GNU General Public License as published by
    00006  * the Free Software Foundation, either version 3 of the License, or
    00007  * (at your option) any later version.
    00008  *
    00009  * This program is distributed in the hope that it will be useful,
    00010  * but WITHOUT ANY WARRANTY; without even the implied warranty of
    00011  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    00012  * GNU General Public License for more details.
    00013  *
    00014  * You should have received a copy of the GNU General Public License
    00015  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
    00016  */
    00017
    00018
    00098 #include <stdio.h>
    00099 #include <stdlib.h>
    00100 #include <malloc.h>
    00101 #include "image.h"                 /* Functions code associated to read/write images*/
    00102 #include "lens_distortion.h"       /* Functions code associated to amy_prototypes.h */
    00103 #include "point2d.h"
    00104
    00105
    00106
    00107 int main(int argc, char *argv[])
    00108 {
    00109     int width;                /* Image width */
    00110     int height;               /* Image height */
    00111     int Na;                   /* Degree of the lens distortion model polynomial */
    00112     int *Np;                  /* Number of points/line */
    00113     int Nl;                   /* Number of lines */
    00114     int optimize_center;      /* To indicate if the center of distortion must be optimized as well */
    00115     int zoom;                 /* To indicate if the zoom strategy applied */
    00116     double *a;                /* Lens distortion model polynom: it corresponds to the "k"
    00117                               coefficients of the lens distortion model polynom */
    00118     double *solution;         /* Vector having the lens distortion polynom and the center of distortion */
    00119     double **x,**y;           /* Coordinates of points (normalized) */
    00120     double **xx,**yy;         /* Coordinates of points (pixels) */
    00121     double x0,y0;             /* x,y center of distortion of the image (pixels) */
    00122     double *trivial;          /* Vector to save the trivial initial solution, Emin, Vmin, D */
    00123     ami::image<unsigned char> img; /* intermediate image to use image libraries */
    00124
    00125     /* AUXILIAR VARIABLES */
    00126     double factor_n;          /* AUXILIAR VARIABLES TO NORMALIZE COORDINATES */
    00127     int i,cont,m;             /* AUXILIAR VARIABLES */
    00128     char filename[300];       /* POINTER TO THE OUTPUT FILE */
    00129     FILE *fp2;
    00130
    00131     /* WE INIT TO NULL ALL POINTERS */
    00132     a=NULL;
    00133     solution=NULL;
    00134     x=y=xx=yy=NULL;
    00135     Np=NULL;
    00136     trivial=NULL;
    00137
    00138
    00139
    00140     /* WE CHECK COMMAND LINE PARAMETES */
    00141     /* WE HAVE SEVERAL OPTIONS */
    00142
    00143     /* FIRST OPTION: NUMBER OF INPUT PARAMETERS IS 5
    00144         "lens_distortion_estimation input_image.bmp output_undistorted_image.bmp
    00145          input_line_primitives.dat output_lens_distortion_models.dat"
    00146         - CENTER OF DISTORTION: CENTER OF THE IMAGE
    00147          (IT IS READ INTERNALLY FROM THE LOADED IMAGE)
    00148         - THE CENTER OF DISTORTION IS NOT OPTIMIZED
    00149     */
    00150
    00151     /* SECOND OPTION: NUMBER OF INPUT PARAMETERS IS 6
    00152         "lens_distortion_estimation input_image.bmp output_undistorted_image.bmp
    00153        input_line_primitives.dat output_lens_distortion_models.dat optimize_center"
    00154        - CENTER OF DISTORTION: CENTER OF THE IMAGE (IT IS READ INTERNALLY FROM THE LOADED IMAGE)
    00155        - THE CENTER OF DISTORTION IS OPTIMIZED
    00156     */
    00157
    00158     /* THIRD OPTION: NUMBER OF INPUT PARAMETERS IS 7
    00159         "lens_distortion_estimation input_image.bmp output_undistorted_image.bmp
    00160        input_line_primitives.dat output_lens_distortion_models.dat x_center y_center"
    00161        - CENTER OF DISTORTION: IT IS INDICATED BY THE USER
    00162        - THE CENTER OF DISTORTION IS NOT OPTIMIZED
    00163     */
    00164
    00165     /* FOURTH OPTION: NUMBER OF INPUT PARAMETERS IS 8
    00166         "lens_distortion_estimation input_image.bmp output_undistorted_image.bmp
    00167        input_line_primitives.dat output_lens_distortion_models.dat x_center y_center optimize_center"
    00168        - CENTER OF DISTORTION: IT IS INDICATED BY THE USER
    00169        - THE CENTER OF DISTORTION IS OPTIMIZED
    00170     */
    00171
    00172     /*
    00173       input_image.bmp:                    image to correct the radial distortion
    00174       output_undistorted_image.bmp        output corrected image
    00175       input_line_primitives.dat           file (ASCII) with the primitives of the image (points & lines)
    00176       output_lens_distortion_models.dat   output file with the lens distortion coefficients and energy values
    00177       center_image                        1 to select the center of the image as center of distortion
    00178       x_center                            if center_image==0, the x_center is taken as x_center of distortion
    00179       y_center                            if center_image==0, the y_center is taken as y_center of distortion
    00180       optimize_center                     1: to indicate to optimize the center of distortion
    00181                                           0: to indicate not to optimize the center of distortion
    00182     */
    00183
    00184
    00185     // WE PRINT ON SCREEN THAT THE PROGRAM STARTS WORKING
    00186     printf("\n****************************************************************************");
    00187     printf("\nLENS DISTORTION PROGRAM STARTS WORKING\n");
    00188
    00189     /* WE CAPTURE THE MULTIPLE CHOICE */
    00190     if(argc==5) {   /* FIRST OPTION */
    00191         /* WE READ INPUT IMAGE FROM DISK */
    00192         strcpy(filename,argv[1]);
    00193         if(img.read(filename)!=0)
    00194         {
    00195             printf("Image can not be loaded\n");
    00196             return(-1);
    00197         }
    00198         width=img.width();
    00199         height=img.height();
    00200         /* WE READ POINT LINE PRIMITES FROM DISK */
    00201         if (read_line_primitives(argv[3],&Nl,&Np,&x,&y)!=0) {
    00202             printf("point line primitives can not be loaded\n");
    00203             return(-1);
    00204         }
    00205         x0=(double)width/2.0;
    00206         y0=(double)height/2.0;
    00207         optimize_center=0;
    00208     }
    00209
    00210     else if(argc==6) {   /* SECOND OPTION */
    00211         /* WE READ INPUT IMAGE FROM DISK */
    00212         strcpy(filename,argv[1]);
    00213         if(img.read(filename)!=0)
    00214         {
    00215             printf("Image can not be loaded\n");
    00216             return(-1);
    00217         }
    00218         width=img.width();
    00219         height=img.height();
    00220         /* WE READ POINT LINE PRIMITES FROM DISK */
    00221         if (read_line_primitives(argv[3],&Nl,&Np,&x,&y)!=0) {
    00222             printf("point line primitives can not be loaded\n");
    00223             return(-1);
    00224         }
    00225         x0=(double)width/2.0;
    00226         y0=(double)height/2.0;
    00227         /* WE CAPTURE THE OPTION TO OPTIMIZE THE CENTER OF DISTORTION */
    00228         optimize_center=atoi(argv[5]);
    00229     }
    00230
    00231     else if(argc==7) {   /* THIRD OPTION */
    00232         /* WE READ INPUT IMAGE FROM DISK */
    00233         strcpy(filename,argv[1]);
    00234         if(img.read(filename)!=0)
    00235         {
    00236             printf("Image can not be loaded\n");
    00237             return(-1);
    00238         }
    00239         width=img.width();
    00240         height=img.height();
    00241         /* WE READ POINT LINE PRIMITES FROM DISK */
    00242         if (read_line_primitives(argv[3],&Nl,&Np,&x,&y)!=0) {
    00243             printf("point line primitives can not be loaded\n");
    00244             return(-1);
    00245         }
    00246         /* WE CAPTURE THE CENTER OF DISTORTION FROM THE GRAPHIC INTERFACE */
    00247         x0=atof(argv[5]);
    00248         y0=atof(argv[6]);
    00249         optimize_center=0;
    00250     }
    00251
    00252     else if(argc==8) {   /* FOURTH OPTION */
    00253         /* WE READ INPUT IMAGE FROM DISK */
    00254         strcpy(filename,argv[1]);
    00255         if(img.read(filename)!=0)
    00256         {
    00257             printf("Image can not be loaded\n");
    00258             return(-1);
    00259         }
    00260         width=img.width();
    00261         height=img.height();
    00262         /* WE READ POINT LINE PRIMITES FROM DISK */
    00263         if (read_line_primitives(argv[3],&Nl,&Np,&x,&y)!=0) {
    00264             printf("point line primitives can not be loaded\n");
    00265             return(-1);
    00266         }
    00267         /* WE CAPTURE THE CENTER OF DISTORTION FROM THE GRAPHIC INTERFACE */
    00268         x0=atof(argv[5]);
    00269         y0=atof(argv[6]);
    00270         /* WE CAPTURE THE OPTION TO OPTIMIZE THE CENTER OF DISTORTION */
    00271         optimize_center=atoi(argv[7]);
    00272     }
    00273
    00274     else {
    00275         printf("lens_distortion_estimation : number of parameter inadequate\n");
    00276         printf("Usage\n");
    00277         printf("lens_distortion_estimation input_image.bmp output_undistorted_image.bmp input_line_primitives.dat output_lens_distortion_models.dat  [x_center y_center] [optimize_center]");
    00278         return(-1);
    00279     }
    00280
    00281     /* WE CHECK THAT THE CENTER OF DISTORTION IS CONSISTENT WITH THE IMAGE */
    00282     if(x0<0.0) {
    00283         printf("\n The x-coordinate for the center of distortion is < 0.");
    00284         printf("\n The x-coordinate should be 0.0 <x < width of the image.");
    00285         return(-1);
    00286     }
    00287     if(x0>width) {
    00288         printf("\n The x-coordinate for the center of distortion is > width of the image.");
    00289         printf("\n The x-coordinate should be 0.0 < x <width of the image.");
    00290         return(-1);
    00291     }
    00292     if(y0<0.0) {
    00293         printf("\n The y-coordinate for the center of distortion is < 0.");
    00294         printf("\n The y-coordinate should be 0.0 <x < height of the image.");
    00295         return(-1);
    00296     }
    00297     if(y0>height) {
    00298         printf("\n The y-coordinate for the center of distortion is > width of the image.");
    00299         printf("\n The y-coordinate should be 0.0 < x < height of the image.");
    00300         return(-1);
    00301     }
    00302
    00303
    00304
    00305     /* WE CHECK THAT THE SELECTED VALUE FOR THE VARIABLE TO ACCOUNT FOR OPTIMIZING
    00306        THE CENTER OF DISTORTION IS CONSISTENT
    00307           optimize_center=0: do not optimize the center of distortion (it is fixed
    00308                            to the indicated (x-center,y_center)
    00309           optimize_center=1: optimize the center of distortion using the patch search
    00310           strategy -pixel precision- and then, it is optimized by the gradient from the algebraic solution
    00311           at sub-pixel precision.
    00312      */
    00313     if(optimize_center!=0 && optimize_center!=1) {
    00314         printf("\n The selection for this input variable is not allowed.");
    00315         printf("\n The selection should be 0 (not optimize); or 1 (optimize).");
    00316         return(-1);
    00317     }
    00318
    00319     /* TO AVOID A BAD SOLUTION WHEN OPTIMIZING THE CENTER OF DISTORTION USING THE A MINIMUM
    00320     NUMBER OF PRIMITIVES, IT IS NOT ALLOWED OPTIMIZING THE CENTER OF DISTORTION FOR THE CASE OF
    00321     USING ONLY ONE LINE. FROM THE EXPERIENCE, WE NOTE THAT THE ALGORITHM WHICH OPTIMIZES THE
    00322     CENTER OF DISTORTION CAN BE TRAPPED INTO A BAD LOCAL MINIMA. AS A RESULT, THE SOLUTION CAN BE
    00323     A BAD ONE, AND HENCE, THIS SITUATION IS SIMPLY AVOIDED */
    00324     if(Nl==1) optimize_center=0;
    00325
    00326     /* WE ALLOCATE MEMORY FOR AUXILARY POINTS AND WE NORMALIZE ORIGINAL POINTS */
    00327     xx=(double**)malloc( sizeof(double*)*Nl);  /*  xx pixels coordinates */
    00328     yy=(double**)malloc( sizeof(double*)*Nl);  /*  yy pixels coordinates */
    00329     for(i=0; i<Nl; i++) {
    00330         xx[i]=(double*)malloc( sizeof(double)*Np[i] );
    00331         yy[i]=(double*)malloc( sizeof(double)*Np[i] );
    00332     }
    00333     printf("****************************************************************************");
    00334     printf("\nInitial Center of distortion:");
    00335     printf("\nx0=%f,y0=%f",x0,y0);
    00336     if(optimize_center==1) printf("\nIs the Center of Distortion going to be optimized?: YES");
    00337     else printf("\nIs the Center of Distortion going to be optimized?: NO");
    00338     printf("\n****************************************************************************");
    00339
    00340     /* INITIALIZE MEMORY FOR IMAGE */
    00341     int size=width*height;
    00342     int size2=2*size;
    00343
    00344     /* WE KEEP A COPY OF THE PIXEL COORDINATES VECTORS */
    00345     for(m=0; m<Nl; m++) {
    00346         for(i=0; i<Np[m]; i++) {
    00347             xx[m][i]=x[m][i];
    00348             yy[m][i]=y[m][i];
    00349         }
    00350     }
    00351
    00352
    00353     /* WE ALLOCATE MEMORY FOR MAXIMUM LENS DISTORTION POLYNOMIAL SIZE AND WE INITIALIZE TO TRIVIAL CASE */
    00354     /* WE USE IT ONLY FOR THE ALGEBRAIC ALONE METHOD */
    00355     ami_calloc1d(a,double,7);
    00356     a[0]=1.0;
    00357     /* MAXIMUM DEGREE OF THE LENS DISTORTION POLYNOM */
    00358     Na=4;
    00359
    00360     /* WE COPY THE CENTER OF DISTORTION TO THE THE 5 AND 6th POSITIONS*/
    00361     a[5]=x0;
    00362     a[6]=y0;
    00363
    00364     /* INITIAL TRIVIAL SOLUTION IS IN VECTOR a=[1,0,0,0,0,x0,y0] */
    00365
    00366     /* WE COMPUTE THE SOLUTIONS AND SAVE IT TO THE INDICATED "output_lens_distortion_models.dat"
    00367     WE PREPARE THE OUTPUT FILE */
    00368     if(!(fp2=fopen(argv[4],"w"))) {
    00369         printf("Error while opening the %s file",filename);
    00370         exit(0);
    00371     }
    00372
    00373     /* FIRST WE COMPUTE THE TRIVIAL SOLUTION AND SAVE IT TO THE OUTPUT FILE */
    00374     /* WE ALLOCATE MEMORY TRIVIAL SOLUTION FOR Emin, Vmin, D */
    00375     ami_calloc1d(trivial,double,3);
    00376     trivial_solution(Nl,Np,a,xx,yy,factor_n,fp2,trivial,optimize_center);
    00377
    00378
    00379    /* SELECT THE BIGGEST VALUE (FOR THE INPUT IMAGE): THE POINT IN A CORNER, TO BE USED IN
    00380     THE LOCATION OF THE ALGEBRAIC ROOTS */
    00381     double xtmp=width-width/2.0;
    00382     double ytmp=height-height/2.0;
    00383     double max_radius=sqrt(pow(xtmp,2.0)+pow(ytmp,2.0));
    00384
    00385     /* OPTIMIZED CENTER: WE SCAN A PATCH CENTERED AT THE GIVEN CENTER OF DISTORTION:
    00386        WE SELECT THE BEST SOLUTION AND THEN APPLIED THE GRADIENT METHOD,
    00387        TO GET THE SOLUTIONS USING ZOOM AND NOT USING ZOOM
    00388       */
    00389     if(optimize_center==1) {
    00390         search_for_best_center(Nl,Np,a,xx,yy,width,height,max_radius);
    00391         x0=a[5];      /* WE CAPTURE THE OPTIMIZED CENTER OF DISTORTION */
    00392         y0=a[6];
    00393     }
    00394
    00395     /* ALGEBRAIC METHOD & GRADIENT METHOD WORK WITH NORMALIZED UNITS */
    00396     factor_n=0.0;
    00397     cont=0;
    00398     for(m=0; m<Nl; m++) {
    00399         for(i=0; i<Np[m]; i++) {
    00400             x[m][i]-=x0;
    00401             y[m][i]-=y0;
    00402             cont++;
    00403             factor_n+=x[m][i]*x[m][i]+y[m][i]*y[m][i];
    00404         }
    00405     }
    00406     factor_n=sqrt(factor_n/cont);
    00407     for(m=0; m<Nl; m++) for(i=0; i<Np[m]; i++) {
    00408             x[m][i]/=factor_n;
    00409             y[m][i]/=factor_n;
    00410     }
    00411
    00412  /* IN NORMALIZED COORDINATES */
    00413    xtmp=xtmp/factor_n;ytmp=ytmp/factor_n;
    00414    max_radius=sqrt(pow(xtmp,2.0)+pow(ytmp,2.0));
    00415
    00416    printf("\nmax_radius=%f",max_radius);
    00417    printf("\n****************************************************************************");
    00418
    00419
    00420     /* ALGEBRAIC METHOD FROM TRIVIAL SOLUTION + GRADIENT TO IMPROVE THE ALGEBRAIC SOLUTION; WITHOUT ZOOM */
    00421     /* ALGEBRAIC METHOD & GRADIENT METHOD WORK WITH NORMALIZED UNITS */
    00422     zoom=0;       /* NO ZOOM APPLIED */
    00423     algebraic_method_pre_gradient(Nl,Np,a,x,y,xx,yy,factor_n,zoom,fp2,optimize_center,max_radius);
    00424
    00425     /* ALGEBRAIC METHOD FROM TRIVIAL SOLUTION + GRADIENT TO IMPROVE THE ALGEBRAIC SOLUTION; WITH ZOOM */
    00426     zoom=1;       /* ZOOM APPLIED */
    00427     algebraic_method_pre_gradient(Nl,Np,a,x,y,xx,yy,factor_n,zoom,fp2,optimize_center,max_radius);
    00428     fprintf(fp2,"\n*******************************END OF FILE*******************************");
    00429     fclose(fp2);  /* CLOSE THE OUTPUT FILE */
    00430
    00431     /* PRINT SOME INFORMATION ON SCREEN */
    00432     printf("\nFinal Center of Distortion:");
    00433     printf("\nx0=%f,y0=%f",a[5],a[6]);
    00434     printf("\n****************************************************************************");
    00435     printf("\nSolution:\n");
    00436     for(i=0; i<=Na; i++) printf("k[%d]=%e\n",i,a[i]);
    00437     printf("****************************************************************************");
    00438     /* WE COMPUTE THE UNDISTORTED IMAGE */
    00439     ami::point2d<double> dc(a[5],a[6]);
    00440
    00441     /* FAST UNDISTORT INVERSE */
    00442     cout << endl << "Undistorting image..." << endl;
    00443     ami::image<unsigned char> img_out3 = undistort_image_inverse_fast(img,Na,a,dc,1.);
    00444     strcpy(filename,argv[2]);
    00445     if(img_out3.write(filename)!=0) printf("Image can not be saved\n");
    00446
    00447     /* WE FREE THE MEMORY */
    00448     if(a!=NULL) free(a);
    00449     if(x!=NULL) {
    00450         for(i=0; i<Nl; i++) {
    00451             if(x[i]!=NULL) free(x[i]);
    00452         }
    00453         free(x);
    00454     }
    00455     if(y!=NULL) {
    00456         for(i=0; i<Nl; i++) {
    00457             if(y[i]!=NULL) free(y[i]);
    00458         }
    00459         free(y);
    00460     }
    00461     if(xx!=NULL) {
    00462         for(i=0; i<Nl; i++) {
    00463             if(xx[i]!=NULL) free(xx[i]);
    00464         }
    00465         free(xx);
    00466     }
    00467     if(yy!=NULL) {
    00468         for(i=0; i<Nl; i++) {
    00469             if(yy[i]!=NULL) free(yy[i]);
    00470         }
    00471         free(yy);
    00472     }
    00473     if(trivial!=NULL) free(trivial);
    00474     if(Np!=NULL) free(Np);
    00475
    00476     // WE PRINT ON SCREEN THAT THE PROGRAM FINISHED WORKING
    00477     printf("\n****************************************************************************");
    00478     printf("\nLENS DISTORTION PROGRAM FINISHED WORKING\n");
    00479     return(0);
    00480
    00481 }
    00482