USGS

Isis 3.0 Object Programmers' Reference

Home

ApolloMetricDistortionMap.cpp

Go to the documentation of this file.
00001 
00020 #include "ApolloMetricDistortionMap.h"
00021 #include "CameraFocalPlaneMap.h"
00022 
00023 using namespace std;
00024 namespace Isis {
00042   ApolloMetricDistortionMap::ApolloMetricDistortionMap(Camera *parent,
00043                                                        double xp, double yp,
00044                                                        double k1, double k2,
00045                                                        double k3, double j1,
00046                                                        double j2, double t0) :
00047                                                        CameraDistortionMap(parent) {
00048     p_xp = xp;
00049     p_yp = yp;
00050     p_k1 = k1;
00051     p_k2 = k2;
00052     p_k3 = k3;
00053     p_j1 = j1;
00054     p_j2 = j2;
00055     p_t0 = t0;
00056   }
00057 
00070   bool ApolloMetricDistortionMap::SetFocalPlane(const double dx, const double dy) {
00071     p_focalPlaneX = dx;
00072     p_focalPlaneY = dy;
00073 
00074     // reducing to principal point offset (xp,yp)
00075     double x = dx - p_xp;
00076     double y = dy - p_yp;
00077 
00078     // r is the distance between the principal point and the measured point on the image
00079     double rr = x * x + y * y;
00080     double rrrr = rr * rr;
00081 
00082     //  dr is the radial distortion contribution
00083     // -dt*sin(p_t0) is the decentering distortion contribution in the x-direction
00084     //  dt*cos(p_t0) is the decentering distortion contribution in the y-direction
00085     double dr = 1 + p_k1 * rr + p_k2 * rrrr + p_k3 * rr * rrrr;
00086     double dt = p_j1 * rr + p_j2 * rrrr;
00087 
00088     // image coordinates corrected for principal point, radial and decentering distortion
00089     p_undistortedFocalPlaneX = dr * x - dt * sin(p_t0);
00090     p_undistortedFocalPlaneY = dr * y + dt * cos(p_t0);
00091 
00092     return true;
00093   }
00094 
00095 
00108   bool ApolloMetricDistortionMap::SetUndistortedFocalPlane(const double ux, const double uy) {
00109     // image coordinates prior to introducing distortion
00110     p_undistortedFocalPlaneX = ux;
00111     p_undistortedFocalPlaneY = uy;
00112 
00113     double xt = ux;
00114     double yt = uy;
00115 
00116     double xx, yy, rr, rrrr, dr;
00117     double xdistortion, ydistortion;
00118     double xdistorted, ydistorted;
00119     double xprevious, yprevious;
00120     //  dr is the radial distortion contribution
00121     // -dt*sin(p_t0) is the decentering distortion contribution in the x-direction
00122     //  dt*cos(p_t0) is the decentering distortion contribution in the y-direction
00123     double dt;
00124 
00125     xprevious = 1000000.0;
00126     yprevious = 1000000.0;
00127 
00128     double tolerance = 0.000001;
00129 
00130     bool bConverged = false;
00131 
00132     // iterating to introduce distortion...
00133     // we stop when the difference between distorted coordinates
00134     // in successive iterations is at or below the given tolerance
00135     for(int i = 0; i < 50; i++) {
00136       xx = xt * xt;
00137       yy = yt * yt;
00138       rr = xx + yy;
00139       rrrr = rr * rr;
00140 
00141       // radial distortion
00142       //  dr is the radial distortion contribution
00143       // -dt*sin(p_t0) is the decentering distortion contribution in the x-direction
00144       //  dt*cos(p_t0) is the decentering distortion contribution in the y-direction
00145       dr = p_k1 * rr + p_k2 * rrrr + p_k3 * rr * rrrr;
00146 
00147       dt = p_j1 * rr + p_j2 * rrrr;
00148 
00149       // distortion at the current point location
00150       xdistortion = xt * dr - dt * sin(p_t0);
00151       ydistortion = yt * dr + dt * cos(p_t0);
00152 
00153       // updated image coordinates
00154       xt = ux - xdistortion;
00155       yt = uy - ydistortion;
00156 
00157       // distorted point corrected for principal point
00158       xdistorted = xt + p_xp;
00159       ydistorted = yt + p_yp;
00160 
00161       // check for convergence
00162       if((fabs(xt - xprevious) <= tolerance) && (fabs(yt - yprevious) <= tolerance)) {
00163         bConverged = true;
00164         break;
00165       }
00166 
00167       xprevious = xt;
00168       yprevious = yt;
00169     }
00170 
00171     if(bConverged) {
00172       p_focalPlaneX = xdistorted;
00173       p_focalPlaneY = ydistorted;
00174     }
00175 
00176     return bConverged;
00177   }
00178 }