ENH: Adding a helper class to hold duplicated code computing p/xs offset for...

ENH: Adding a helper class to hold duplicated code computing p/xs offset for Pleiades sensor bundles

Code/Projections/otbPleiadesPToXSAffineTransformCalculator.cxx

+/*=========================================================================
+
+  Program:   ORFEO Toolbox
+  Language:  C++
+  Date:      $Date$
+  Version:   $Revision$
+
+
+  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+  See OTBCopyright.txt for details.
+
+
+     This software is distributed WITHOUT ANY WARRANTY; without even
+     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+     PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbPleiadesPToXSAffineTransformCalculator__cxx
+#define __otbPleiadesPToXSAffineTransformCalculator__cxx
+
+#include "otbPleiadesPToXSAffineTransformCalculator.h"
+#include "otbPleiadesImageMetadataInterface.h"
+#include "otbDateTimeAdapter.h"
+#include "otbImageKeywordlist.h"
+#include "itkMetaDataObject.h"
+
+namespace otb {
+
+bool 
+PleiadesPToXSAffineTransformCalculator
+::CanCompute(const itk::ImageBase<2> * panchromaticImage, const itk::ImageBase<2> * xsImage)
+{
+  bool isPanPHR = false;
+  bool isXSPHR = false;
+      
+  otb::PleiadesImageMetadataInterface::Pointer phrIMI =
+    otb::PleiadesImageMetadataInterface::New();
+  phrIMI->SetMetaDataDictionary(panchromaticImage->GetMetaDataDictionary());
+  isPanPHR = phrIMI->CanRead();
+      
+  phrIMI->SetMetaDataDictionary(xsImage->GetMetaDataDictionary());
+  isXSPHR = phrIMI->CanRead();
+      
+  if (isPanPHR && isXSPHR)
+    {
+    ImageKeywordlist kwlPan;
+    ImageKeywordlist kwlXS;
+        
+    itk::ExposeMetaData<ImageKeywordlist>(
+      panchromaticImage->GetMetaDataDictionary(),
+      MetaDataKey::OSSIMKeywordlistKey,
+      kwlPan);
+        
+    itk::ExposeMetaData<ImageKeywordlist>(
+      xsImage->GetMetaDataDictionary(),
+      MetaDataKey::OSSIMKeywordlistKey,
+      kwlXS);
+    
+    // Get geometric processing
+    std::string panProcessing = kwlPan.GetMetadataByKey("support_data.processing_level");
+    std::string xsProcessing = kwlXS.GetMetadataByKey("support_data.processing_level");
+        
+    if (panProcessing.compare("SENSOR") == 0 &&
+        xsProcessing.compare("SENSOR") == 0)
+      {
+      
+      std::string pid = kwlPan.GetMetadataByKey("image_id");
+      std::string xsid = kwlXS.GetMetadataByKey("image_id");
+      
+      pid  = pid.substr(0,pid.size()-4);
+      xsid = xsid.substr(0,xsid.size()-4);
+
+      if(pid == xsid)
+        {
+        return true;
+        }
+      }
+    }
+  return false;
+}
+
+PleiadesPToXSAffineTransformCalculator
+::TransformType::Pointer
+PleiadesPToXSAffineTransformCalculator
+::Compute(const itk::ImageBase<2> * panchromaticImage, const itk::ImageBase<2> * xsImage)
+{
+  if(!CanCompute(panchromaticImage,xsImage))
+    {
+    itkGenericExceptionMacro("Can not compute affine transform between images, they do not correspond to Pleiades sensor bundle.");
+    }
+
+    ImageKeywordlist kwlPan;
+    ImageKeywordlist kwlXS;
+    
+    itk::ExposeMetaData<ImageKeywordlist>(
+      panchromaticImage->GetMetaDataDictionary(),
+      MetaDataKey::OSSIMKeywordlistKey,
+      kwlPan);
+        
+    itk::ExposeMetaData<ImageKeywordlist>(
+      xsImage->GetMetaDataDictionary(),
+      MetaDataKey::OSSIMKeywordlistKey,
+      kwlXS);
+
+    // Compute time delta
+    std::string strStartTimePan = kwlPan.GetMetadataByKey("support_data.time_range_start");
+    std::string strStartTimeXS = kwlXS.GetMetadataByKey("support_data.time_range_start");
+      
+    DateTimeAdapter::Pointer startTimePan = DateTimeAdapter::New();
+    DateTimeAdapter::Pointer startTimeXS = DateTimeAdapter::New();
+      
+    startTimePan->SetFromIso8601(strStartTimePan);
+    startTimeXS->SetFromIso8601(strStartTimeXS);
+      
+    double timeDelta = startTimeXS->GetDeltaInSeconds(startTimePan);
+      
+    // Retrieve line period in Pan
+    std::string tmpStr = kwlPan.GetMetadataByKey("support_data.line_period");
+    double linePeriodPan = atof(tmpStr.c_str());
+      
+    // Retrieve column start
+    tmpStr = kwlPan.GetMetadataByKey("support_data.swath_first_col");
+    int colStartPan = atoi(tmpStr.c_str());
+    tmpStr = kwlXS.GetMetadataByKey("support_data.swath_first_col");
+    int colStartXS = atoi(tmpStr.c_str());
+      
+    // Compute shift between MS and P (in Pan pixels)
+    // in order to keep the top left corners unchanged, apply an
+    // additional offset of (3/2) panchro pixels, or 0.375 xs pixels
+    int lineShift_MS_P =int(vcl_floor((timeDelta/(linePeriodPan/1000))  + 0.5)) + 0.375;
+    int colShift_MS_P = colStartXS*4 - colStartPan-4 + 0.375;
+     
+    // Apply the scaling
+    typedef itk::ScalableAffineTransform<double, 2>  TransformType;
+    TransformType::Pointer transform = TransformType::New();
+    transform->Scale(4.0);
+
+    // Apply the offset
+    TransformType::OutputVectorType offset;
+    offset[0] = static_cast<double>(colShift_MS_P);
+    offset[1] = static_cast<double>(lineShift_MS_P);
+    transform->Translate(offset);
+    
+    // Invert the transform to get the P to XS transform
+    TransformType::Pointer realTransform = TransformType::New();
+    transform->GetInverse(realTransform);
+    
+    return realTransform;
+    }
+
+} // End namespace otb
+
+#endif
+
+

Code/Projections/otbPleiadesPToXSAffineTransformCalculator.h

+/*=========================================================================
+
+  Program:   ORFEO Toolbox
+  Language:  C++
+  Date:      $Date$
+  Version:   $Revision$
+
+
+  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+  See OTBCopyright.txt for details.
+
+
+     This software is distributed WITHOUT ANY WARRANTY; without even
+     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+     PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbPleiadesPToXSAffineTransformCalculator__h
+#define __otbPleiadesPToXSAffineTransformCalculator__h
+
+#include "itkLightObject.h"
+#include "itkObjectFactory.h"
+#include "itkScalableAffineTransform.h"
+#include "itkImageBase.h"
+
+namespace otb {
+/**
+ * \class PleiadesPToXSAffineTransformCalculator
+ * \brief Compute the affine transform linking P and XS pixel position for Pleiades sensor bundles
+ * 
+ * Pleiades sensor bundle are exactly homotetic, it is therefore
+ * possible to corregister the pancrhomatic and multispectral images
+ * with a simple affine transform without using any sensor
+ * modelling. This yelds a very accurate corregistration and avoid the
+ * use of a DSM which may cause registration errors due to height errors.
+ * 
+ * This calculator is a helper class to build the affine transform. It
+ * consists in only two static methods: one to check if the transform
+ * calculation applies to the given pair of images, the other to
+ * actually estimate the transfrom.
+ * 
+ * The estimated transform returned by the latter transforms
+ * pancrhomatic image positions to multispectral image positions. If
+ * the inverse transform is needed, one can call the GetInverse()
+ * method of the transform to retrieve it.
+ * 
+ */
+class ITK_EXPORT PleiadesPToXSAffineTransformCalculator
+  : public itk::LightObject
+{
+public:
+  typedef PleiadesPToXSAffineTransformCalculator Self;
+  typedef itk::LightObject                       Superclass;
+  typedef itk::SmartPointer<Self>                Pointer;
+  typedef itk::SmartPointer<const Self>          ConstPointer;
+
+  itkTypeMacro(PleiadesPToXSAffineTransformCalculator,itk::LightObject);
+  
+  typedef itk::ScalableAffineTransform<double,2> TransformType;
+  
+  /** 
+   * This function checks if the transform calculation applies to the
+   * given pair of images. Checked items are:
+   * - Both images are sucessfully undertood by OTB as Pléiades images,
+   * - Both images processing level is SENSOR",
+   * - XS and Pan ids (except last 4 letters) are equal.
+   * \return True if the calculation applies
+   */
+  static bool CanCompute(const itk::ImageBase<2> * panchromaticImage, const itk::ImageBase<2> * xsImage);
+     
+
+  /**
+   * This function computes the transform for a pair of image. Note
+   * that the CanCompute() method is first called, and that an
+   * exception will be raised if computation can not be done.
+   * 
+   * This function reads both images support data and builds a
+   * transform that will exactly coregister the images.
+   * 
+   * \return The computed transform
+   */
+  static TransformType::Pointer Compute(const itk::ImageBase<2> * panchromaticImage, const itk::ImageBase<2> * xsImage);
+                  
+private:
+  PleiadesPToXSAffineTransformCalculator(); // purposely not implemented
+  PleiadesPToXSAffineTransformCalculator(const Self &); // purposely not implemented
+  void operator =(const Self&); // purposely not implemented
+};
+
+
+} // End namespace otb
+
+
+#endif