contrib/mul/mbl/mbl_thin_plate_spline_3d.h
Go to the documentation of this file.
00001 // This is mul/mbl/mbl_thin_plate_spline_3d.h
00002 #ifndef mbl_thin_plate_spline_3d_h_
00003 #define mbl_thin_plate_spline_3d_h_
00004 //:
00005 // \file
00006 // \brief Construct thin plate spline to map 3D to 3D
00007 // \author Tim Cootes
00008 
00009 #include <vgl/vgl_point_3d.h>
00010 #include <vnl/vnl_vector.h>
00011 #include <vnl/vnl_matrix.h>
00012 #include <vsl/vsl_binary_io.h>
00013 #include <vcl_vector.h>
00014 #include <vcl_iosfwd.h>
00015 
00016 //=======================================================================
00017 //: Construct thin plate spline to map 3D to 3D.
00018 // I.e. does some mapping (x',y',z') = f(x,y,z). (See Booksteins work, e.g. IPMI 1993)
00019 // The warp is `guided' by a set of
00020 // landmarks p(0) .. p(n-1) in the source plane which are to be
00021 // mapped to a (possibly deformed) set q(0)..q(n-1) in the destination.
00022 // Thus the mapping is constrained so that f(p(i)) = q(i) for i = 0..n-1.
00023 // The points are given to the build() function to set up the object.
00024 //
00025 // If one wishes to map a set of source points to multiple target points,
00026 // use set_source_pts(src_pts);  then build(target_pts); for each target set.
00027 //
00028 // \code
00029 // vcl_vector<vgl_point_3d<double> > src_pts(n_points),dest_pts(n_points);
00030 //
00031 // // Fill src_pts and dest_pts
00032 // .....
00033 //
00034 // // Construct spline object
00035 // mbl_thin_plate_spline_3d tps;
00036 // tps.build(src_pts,dest_pts);
00037 //
00038 // // Apply to point p:
00039 // vgl_point_3d<double> p(1,2,3);
00040 // vgl_point_3d<double> new_p = tps(p);
00041 // \endcode
00042 class mbl_thin_plate_spline_3d
00043 {
00044   vnl_vector<double> Wx_,Wy_,Wz_;
00045   double Ax0_, AxX_, AxY_, AxZ_;
00046   double Ay0_, AyX_, AyY_, AyZ_;
00047   double Az0_, AzX_, AzY_, AzZ_;
00048   double energy_x_,energy_y_,energy_z_;
00049 
00050   vcl_vector<vgl_point_3d<double> > src_pts_;
00051 
00052   //: Used to estimate weights in set_source_points()
00053   vnl_matrix<double> L_inv_;
00054 
00055   //: Build from small number of points
00056   void build_pure_affine(const vcl_vector<vgl_point_3d<double> >& source_pts,
00057                          const vcl_vector<vgl_point_3d<double> >& dest_pts);
00058 
00059   //: Set parameters from vectors
00060   void set_params(const vnl_vector<double>& W1,
00061                   const vnl_vector<double>& W2,
00062                   const vnl_vector<double>& W3);
00063 
00064   void set_up_rhs(vnl_vector<double>& Bx,
00065                   vnl_vector<double>& By,
00066                   vnl_vector<double>& Bz,
00067                   const vcl_vector<vgl_point_3d<double> >& dest_pts);
00068 
00069   //: Compute spline-bending energy
00070   void compute_energy(vnl_vector<double>& W1,
00071                       vnl_vector<double>& W2,
00072                       vnl_vector<double>& W3,
00073                       const vnl_matrix<double>& L);
00074 
00075  public:
00076 
00077   //: Dflt ctor
00078   mbl_thin_plate_spline_3d();
00079 
00080   //: Destructor
00081   virtual ~mbl_thin_plate_spline_3d();
00082 
00083   //: Sets up internal transformation to map source_pts onto dest_pts
00084   void build(const vcl_vector<vgl_point_3d<double> >& source_pts,
00085              const vcl_vector<vgl_point_3d<double> >& dest_pts,
00086              bool compute_the_energy=false);
00087 
00088   //: Define source point positions
00089   //  Performs pre-computations so that build(dest_points) can be
00090   //  called multiple times efficiently
00091   void set_source_pts(const vcl_vector<vgl_point_3d<double> >& source_pts);
00092 
00093   //: Sets up internal transformation to map source_pts onto dest_pts
00094   void build(const vcl_vector<vgl_point_3d<double> >& dest_pts);
00095 
00096   //: Return transformed version of (x,y,z)
00097   vgl_point_3d<double>  operator()(double x, double y, double z) const;
00098 
00099   //: Return transformed version of (x,y,z)
00100   vgl_point_3d<double>  operator()(const vgl_point_3d<double>&  p) const
00101   { return operator()(p.x(),p.y(),p.z()); }
00102 
00103   //: Bending energy of X component (zero for pure affine)
00104   //  A measure of total amount of non-linear deformation
00105   double bendingEnergyX() const { return energy_x_; }
00106 
00107   //: Bending energy of Y component (zero for pure affine)
00108   //  A measure of total amount of non-linear deformation
00109   double bendingEnergyY() const { return energy_y_; }
00110 
00111   //: Bending energy of Z component (zero for pure affine)
00112   //  A measure of total amount of non-linear deformation
00113   double bendingEnergyZ() const { return energy_z_; }
00114 
00115   //: Version number for I/O
00116   short version_no() const;
00117 
00118   //: Print class to os
00119   void print_summary(vcl_ostream& os) const;
00120 
00121   //: Save class to binary file stream
00122   void b_write(vsl_b_ostream& bfs) const;
00123 
00124   //: Load class from binary file stream
00125   void b_read(vsl_b_istream& bfs);
00126 
00127   //: Comparison operator
00128   bool operator==(const mbl_thin_plate_spline_3d& tps) const;
00129 };
00130 
00131 //: Binary file stream output operator for class reference
00132 void vsl_b_write(vsl_b_ostream& bfs, const mbl_thin_plate_spline_3d& b);
00133 
00134 //: Binary file stream input operator for class reference
00135 void vsl_b_read(vsl_b_istream& bfs, mbl_thin_plate_spline_3d& b);
00136 
00137 //: Stream output operator for class reference
00138 vcl_ostream& operator<<(vcl_ostream& os,const mbl_thin_plate_spline_3d& b);
00139 
00140 #endif
00141 
00142