[an error occurred while processing this directive]
Main Page | Modules | Namespace List | Class Hierarchy | Class List | File List | Namespace Members | Class Members | File Members | Related Pages

lens_calibrate/lensFunc.h

Go to the documentation of this file.
00001 /* lensFunc.h           for calibrate_lens, 28July09 TKS
00002 
00003   This is a C++ header, not to be used in libpano source.
00004   In that context, use "lensFunc_glue.h"
00005 
00006   lensFunc encapsulates calibrated lens projection functions 
00007   in a form compatible with the PanoTools remapping function 
00008   stack and parameter optimizer.
00009 
00010   By providing what is essentially a single remapping function 
00011   for all lenses, and only for lenses, lensFunc can support a 
00012   more comprehensive set of lens models and a practical lens 
00013   calibration facility, both of which are sorely needed.  
00014 
00015   Remapping with lensFunc is also faster than with existing
00016   stack functions because it tabulates the radial part of the
00017   mapping as a cubic spine.  There are no evaluations of trig 
00018   functions or a radial polynomial during remapping, and
00019   forward and inverse mappings are equally fast.
00020 
00021   There are also functions to map image coordinates to/from
00022   points on the unit Cartesian sphere. The projection center 
00023   is the center of the sphere, and the image center is tangent 
00024   to the sphere at (0,0,1).  The sphere coordinates are right
00025   handed if image x runs right and y runs up -- caution: in
00026   most image formats, y runs down, which gives left handed
00027   3D coordinates.  These functions use the radial spline, but 
00028   do have to evaluate sin & cos, or arccos, as well.
00029  
00030   The C interface to libpano is very thin: two new remapping
00031   stack functions; two to copy lensFunc parameters to/from 
00032   an array of optimizer parameters, and one to select which set 
00033   of parameters is copied; one new format code; and a new void * 
00034   parameter table entry as a handle to an instance of lensFunc.
00035 
00036   The lensFunc parameters are expected to come from calibration
00037   procedures, via new API's separate from those used for the 
00038   traditional PT corrections.  It is intended that those APIs
00039   be developed, not inside libpano, but in the C++ contexts of 
00040   Hugin and other programs that use libpano.  
00041   
00042   The API for parameter optimization is compatible with libpano's
00043   but a bit more flexible.  The lensFunc basically decides what 
00044   parameters to make available to the optimizer, based on the 
00045   declared purpose of the optimization (lens calibration, 
00046   pano alignment, ...).  Two functions then copy that set of 
00047   values to and from an external double array during optimization.  
00048 
00049   All lens parameters and the center shifts can be adjusted via the 
00050   optimizer API.  Calibrated parameters are protected during routine 
00051   panorama alignment; however apparent fov and distortion can be 
00052   optimized via secondary parameters.
00053 
00054   During optimization, certain parameters need to be kept within
00055   limited ranges.  There is a flavor of levmar that can do that, 
00056   however it runs less than half as fast as the unconstrained
00057   version.  So lensFunc presents unrestricted versions of those 
00058   parameters to the optimizer, using the tangent function to map
00059   the limited range to [-infinity:infinity], and you can call the
00060   normal levmar-diff().
00061    
00062   The traditional radius correction polynomial causes all kinds
00063   of trouble for optimization, so lensFunc does not use one.  
00064   Instead it uses an adjustable spline, that is guaranteed to be 
00065   monotonic no matter what parameter values the optimizer chooses.
00066   The number of radius correction points can be 0, 1, 3, 7, or 15; 
00067   each is one optimizer parameter. Unlike polynomial coefficients 
00068   these parameters are stable under optimization.  Using more of 
00069   them improves detail but does not encourage wild excursions.
00070 
00071   The radius function is normalized to the maximum possible pixel
00072   radius, which is half the diagonal (PT's constant radius is half 
00073   the larger image dimension).
00074 
00075   lensFunc implements five "ideal" lens models: rectilinear, 
00076   stereographic, equal-angle, equal-area, orthographic; and two
00077   "generic" models: fisheye and universal.  
00078   
00079   The ideal models are not adjustable, and are only valid for 
00080   certain combinations of focal length, image diagonal, and field
00081   of view.  The generic models have adjustable scale limits, that 
00082   are set such that no illegal trig trig function arguments can be 
00083   generated. Those limits can change due to optimization of focal 
00084   length scale, otherwise they are fixed by the design parameters 
00085   given to setCamLens().
00086 
00087   The only optional setup parameter, the crop radius in pixels,
00088   defaults to half the diagonal of the image.  It needs to be
00089   specified when the true image circle is smaller than the 
00090   sensor diagonal, to ensure that the full angular field of the 
00091   lens can be mapped.  CAUTION: image coordinates outside the
00092   crop circle, or sphere coordinates outside the corresponding
00093   view cone, will likely cause numerical errors and must never
00094   be passed in for remapping.  Note in case of non-square
00095   pixels, give crop radius in vertical pixels.
00096 
00097   The generic fisheye model is a weighted average of two fixed 
00098   functions, nominally the stereographic and equal-area projections, 
00099   but with the angle scale set initially according to the ratio of 
00100   image diagonal to focal length.  Its parameter [0:1] adjusts the 
00101   fraction of the stereographic function in the mixture.
00102 
00103   The universal model is the same as the generic fisheye model
00104   but with the angle scale factor also adjustable.  It can thus
00105   approximate the rectilinear, orthographic and equal-angle as well 
00106   as the generic fisheye projections.  
00107   
00108   Other features of the calibration API: 
00109 
00110   Physical parameters are kept in standard physical units: mm for 
00111   lengths, degrees for angles.
00112 
00113   Wherever possible parameters are defined so that the value zero
00114   is either a reasonable default, or means "ignore this parameter".
00115 
00116   Lens and camera parameters are specified separately.  The basic
00117   lens parameters are focal length (mm), nominal projection
00118   function, and radial correction factors including TCA.  Basic 
00119   camera parameters are sensor size and resolution (pixels per mm, 
00120   each way) and center point shifts (mm).  This separation reflects 
00121   the reality that one lens might be used on several cameras, or that 
00122   an image might have been digitized from film.  We also cater for 
00123   cameras having non-square pixels, and slit-cameras having a lens 
00124   projection only on the vertical axis.
00125 
00126   Following the libpano termninology, we call the projection
00127   in which coordinates are given the "destination" and the one
00128   whose coordinates are to be computed the "source"; and name
00129   transformation functions <destination>_<source>.  But function 
00130   arguments are named after the projection they belong to.
00131    
00132 
00133 Notes on the Derivation of a Universal Lens Model
00134 
00135         Photographic lenses have projection functions that lie between
00136         r = tan( a ) (rectilinear) and r = sin( a ) (orthographic). See
00137         http://www.bobatkins.com/photography/technical/field_of_view.html
00138 
00139         So we should be able to compose a universal lens model out of
00140         those two functions.  The simplest such model is a weighted
00141         average: r = p * tan(a) + (1-p) * sin(a).  But this is limited
00142         to a < 90 degrees, while some lenses have fov > 180 degrees.  To
00143         increase the fov we can scale the angle by a factor 1/k:
00144           r = p * k * tan(a/k) + (1-p) * k * sin(a/k)
00145     Multiplying the radius by k normalizes the slope to 1 at a = 0
00146         which effectively holds the focal length constant.  When k = 2, 
00147         this is a weighted average of the stereographic and     azimuthal 
00148         equal-area projections, which is fine for most fish eye lenses.  
00149         But k must be 1 to fit rectilinear and orthographic lenses; and
00150         for the equal-angle projection, k must approach infinity.  As k 
00151         grows large, p loses effect, because the initial parts of the 
00152         sin and tan functions are very similar. 
00153 
00154         For fitting image data one uses the inverse model:
00155           a = p * k * atan(r/k) + (1-p) * k * asin(r/k)
00156         In the sine function, k must not be allowed to get smaller
00157         than the actual maximum radius, as asin(r/k) is undefined 
00158         when r > k.  So maybe we need two different k values:
00159           a = p * kt * atan(r/kt) + (1-p) * ks * asin(r/ks)
00160     with ks held >= rmax.  Does this mean we need 3 model 
00161         parameters?
00162 
00163         The answer has to be no.  Each of p, ks, kt has a "dead zone"
00164         where is has no effect: ks is dead when p = 1, kt is dead when
00165         p = 0, and p is dead when ks and kt are both large.  So they 
00166         are rather like the triangular coordinates of a 2D point.
00167         We can control them with 2 parameters, u and v, that range from 
00168         0 to 1:  p = 1/2 + u * (v - 1/2), ks = R/u, kt = 1/u.  Of course
00169         the lower limit on u must be just a hair positive to keep ks 
00170         and kt finite.
00171           
00172         A big problem is how to estimate the actual field of view.
00173         Let A be the maximum angle of incidence = fov/2, and R be the
00174         corresponding maximum radius in the focal plane, divided by the
00175         focal length.  There is clearly an upper limit on R, namely 
00176         half the diagonal of the image sensor / FL. The model gives
00177           A = p * k * atan( R/k ) + (1-p) * k * asin( R/k )
00178         But we can't compute A until we have fixed p and k by fitting
00179         the model to some data.  Moreover, the actual image circle may
00180         be smaller than the sensor diagonal.
00181 
00182 */
00183 
00184 /***************************************************************************
00185  *   Copyright (C) 2009 Thomas K Sharpless                                 *
00186  *   tksharpless@gmail.com                                                 *
00187  *                                                                         *
00188  *   This program is free software; you can redistribute it and/or modify  *
00189  *   it under the terms of the GNU General Public License as published by  *
00190  *   the Free Software Foundation; either version 2 of the License, or     *
00191  *   (at your option) any later version.                                   *
00192  *                                                                         *
00193  *   This program is distributed in the hope that it will be useful,       *
00194  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00195  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00196  *   GNU General Public License for more details.                          *
00197  *                                                                         *
00198  *   You should have received a copy of the GNU General Public License     *
00199  *   along with this program; if not, write to the                         *
00200  *   Free Software Foundation, Inc.,                                       *
00201  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
00202  ***************************************************************************/
00203 
00204 #ifndef _LENSFUNC_H
00205 #define _LENSFUNC_H
00206 
00207 /* ssifn -- 
00208   A simple strictly increasing adjustable function.
00209 
00210   Spans the range [0:1] in both x and y and can be
00211   read as y(x) or as x(y).
00212 
00213   The fn is a cubic spline curve, defined recursively in 
00214   terms of the number of levels, k. There are 2^^k 
00215   segments and 2^^k - 1 parameters. For k = 1 there is 1 
00216   parameter, that locates a point on the diagonal from 
00217   (0,1) to (1,0).  That point is the mutual corner of two 
00218   new boxes on whose diagonals two level 2 parameters
00219   locate new points.  That creates 4 boxes at level 3, and 
00220   so on. 
00221 
00222   The corner points are tabulated as x and y arrays and
00223   converted to a pair of monotonic Hermite splines that
00224   can be read as y(x) and x(y).
00225 
00226   The basic parameters range from -1/2 at the lower right
00227   corner of a box to 1/2 at the upper left corner.  For
00228   convenience in dealing with unconstrained optimizers,
00229   there is an interface that maps this range to the 
00230   tangent of an angle (-infinity at lwr rgt, +infinity
00231   at upper left).
00232 
00233   Initially all parameters are 0 (straight line).
00234 
00235   The fns to copy params in and out take an optional argument
00236   k to select a "leading subset", e.g. k = 1 just copies the
00237   center point parameter.  This k must be <= the k passed to
00238   c'tor and defaults to that value.  Note the params are
00239   copied in left to right order in any case.
00240 
00241 */
00242 
00243 class ssifn {
00244 public:
00245         ssifn( int k = 4 );     // default 15 adjustable points
00246         ~ssifn();
00247         int getParamCount( int k = 0 ); // returns number of params
00248         int getParams( double * p, int k = 0 ); // returns count
00249         int setParams( double * p, int k = 0 ); // returns 1 ok 0 fail
00250         int getInfParams( double * p, int k = 0 );      // returns count
00251         int setInfParams( double * p, int k = 0 );      // returns 1 ok 0 fail
00252         double yofx( double x );
00253         double xofy( double y );
00254 protected:
00255         int kmax;
00256         int npts;
00257         double * para;
00258         double * xval, * yval;
00259         double * Cxy, * Cyx;
00260         int setSpline();
00261         void sssub(int l, int r);
00262 };
00263 
00264 /* A libpano stack function transforms 'dest' coordinates to 'source' 
00265    coordinates according to an unspecified set of parameters.
00266    These definitions are from filter.h
00267 */
00268 typedef int (*trfn)( double x_dest,double  y_dest, double* x_src, double* y_src, void* params );
00269 // Function descriptor to be executed by exec_function
00270 struct fDesc {
00271         trfn    func;                   // The function to be called
00272         void    *param;                 // The parameters to be used
00273 };              
00274 
00275 class lensFunc {
00276 public:
00277         typedef enum{
00278                 no_lens = 0,
00279                 rectilinear,
00280                 equal_angle,
00281                 equal_area,
00282                 stereographic,
00283                 orthographic,
00284                 fisheye_model,
00285                 universal_model
00286         } lens_kind;
00287 
00288         typedef enum{
00289                 no_cam = 0,
00290                 digicam,
00291                 slitcam,
00292                 scanner
00293         } camera_kind;
00294 
00295         typedef enum{
00296                 align = 0,
00297                 calibrate
00298         } opt_purpose;
00299 
00300         lensFunc();
00301         lensFunc( camera_kind c_type,
00302                           int hPixels, int vPixels,             // sensor dimensions
00303                       double hPixmm, double vPixmm, // pixel spacings
00304                           lens_kind l_type, double FLmm,
00305                           double cropRadPix = 0
00306                         );
00307 
00308         ~lensFunc(){}
00309 
00310 /* coordinate mapping API
00311    CAUTION: all image coordinates must be in a circle whose diameter
00312    is the diagonal of the sensor given to the c'tor or setCamLens().  
00313    Points outside that circle may cause errors.
00314 */
00315 
00316   /* put function and parameter pointers in a libPano fDesc
00317      Note these remain valid so long as the lens model is
00318          not changed via the calibration API  */
00319         int fD_photo_sphere( fDesc * pfd );     // to real 
00320         int fD_sphere_photo( fDesc * pfd );     // to ideal
00321   /* compute the image <=> 2d sphere (equal angle fisheye) mappings */
00322         int photo_sphere( double x_sphr, double  y_sphr, double* x_phot, double* y_phot ); 
00323         int sphere_photo( double x_phot, double  y_phot, double* x_sphr, double* y_sphr );
00324   /* compute image <=> Cartesian unit sphere mappings */
00325         int photo_cart3d( double xyz_sphr[3], double* x_phot, double* y_phot ); 
00326         int cart3d_photo( double x_phot, double  y_phot, double xyz_sphr[3] );
00327 /* set the color channel to be mapped (selects a tca parameter) */
00328         bool setColorChannel( int rgb );
00329   /* parameter optimization API
00330     setOptPurpose() selects a set of parameters and returns their number
00331           polyN:  0 -- no radial polynomial
00332                   1 -- default for lens model
00333                           2 -- even-order
00334                           3 -- cubic
00335                           4 -- quintic
00336           optCen: true to optimize center offsets
00337           optScl: true to optimize h and v scale factors
00338           optTca: true to optimize tca factor for color
00339     get/setOptParams() copy the selected values out/in
00340   */
00341         int setOptPurpose( opt_purpose p,
00342                                int optRad = 2,  // detail level of radius correction
00343                                            bool optCen = false, // optimize center shifts
00344                                            bool optScl = false, // optimize fov
00345                                            bool optTca = false );       // optimize a tca param
00346         int getOptParams( double * ppa );       // copy to array
00347         int setOptParams( double * ppa );       // copy from array
00348 
00349 /* Calibration API
00350 
00351     setCamLens() must be called first (typically from c'tor).  It
00352         returns a pointer to a text error message, or null if no error.
00353         
00354         Then you can either optimize the model, or reload optimized primary 
00355         parameters.  Don't load "made up" parameters.  
00356         
00357         After an optimization you can read out the working optimized values 
00358         and store them as primary if OK.
00359 
00360         Each "set" step yields a usable mapping function unless the
00361         passed parameter values are illegal, when these fns return false
00362         with all old values retained.
00363 
00364   */
00365         const char *
00366                 setCamLens(  camera_kind c_type,
00367                                           int hPixels, int vPixels,
00368                                           double hPixmm, double vPixmm,
00369                                           lens_kind l_type, double FLmm,
00370                                           double cropRadPix = 0
00371                                         );
00372 
00373   // inquiry functions
00374         double getFL_mm(){ return FL_mm; }
00375         double avgFL_pix(){ return 0.5 * (hFLpix + vFLpix); }
00376         const char * errMsg(){ return errmsg; } // null if no error
00377         const char * camDesc(); 
00378         const char * lensDesc();
00379 
00380   // translate internal radius function to/from external polynomial
00381 //      bool getRadiusPoly( double * pc, int polyN, bool fwd = false );
00382 //      bool setRadiusPoly( double * pc, int polyN, bool fwd = false );
00383 
00384   // set primary and working parameter values
00385         bool set_f_scale( double h_scl, double v_scl );
00386         bool set_c_shift( double h_shfmm, double y_shfmm );
00387   // use set_proj_params() only to reload optimized values!
00388         bool set_proj_params( double sina, double tana, double tanf );
00389         bool set_rad_params( double pp[7] );
00390         bool set_tca_params( double RperG, double BperG );
00391   // get primary parameter values
00392         void get_f_scale( double& h_scl, double& v_scl );
00393         void get_c_shift( double& h_shfmm, double& v_shfmm );
00394         void get_proj_params( double & sina, double & tana, double & tanf );
00395         void get_rad_params( double pp[7] );
00396         void get_tca_params( double& RperG, double& BperG );
00397   // get working (optimizable) parameter vlaues
00398         void get_w_f_scale( double& h_scl, double& v_scl );
00399         void get_w_c_shift( double& h_shfmm, double& v_shfmm );
00400         void get_w_proj_params( double & sina, double & tana, double & tanf );
00401         void get_w_rad_params( double pp[7] );
00402         void get_w_tca_params( double& RperG, double& BperG );
00403 
00404 /* round-trip tests for verifying internal computations
00405    return inv(fwd(arg)) - arg or fwd(inv(arg)) - arg.
00406 */
00407   // currently valid max arguments
00408         double test_Tr_max();
00409         double test_Ta_max();
00410         double test_Sr_max();
00411         double test_Sa_max();
00412   // trig fns
00413         double test_Trar( double v );
00414         double test_Tara( double v );
00415   // spline tables
00416         double test_Srar( double v );
00417         double test_Sara( double v );
00418 
00419  private:
00420   /* primary design parameters */
00421         lens_kind lens; // selects a model
00422         camera_kind camera;     // selects a model
00423         double FL_mm;           // nominal focal length
00424         double FOV_deg;         // nominal view cone angle
00425         int hPix_num, vPix_num;  // nominal image dimensions
00426         double hPix_mm, vPix_mm; // horizontal & vertical resolutions
00427     /* derived parameters */
00428         double hFLpix, vFLpix;  // focal lengths in pixels
00429     double Rmaxrad;    // largest valid R = radius/flpix
00430         /* derived limits */
00431         double min_sina;
00432         double min_tana;
00433         double hc_s_lim, vc_s_lim;
00434 
00435 /* calibrated parameters   */
00436         double f_scale[2];      // h, v multipliers for FLmm
00437         double c_shift[2];      // optical center offsets
00438         double rad_params[7];   // spline params
00439         double tca_params[2];   // red/grn, blu/grn radius
00440         double proj_sina;       // a in A = a * sin( R / a )
00441         double proj_tana;       // a in A = a * tan( R / a )
00442         double proj_tanf;       // fraction tan in model
00443   /* volatile working parameters (optimizable) 
00444      note real working rad params are inside radfn
00445          note sina, tana, tanf are optimized indirectly
00446                 via w_proj_u and w_proj_v.
00447   */
00448         double w_hFLpix, w_vFLpix;      // scaled focal lengths
00449         double w_maxrad;
00450         double w_f_scale[2];
00451         double w_c_shift[2];
00452         double w_rad_params[7]; 
00453         double w_tca_params[2]; 
00454         double w_proj_sina;     
00455         double w_proj_tana;     
00456         double w_proj_tanf;     
00457   // actual projection params 
00458         double w_proj_u, w_proj_v;
00459         void setup_proj();  // applies w_proj_u, v
00460 
00461   /* state */
00462         const char * errmsg;
00463         bool model_valid;
00464         bool spline_valid;
00465         int color;      // 0: red, 1: green, 2: blue
00466   // these deternine which params get optimized...
00467         int opt_rad_k;  // number of radius fn levels, 0:4
00468         bool opt_center;        
00469         bool opt_scale;
00470         bool opt_tca;
00471   /* cubic spline tables */
00472 #define NSpts 110
00473         double r_val[ NSpts ], a_val[ NSpts ];
00474         double Ca2r[ 4 * NSpts ], 
00475                    Cr2a[ 4 * NSpts ];
00476   // validate model and fill the tables
00477         bool setSpline();
00478   // partial radial functions
00479   // ideal fns computed direct from model parameters
00480         double AofR( double R );
00481         double RofA( double A );        
00482   /* internal radius correction function
00483     Note yofx() is real(ideal), xofy() is ideal(real); 
00484         they cover 0:1 in both x and y
00485   */
00486         ssifn radfn;
00487 };
00488 
00489 #endif  //ndef _LENSFUNC_H
00490 

Generated on Mon Sep 20 01:01:27 2010 for Hugintrunk by doxygen 1.3.9.1