StraightenPanorama.cpp

Go to the documentation of this file.
00001 // -*- c-basic-offset: 4 -*-
00026 #include "StraightenPanorama.h"
00027 
00028 #include <hugin_math/eig_jacobi.h>
00029 
00030 namespace HuginBase {
00031 
00032 Matrix3 StraightenPanorama::calcStraighteningRotation(const PanoramaData& panorama)
00033 {
00034     // landscape/non rotated portrait detection is not working correctly
00035     // should use the exif rotation tag but thats not stored anywhere currently...
00036     // 1: use y axis (image x axis), for normal image
00037     // 0: use z axis (image y axis), for non rotated portrait images
00038     //    (usually rotation is just stored in EXIF tag)
00039     std::vector<int> coord_idx;
00040 
00041     for (unsigned int i = 0; i < panorama.getNrOfImages(); i++) {
00042         SrcPanoImage img = panorama.getSrcImage(i);
00043         // BUG: need to read exif data here, since exif orientation is not
00044         // stored in Panorama data model
00045         double fl=0;
00046         double crop=0;
00047         img.readEXIF(fl, crop, false, false);
00048         double roll = img.getExifOrientation();
00049         if (roll == 90 || roll == 270 ) {
00050             coord_idx.push_back(2);
00051         } else {
00052             coord_idx.push_back(1);
00053         }
00054     }
00055 
00056     // build covariance matrix of X
00057     Matrix3 cov;
00058     unsigned int nrOfVariableImages=0;
00059 
00060     for (unsigned int i = 0; i < panorama.getNrOfImages(); i++) 
00061     {
00062         const SrcPanoImage & img=panorama.getImage(i);
00063         if(img.YawisLinked())
00064         {
00065             //only consider images which are not linked with the previous ones
00066             bool consider=true;
00067             for(unsigned int j=0; j<i; j++)
00068             {
00069                 if(img.YawisLinkedWith(panorama.getImage(j)))
00070                 {
00071                     consider=false;
00072                     break;
00073                 };
00074             };
00075             if(!consider)
00076                 continue;
00077         };
00078         double y = const_map_get(panorama.getImageVariables(i), "y").getValue();
00079         double p = const_map_get(panorama.getImageVariables(i), "p").getValue();
00080         double r = const_map_get(panorama.getImageVariables(i), "r").getValue();
00081         Matrix3 mat;
00082         mat.SetRotationPT(DEG_TO_RAD(y), DEG_TO_RAD(p), DEG_TO_RAD(r));
00083         nrOfVariableImages++;
00084         DEBUG_DEBUG("mat = " << mat);
00085         for (int j=0; j<3; j++) {
00086             for (int k=0; k<3; k++) {
00087                 cov.m[j][k] += mat.m[j][coord_idx[i]] * mat.m[k][coord_idx[i]];
00088             }
00089         }
00090     }
00091     cov /= nrOfVariableImages;
00092     DEBUG_DEBUG("cov = " << cov);
00093 
00094     // calculate eigenvalues and vectors
00095     Matrix3 eigvectors;
00096     double eigval[3];
00097     int eigvalIdx[3];
00098     int maxsweep = 100;
00099     int maxannil = 0;
00100     double eps = 1e-16;
00101     
00102     hugin_utils::eig_jacobi(3, cov.m, eigvectors.m, eigval, eigvalIdx, &maxsweep, &maxannil, &eps);
00103     
00104     DEBUG_DEBUG("Eigenvectors & eigenvalues:" << std::endl
00105                 << "V = " << eigvectors << std::endl
00106                 << "D = [" << eigval[0] << ", " << eigval[1] << ", " << eigval[2] << " ]"
00107                 << "idx = [" << eigvalIdx[0] << ", " << eigvalIdx[1] << ", " << eigvalIdx[2] << " ]");
00108     
00109     // get up vector, eigenvector with smallest eigenvalue
00110     Vector3 up;
00111     up.x = eigvectors.m[eigvalIdx[2]][0];
00112     up.y = eigvectors.m[eigvalIdx[2]][1];
00113     up.z = eigvectors.m[eigvalIdx[2]][2];
00114     
00115     // normalize vector
00116     up.Normalize();
00117     DEBUG_DEBUG("Up vector: up = " << up );
00118     
00119     double rotAngle = acos(up.Dot(Vector3(0,0,1)));
00120     if (rotAngle > M_PI/2) {
00121         // turn in shorter direction
00122         up *= -1;
00123         rotAngle = acos(up.Dot(Vector3(0,0,1)));
00124     }
00125     DEBUG_DEBUG("rotation Angle: " << rotAngle);
00126     
00127     // get rotation axis
00128     Vector3 rotAxis = up.Cross(Vector3(0,0,1));
00129     DEBUG_DEBUG("rotAxis = " << rotAngle);
00130     
00131     // calculate rotation matrix
00132     Matrix3 rotMat = GetRotationAroundU(rotAxis, -rotAngle);
00133     DEBUG_DEBUG("rotMat = " << rotMat);
00134 
00135     return rotMat;
00136 }
00137 
00138 } // namespace

Generated on Tue Sep 30 01:25:37 2014 for Hugintrunk by  doxygen 1.3.9.1