contrib/oxl/mvl/TriTensor.h Source File

Go to the documentation of this file.
00001 // This is oxl/mvl/TriTensor.h
00002 #ifndef TriTensor_h_
00003 #define TriTensor_h_
00004 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00005 #pragma interface
00006 #endif
00007 //:
00008 // \file
00009 // \brief The trifocal tensor
00010 //
00011 // A class to hold a Trifocal Tensor and perform common operations, such as
00012 // point and line transfer, coordinate-frame transformation and I/O.
00013 //
00014 // \author
00015 //             Paul Beardsley, 29.03.96
00016 //             Oxford University, UK
00017 //
00018 // \verbatim
00019 //  Modifications:
00020 //   AWF - Added composition, transformation, homography generation.
00021 //   Peter Vanroose - 11 Mar 97 - added operator==
00022 //   Peter Vanroose - 22 Jun 03 - added vgl interface
00023 // \endverbatim
00024 //
00025 //------------------------------------------------------------------------------
00026 
00027 #include <vcl_vector.h>
00028 #include <vcl_iosfwd.h>
00029 
00030 #include <vbl/vbl_array_3d.h>
00031 
00032 #include <vnl/vnl_matrix.h>
00033 #include <vnl/vnl_vector.h>
00034 #include <vnl/vnl_double_3.h>
00035 #include <vnl/vnl_double_3x3.h>
00036 
00037 #include <vgl/vgl_fwd.h>
00038 #include <vgl/algo/vgl_algo_fwd.h>
00039 #include <vgl/algo/vgl_p_matrix.h>
00040 
00041 #include <mvl/HomgLine2D.h>
00042 #include <mvl/HomgLineSeg2D.h>
00043 #include <mvl/HomgPoint2D.h>
00044 
00045 #include <mvl/PMatrix.h>
00046 class HMatrix2D;
00047 class FMatrix;
00048 class PMatrix;
00049 class FManifoldProject;
00050 
00051 class TriTensor
00052 {
00053   // Data Members------------------------------------------------------------
00054 
00055   vbl_array_3d<double> T;
00056 
00057   // Caches for various computed quantities
00058   mutable const HomgPoint2D* e12_;
00059   mutable const HomgPoint2D* e13_;
00060 
00061   mutable const FManifoldProject* fmp12_;
00062   mutable const FManifoldProject* fmp13_;
00063   mutable const FManifoldProject* fmp23_;
00064 
00065  public:
00066 
00067   // Constructors/Initializers/Destructors-----------------------------------
00068 
00069   TriTensor();
00070   TriTensor(const TriTensor&);
00071   // Construct from 27-element vector
00072   TriTensor(const double *tritensor_array);
00073   TriTensor(const PMatrix& P1, const PMatrix& P2, const PMatrix& P3);
00074   TriTensor(const PMatrix& P2, const PMatrix& P3);
00075   TriTensor(const vnl_matrix<double>& T1, const vnl_matrix<double>& P2, const vnl_matrix<double>& P3);
00076  ~TriTensor();
00077 
00078   // Data Access-------------------------------------------------------------
00079 
00080   TriTensor& operator=(const TriTensor&);
00081   bool operator==(TriTensor const& p) const {
00082     for (int i=0;i<3;++i) for (int j=0;j<3;++j) for (int k=0;k<3;++k) if (p(i,j,k)!=T(i,j,k)) return false;
00083     return true; }
00084   double& operator() (unsigned int i1, unsigned int i2, unsigned int i3) { return T(i1,i2,i3); }
00085   double operator() (unsigned int i1, unsigned int i2, unsigned int i3) const { return T(i1,i2,i3); }
00086 
00087   void set(unsigned int i1, unsigned int i2, unsigned int i3, double value);
00088 
00089   void set(const double* vec);
00090   void set(const vnl_matrix<double>& tvector); // 27x1 matrix
00091   void convert_to_vector(vnl_matrix<double>* tvector) const; // 27x1 matrix
00092 
00093   void set(const PMatrix& P1, const PMatrix& P2, const PMatrix& P3);
00094   void set(const PMatrix& P2, const PMatrix& P3);
00095   void set(const vnl_matrix<double>& T1, const vnl_matrix<double>& T2, const vnl_matrix<double>& T3);
00096 
00097   // Data Control------------------------------------------------------------
00098 
00099   vgl_homg_point_2d<double> image1_transfer(vgl_homg_point_2d<double> const& point2,
00100                                             vgl_homg_point_2d<double> const& point3,
00101                                             vgl_homg_point_2d<double> corrected[] = 0) const;
00102   vgl_homg_point_2d<double> image2_transfer(vgl_homg_point_2d<double> const& point1,
00103                                             vgl_homg_point_2d<double> const& point3,
00104                                             vgl_homg_point_2d<double> corrected[] = 0) const;
00105   vgl_homg_point_2d<double> image3_transfer(vgl_homg_point_2d<double> const& point1,
00106                                             vgl_homg_point_2d<double> const& point2,
00107                                             vgl_homg_point_2d<double> corrected[] = 0) const;
00108 
00109   HomgPoint2D image1_transfer(HomgPoint2D const& point2,
00110                               HomgPoint2D const& point3,
00111                               HomgPoint2D corrected[] = 0) const;
00112   HomgPoint2D image2_transfer(HomgPoint2D const& point1,
00113                               HomgPoint2D const& point3,
00114                               HomgPoint2D corrected[] = 0) const;
00115   HomgPoint2D image3_transfer(HomgPoint2D const& point1,
00116                               HomgPoint2D const& point2,
00117                               HomgPoint2D corrected[] = 0) const;
00118 
00119   vgl_homg_point_2d<double> image1_transfer_qd(vgl_homg_point_2d<double> const& point2,
00120                                                vgl_homg_point_2d<double> const& point3) const;
00121   vgl_homg_point_2d<double> image2_transfer_qd(vgl_homg_point_2d<double> const& point1,
00122                                                vgl_homg_point_2d<double> const& point3) const;
00123   vgl_homg_point_2d<double> image3_transfer_qd(vgl_homg_point_2d<double> const& point1,
00124                                                vgl_homg_point_2d<double> const& point2) const;
00125 
00126   HomgPoint2D image1_transfer_qd(const HomgPoint2D& point2, const HomgPoint2D& point3) const;
00127   HomgPoint2D image2_transfer_qd(const HomgPoint2D& point1, const HomgPoint2D& point3) const;
00128   HomgPoint2D image3_transfer_qd(const HomgPoint2D& point1, const HomgPoint2D& point2) const;
00129 
00130   vgl_homg_point_2d<double> image1_transfer(vgl_homg_point_2d<double> const& point1,
00131                                             vgl_homg_line_2d <double> const& line2) const;
00132   vgl_homg_point_2d<double> image2_transfer(vgl_homg_point_2d<double> const& point1,
00133                                             vgl_homg_line_2d <double> const& line3) const;
00134   vgl_homg_point_2d<double> image3_transfer(vgl_homg_point_2d<double> const& point2,
00135                                             vgl_homg_line_2d <double> const& line3) const;
00136 
00137   HomgPoint2D image1_transfer(HomgPoint2D const& point1,
00138                               HomgLine2D  const& line2) const;
00139   HomgPoint2D image2_transfer(HomgPoint2D const& point1,
00140                               HomgLine2D  const& line3) const;
00141   HomgPoint2D image3_transfer(HomgPoint2D const& point2,
00142                               HomgLine2D  const& line3) const;
00143 
00144   vgl_homg_line_2d<double> image1_transfer(vgl_homg_line_2d<double> const& line2,
00145                                            vgl_homg_line_2d<double> const& line3) const;
00146   vgl_homg_line_2d<double> image2_transfer(vgl_homg_line_2d<double> const& line2,
00147                                            vgl_homg_line_2d<double> const& line3) const;
00148   vgl_homg_line_2d<double> image3_transfer(vgl_homg_line_2d<double> const& line2,
00149                                            vgl_homg_line_2d<double> const& line3) const;
00150 
00151   HomgLine2D image1_transfer(const HomgLine2D& line2, const HomgLine2D& line3) const;
00152   HomgLine2D image2_transfer(const HomgLine2D& line2, const HomgLine2D& line3) const;
00153   HomgLine2D image3_transfer(const HomgLine2D& line2, const HomgLine2D& line3) const;
00154 
00155   vgl_h_matrix_2d<double> get_hmatrix_31(vgl_homg_line_2d<double> const& line2) const;
00156   vgl_h_matrix_2d<double> get_hmatrix_21(vgl_homg_line_2d<double> const& line3) const;
00157 
00158   HMatrix2D get_hmatrix_31(const HomgLine2D& line2) const;
00159   HMatrix2D get_hmatrix_21(const HomgLine2D& line3) const;
00160 
00161   bool get_epipoles(vgl_homg_point_2d<double>& e2, vgl_homg_point_2d<double>& e3) const;
00162   bool get_epipoles(HomgPoint2D* e2, HomgPoint2D* e3) const;
00163   bool compute_epipoles() const; // mutable const
00164 
00165   HomgPoint2D get_epipole_12() const;
00166   HomgPoint2D get_epipole_13() const;
00167 
00168   FMatrix get_fmatrix_13() const;
00169   FMatrix get_fmatrix_12() const;
00170 
00171   FMatrix compute_fmatrix_23() const;
00172 
00173   const FManifoldProject* get_fmp12() const;
00174   const FManifoldProject* get_fmp23() const;
00175   const FManifoldProject* get_fmp13() const;
00176 
00177   void compute_P_matrices(const vnl_vector<double>& x, double alpha, double beta, PMatrix* P2, PMatrix* P3) const;
00178   void compute_P_matrices(const vnl_vector<double>& x, double alpha, PMatrix* P2, PMatrix* P3) const {
00179     compute_P_matrices(x,alpha,alpha, P2, P3);
00180   }
00181   void compute_P_matrices(const vnl_vector<double>& x, PMatrix* P2, PMatrix* P3) const {
00182     compute_P_matrices(x, 1, 1, P2, P3);
00183   }
00184   void compute_P_matrices(PMatrix* P2, PMatrix* P3) const {
00185     compute_P_matrices(vnl_double_3(1,1,1), 1, 1, P2, P3);
00186   }
00187   void compute_P_matrices(vgl_p_matrix<double> &P2, vgl_p_matrix<double> &P3) const {
00188     PMatrix Ptemp2, Ptemp3;
00189     compute_P_matrices(vnl_double_3(1,1,1), 1, 1, &Ptemp2, &Ptemp3);
00190     P2.set(Ptemp2.get_matrix());
00191     P3.set(Ptemp3.get_matrix());
00192   }
00193   void compute_caches();
00194   void clear_caches();
00195 
00196   // Utility Methods---------------------------------------------------------
00197   void get_constraint_lines_image3(vgl_homg_point_2d<double> const& p1,
00198                                    vgl_homg_point_2d<double> const& p2,
00199                                    vcl_vector<vgl_homg_line_2d<double> >& lines) const;
00200   void get_constraint_lines_image2(vgl_homg_point_2d<double> const& p1,
00201                                    vgl_homg_point_2d<double> const& p3,
00202                                    vcl_vector<vgl_homg_line_2d<double> >& lines) const;
00203   void get_constraint_lines_image1(vgl_homg_point_2d<double> const& p2,
00204                                    vgl_homg_point_2d<double> const& p3,
00205                                    vcl_vector<vgl_homg_line_2d<double> >& lines) const;
00206 
00207   void get_constraint_lines_image3(HomgPoint2D const& p1,
00208                                    HomgPoint2D const& p2,
00209                                    vcl_vector<HomgLine2D>* lines) const;
00210   void get_constraint_lines_image2(HomgPoint2D const& p1,
00211                                    HomgPoint2D const& p3,
00212                                    vcl_vector<HomgLine2D>* lines) const;
00213   void get_constraint_lines_image1(HomgPoint2D const& p2,
00214                                    HomgPoint2D const& p3,
00215                                    vcl_vector<HomgLine2D>* lines) const;
00216 
00217   //: Contract Tensor axis tensor_axis with first component of Matrix M.
00218   // That is:
00219   // For tensor_axis = 1,  Compute T_ijk = T_pjk M_pi
00220   // For tensor_axis = 2,  Compute T_ijk = T_ipk M_pj
00221   // For tensor_axis = 3,  Compute T_ijk = T_ijp M_pk
00222   TriTensor postmultiply(unsigned tensor_axis, const vnl_matrix<double>& M) const;
00223 
00224   //: Contract Tensor axis tensor_axis with second component of Matrix M.
00225   // That is:
00226   // For tensor_axis = 1,  Compute T_ijk = M_ip T_pjk
00227   // For tensor_axis = 2,  Compute T_ijk = M_jp T_ipk
00228   // For tensor_axis = 3,  Compute T_ijk = M_kp T_ijp
00229   TriTensor premultiply(unsigned tensor_axis, const vnl_matrix<double>& M) const;
00230 
00231   TriTensor postmultiply1(const vnl_matrix<double>& M) const;
00232   TriTensor postmultiply2(const vnl_matrix<double>& M) const;
00233   TriTensor postmultiply3(const vnl_matrix<double>& M) const;
00234 
00235   TriTensor premultiply1(const vnl_matrix<double>& M) const;
00236   TriTensor premultiply2(const vnl_matrix<double>& M) const;
00237   TriTensor premultiply3(const vnl_matrix<double>& M) const;
00238 
00239   vnl_double_3x3 dot1(const vnl_double_3& v) const;
00240   vnl_double_3x3 dot2(const vnl_double_3& v) const;
00241   vnl_double_3x3 dot3(const vnl_double_3& v) const;
00242   vnl_double_3x3 dot1t(const vnl_double_3& v) const;
00243   vnl_double_3x3 dot2t(const vnl_double_3& v) const;
00244   vnl_double_3x3 dot3t(const vnl_double_3& v) const;
00245 
00246   bool check_equal_up_to_scale(const TriTensor& that) const;
00247 
00248   // INTERNALS---------------------------------------------------------------
00249 
00250   //: C123 are line conditioning matrices.
00251   // If C * l = lhat, and l1 = T l2 l3, then lhat1 = That lhat2 lhat3
00252   TriTensor condition(const vnl_matrix<double>& line_1_denorm,
00253                       const vnl_matrix<double>& line_2_norm,
00254                       const vnl_matrix<double>& line_3_norm) const;
00255 
00256   TriTensor decondition(const vnl_matrix<double>& line_1_norm,
00257                         const vnl_matrix<double>& line_2_denorm,
00258                         const vnl_matrix<double>& line_3_denorm) const;
00259  private:
00260   void delete_caches() const; // mutable const
00261 };
00262 
00263 vcl_ostream& operator << (vcl_ostream&, const TriTensor& T);
00264 vcl_istream& operator >> (vcl_istream&, TriTensor& T);
00265 
00266 #endif // TriTensor_h_