Homography.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.7.0
// ----------------------------------------------------------------------------
// pcl/Homography.h - Released 2024-06-18T15:48:54Z
// ----------------------------------------------------------------------------
// This file is part of the PixInsight Class Library (PCL).
// PCL is a multiplatform C++ framework for development of PixInsight modules.
//
// Copyright (c) 2003-2024 Pleiades Astrophoto S.L. All Rights Reserved.
//
// Redistribution and use in both source and binary forms, with or without
// modification, is permitted provided that the following conditions are met:
//
// 1. All redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//
// 2. All redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//
// 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
//    of their contributors, may be used to endorse or promote products derived
//    from this software without specific prior written permission. For written
//    permission, please contact info@pixinsight.com.
//
// 4. All products derived from this software, in any form whatsoever, must
//    reproduce the following acknowledgment in the end-user documentation
//    and/or other materials provided with the product:
//
//    "This product is based on software from the PixInsight project, developed
//    by Pleiades Astrophoto and its contributors (https://pixinsight.com/)."
//
//    Alternatively, if that is where third-party acknowledgments normally
//    appear, this acknowledgment must be reproduced in the product itself.
//
// THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
// INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
// DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
// ----------------------------------------------------------------------------
 
#ifndef __PCL_Homography_h
#define __PCL_Homography_h
 
 
#include <pcl/Defs.h>
#include <pcl/Diagnostics.h>
 
#include <pcl/Algebra.h>
#include <pcl/Array.h>
#include <pcl/Matrix.h>
#include <pcl/Point.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
class PCL_CLASS Homography
{
public:
 
   Homography()
      : m_H( Matrix::UnitMatrix( 3 ) )
   {
   }
 
   Homography( const Matrix& H )
      : m_H( H )
   {
   }
 
   template <class point_list1, class point_list2>
   Homography( const point_list1& P1, const point_list2& P2 )
      : m_H( DLT( P1, P2 ) )
   {
   }
 
   Homography( const Homography& ) = default;
 
   Homography& operator =( const Homography& ) = default;
 
   Homography( Homography&& ) = default;
 
   Homography& operator =( Homography&& ) = default;
 
   template <typename T1, typename T2>
   void Apply( T1& x, T2& y ) const
   {
      double fx = double( x );
      double fy = double( y );
      double w = m_H[2][0]*fx + m_H[2][1]*fy + m_H[2][2];
      PCL_CHECK( 1 + w != 1 )
      x = T1( (m_H[0][0]*fx + m_H[0][1]*fy + m_H[0][2])/w );
      y = T2( (m_H[1][0]*fx + m_H[1][1]*fy + m_H[1][2])/w );
   }
 
   template <typename T>
   GenericPoint<T>& Apply( GenericPoint<T>& p ) const
   {
      Apply( p.x, p.y );
      return p;
   }
 
   DPoint operator ()( double x, double y ) const
   {
      DPoint p( x, y );
      Apply( p.x, p.y );
      return p;
   }
 
   template <typename T>
   DPoint operator ()( const GenericPoint<T>& p ) const
   {
      return operator ()( double( p.x ), double( p.y ) );
   }
 
   Homography Inverse() const
   {
      return Homography( m_H.Inverse() );
   }
 
   const Matrix& TransformationMatrix() const
   {
      return m_H;
   }
 
   operator const Matrix&() const
   {
      return TransformationMatrix();
   }
 
   bool IsValid() const
   {
      return !m_H.IsEmpty();
   }
 
   bool IsAffine() const
   {
      return     IsValid()
          &&     Abs( m_H[2][0] ) <= std::numeric_limits<double>::epsilon()
          &&     Abs( m_H[2][1] ) <= std::numeric_limits<double>::epsilon()
          && 1 - Abs( m_H[2][2] ) <= std::numeric_limits<double>::epsilon();
   }
 
   void EnsureUnique()
   {
      m_H.EnsureUnique();
   }
 
private:
 
   Matrix m_H;
 
   /*
    * Normalization of reference points.
    */
   struct NormalizedPoints
   {
      Array<DPoint> N; // the normalized points
      Matrix        T; // 3x3 normalization matrix
 
