core/vgl/algo/vgl_norm_trans_3d.txx
Go to the documentation of this file.
00001 // This is core/vgl/algo/vgl_norm_trans_3d.txx
00002 #ifndef vgl_norm_trans_3d_txx_
00003 #define vgl_norm_trans_3d_txx_
00004 //:
00005 // \file
00006 
00007 #include "vgl_norm_trans_3d.h"
00008 #include <vgl/vgl_point_3d.h>
00009 #include <vnl/vnl_vector_fixed.h>
00010 #include <vcl_iostream.h>
00011 
00012 //--------------------------------------------------------------
00013 //
00014 
00015 //: Default constructor
00016 template <class T>
00017 vgl_norm_trans_3d<T>::vgl_norm_trans_3d() : vgl_h_matrix_3d<T>()
00018 {
00019 }
00020 
00021 //: Copy constructor
00022 template <class T>
00023 vgl_norm_trans_3d<T>::vgl_norm_trans_3d(const vgl_norm_trans_3d<T>& M)
00024 : vgl_h_matrix_3d<T>(M)
00025 {
00026 }
00027 
00028 
00029 //: Constructor from vcl_istream
00030 template <class T>
00031 vgl_norm_trans_3d<T>::vgl_norm_trans_3d(vcl_istream& s)
00032 : vgl_h_matrix_3d<T>(s)
00033 {
00034 }
00035 
00036 //: Constructor from file
00037 template <class T>
00038 vgl_norm_trans_3d<T>::vgl_norm_trans_3d(char const* filename)
00039 : vgl_h_matrix_3d<T>(filename)
00040 {
00041 }
00042 
00043 //--------------------------------------------------------------
00044 //: Constructor
00045 template <class T>
00046 vgl_norm_trans_3d<T>::vgl_norm_trans_3d(vnl_matrix_fixed<T,4,4> const& M)
00047 : vgl_h_matrix_3d<T>(M)
00048 {
00049 }
00050 
00051 //--------------------------------------------------------------
00052 //: Constructor
00053 template <class T>
00054 vgl_norm_trans_3d<T>::vgl_norm_trans_3d(const T* H)
00055 : vgl_h_matrix_3d<T>(H)
00056 {
00057 }
00058 
00059 //: Destructor
00060 template <class T>
00061 vgl_norm_trans_3d<T>::~vgl_norm_trans_3d()
00062 {
00063 }
00064 
00065 // == OPERATIONS ==
00066 //----------------------------------------------------------------
00067 //  Get the normalizing transform for a set of points
00068 // 1) Compute the center of gravity and form the normalizing
00069 //    transformation matrix
00070 // 2) Transform the point set to a temporary collection
00071 // 3) Compute the average point radius
00072 // 4) Complete the normalizing transform
00073 template <class T>
00074 bool vgl_norm_trans_3d<T>::
00075 compute_from_points(vcl_vector<vgl_homg_point_3d<T> > const& points)
00076 {
00077   T cx, cy, cz, radius;
00078   this->center_of_mass(points, cx, cy, cz);
00079   this->t12_matrix_.set_identity();
00080   this->t12_matrix_.put(0,3, -cx);
00081   this->t12_matrix_.put(1,3, -cy);
00082   this->t12_matrix_.put(2,3, -cz);
00083   vcl_vector<vgl_homg_point_3d<T> > temp;
00084   for (typename vcl_vector<vgl_homg_point_3d<T> >::const_iterator
00085        pit = points.begin(); pit != points.end(); pit++)
00086   {
00087     vgl_homg_point_3d<T> p((*this)(*pit));
00088     temp.push_back(p);
00089   }
00090   //Points might be coincident
00091   if (!this->scale_xyzroot2(temp, radius))
00092     return false;
00093   T scale = 1/radius;
00094   this->t12_matrix_.put(0,0, scale);
00095   this->t12_matrix_.put(1,1, scale);
00096   this->t12_matrix_.put(2,2, scale);
00097   this->t12_matrix_.put(0,3, -cx*scale);
00098   this->t12_matrix_.put(1,3, -cy*scale);
00099   this->t12_matrix_.put(2,3, -cz*scale);
00100   return true;
00101 }
00102 
00103 //-------------------------------------------------------------------
00104 // Find the center of a point cloud
00105 //
00106 template <class T>
00107 void vgl_norm_trans_3d<T>::
00108 center_of_mass(vcl_vector<vgl_homg_point_3d<T> > const& in,
00109                T& cx, T& cy, T& cz)
00110 {
00111   T cog_x = 0;
00112   T cog_y = 0;
00113   T cog_z = 0;
00114   T cog_count = 0.0;
00115   T tol = 1e-06;
00116   unsigned n = in.size();
00117   for (unsigned i = 0; i < n; ++i)
00118   {
00119     if (in[i].ideal(tol))
00120       continue;
00121     vgl_point_3d<T> p(in[i]);
00122     T x = p.x();
00123     T y = p.y();
00124     T z = p.z();
00125     cog_x += x;
00126     cog_y += y;
00127     cog_z += z;
00128     ++cog_count;
00129   }
00130   if (cog_count > 0)
00131   {
00132     cog_x /= cog_count;
00133     cog_y /= cog_count;
00134     cog_z /= cog_count;
00135   }
00136   //output result
00137   cx = cog_x;
00138   cy = cog_y;
00139   cz = cog_z;
00140 }
00141 
00142 //-------------------------------------------------------------------
00143 // Find the mean distance of the input pointset. Assumed to have zero mean
00144 //
00145 template <class T>
00146 bool vgl_norm_trans_3d<T>::
00147 scale_xyzroot2(vcl_vector<vgl_homg_point_3d<T> > const& in, T& radius)
00148 {
00149   T magnitude = T(0);
00150   int numfinite = 0;
00151   T tol = T(1e-06);
00152   radius = T(0);
00153   for (unsigned i = 0; i < in.size(); ++i)
00154   {
00155     if (in[i].ideal(tol))
00156       continue;
00157     vgl_point_3d<T> p(in[i]);
00158     vnl_vector_fixed<T, 3> v(p.x(), p.y(), p.z());
00159     magnitude += v.magnitude();
00160     ++numfinite;
00161   }
00162 
00163   if (numfinite > 0)
00164   {
00165     radius = magnitude / numfinite;
00166     return radius>=tol;
00167   }
00168   return false;
00169 }
00170 
00171 //----------------------------------------------------------------------------
00172 #undef VGL_NORM_TRANS_3D_INSTANTIATE
00173 #define VGL_NORM_TRANS_3D_INSTANTIATE(T) \
00174 template class vgl_norm_trans_3d<T >
00175 
00176 #endif // vgl_norm_trans_3d_txx_