      template <class point_list>
      NormalizedPoints( const point_list& points )
      {
         /*
          * Calculate the centroid of the input set of points.
          */
         DPoint centroid( 0 );
         for ( const auto& p : points )
            centroid += DPoint( p.x, p.y );
         centroid /= points.Length();
 
         /*
          * Calculate mean centroid distance and move origin to the centroid.
          */
         double d0 = 0;
         for ( const auto& p : points )
         {
            double dx = double( p.x ) - centroid.x;
            double dy = double( p.y ) - centroid.y;
            N << DPoint( dx, dy );
            d0 += Sqrt( dx*dx + dy*dy );
         }
         d0 /= points.Length();
 
         /*
          * Scale point coordinates.
          */
         double scale = Const<double>::sqrt2()/d0;
         for ( auto& p : N )
            p *= scale;
 
         /*
          * Form the normalization matrix.
          */
         T = Matrix( scale, 0.0,   -scale*centroid.x,
                     0.0,   scale, -scale*centroid.y,
                     0.0,   0.0,   1.0 );
      }
   };
 
   /*
    * Implementation of the Direct Linear Transformation (DLT) method to
    * compute a normalized homography matrix.
    */
   template <class point_list1, class point_list2>
   static Matrix DLT( const point_list1& P1, const point_list2& P2 )
   {
      int n = int( Min( P1.Length(), P2.Length() ) );
      if ( n < 4 )
         throw Error( "Homography::DLT(): Less than four points specified." );
 
      /*
       * Normalize all points.
       */
      NormalizedPoints p1( P1 );
      NormalizedPoints p2( P2 );
 
      /*
       * Setup cross product matrix A.
       */
      Matrix A( 2*n, 9 );
      for ( int i = 0; i < n; ++i )
      {
         double  x1 = p1.N[i].x;
         double  y1 = p1.N[i].y;
         double  x2 = p2.N[i].x;
         double  y2 = p2.N[i].y;
         double* A0 = A[2*i];
         double* A1 = A[2*i + 1];
         //double* A2 = A[3*i + 2]; <-- linearly dependent eqn.
 
         A0[0] = A0[1] = A0[2] = 0;
         A0[3] = -x1;
         A0[4] = -y1;
         A0[5] = -1;
         A0[6] =  y2*x1;
         A0[7] =  y2*y1;
         A0[8] =  y2;
 
         A1[0] =  x1;
         A1[1] =  y1;
         A1[2] =  1;
         A1[3] = A1[4] = A1[5] = 0;
         A1[6] = -x2*x1;
         A1[7] = -x2*y1;
         A1[8] = -x2;
 
         /*                         <-- linearly dependent eqn.
         A2[0] = -y2*x1;
         A2[1] = -y2*y1;
         A2[2] = -y2;
         A2[3] =  x2*x1;
         A2[4] =  x2*y1;
         A2[5] =  x2;
         A2[6] = A2[7] = A2[8] = 0;
         */
      }
 
      /*
       * SVD of cross product matrix.
       */
      InPlaceSVD svd( A );
 
      /*
       * For sanity, set to zero all insignificant singular values.
       */
      for ( int i = 0; i < svd.W.Length(); ++i )
         if ( Abs( svd.W[i] ) <= std::numeric_limits<double>::epsilon() )
            svd.W[i] = 0;
 
      /*
       * Locate the smallest nonzero singular value.
       */
      int i = svd.IndexOfSmallestSingularValue();
 
      /*
       * The components of the homography matrix are those of the smallest
       * eigenvector, i.e. the column vector of V corresponding to the smallest
       * singular value.
       */
      Matrix H( svd.V[0][i], svd.V[1][i], svd.V[2][i],
                svd.V[3][i], svd.V[4][i], svd.V[5][i],
                svd.V[6][i], svd.V[7][i], svd.V[8][i] );
 
      if ( Abs( H[2][2] ) <= std::numeric_limits<double>::epsilon() )
         throw Error( "Homography::DLT(): Singular matrix." );
 
      /*
       * Denormalize matrix components.
       */
      H = p2.T.Inverse() * H * p1.T;
 
      /*
       * Normalize matrix so that H[2][2] = 1.
       */
      return H *= 1/H[2][2];
   }
};
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif   // __PCL_Homography_h
 
// ----------------------------------------------------------------------------
// EOF pcl/Homography.h - Released 2024-06-18T15:48:54Z