From bfb10a47adda5f37cd3040701de66e294dbf3446 Mon Sep 17 00:00:00 2001
From: Otmane Lahlou <otmane.lahlou@c-s.fr>
Date: Tue, 21 Sep 2010 17:59:29 +0200
Subject: [PATCH] ENH : use the new orhtorectification framework

---
 Projections/otbBundleToPerfectSensor.cxx | 105 ++++++-----------------
 1 file changed, 28 insertions(+), 77 deletions(-)

diff --git a/Projections/otbBundleToPerfectSensor.cxx b/Projections/otbBundleToPerfectSensor.cxx
index 31da097bba..38bfec79d7 100644
--- a/Projections/otbBundleToPerfectSensor.cxx
+++ b/Projections/otbBundleToPerfectSensor.cxx
@@ -26,29 +26,17 @@ See OTBCopyright.txt for details.
 
 #include <iostream>
 
-#include "otbMetaDataKey.h"
 #include "otbCommandLineArgumentParser.h"
 #include "otbImageFileReader.h"
 #include "otbStreamingImageFileWriter.h"
-#include "otbImageFileWriter.h"
-#include "otbGenericRSTransform.h"
 #include "otbImage.h"
-#include "otbDEMHandler.h"
 #include "otbVectorImage.h"
-#include "otbMacro.h"
+#include "otbGenericRSResampleImageFilter.h"
 
 #include "itkExceptionObject.h"
-#include "itkMacro.h"
-#include "itkTransform.h"
 
 #include "otbStandardFilterWatcher.h"
-#include "itkTransformToDeformationFieldSource.h"
-#include "otbStreamingWarpImageFilter.h"
-#include "itkLinearInterpolateImageFunction.h"
 #include "otbSimpleRcsPanSharpeningFusionImageFilter.h"
-
-#include "otbMultiChannelExtractROI.h"
-
 #include "itkPixelBuilder.h"
 
 #include "init/ossimInit.h"
@@ -67,7 +55,7 @@ int main(int argc, char* argv[])
 
     parser->SetProgramDescription("Using available image metadata to determine the sensor model, computes a cartographic projection of the image");
     parser->AddOutputImage();
-//    parser->AddOption("--DEMDirectory","Directory were to find the DEM tiles","-dem",1,false);
+    parser->AddOption("--DEMDirectory","Directory were to find the DEM tiles","-dem",1,false);
     parser->AddOption("--NumStreamDivisions","Number of streaming divisions (optional)","-stream",1,false);
     parser->AddOption("--LocMapSpacing","Generate a coarser deformation field with the given spacing.","-lmSpacing",1,false);
     parser->AddOption("--InputPanchro","The input panchromatic image","-inP", 1,true);
@@ -103,16 +91,8 @@ int main(int argc, char* argv[])
     typedef otb::ImageFileReader<XsImageType>          XsReaderType;
     typedef otb::ImageFileReader<PanImageType>         PanReaderType;
     typedef otb::StreamingImageFileWriter<XsImageType> WriterType;
-    typedef otb::GenericRSTransform<>                  RSTransformType;
-    typedef RSTransformType::InputPointType   PointType;
-
-    typedef itk::FixedArray<double,2> DeformationValueType;
-    typedef otb::Image<DeformationValueType> DeformationFieldType;
-    typedef otb::StreamingImageFileWriter<DeformationFieldType> DeformationFieldWriterType;
 
-    typedef itk::TransformToDeformationFieldSource<DeformationFieldType,double> DeformationFieldGeneratorType;
-    typedef otb::StreamingWarpImageFilter<XsImageType,XsImageType,DeformationFieldType> WarpFilterType;
-    typedef itk::LinearInterpolateImageFunction<XsImageType,double> InterpolatorType;
+    typedef otb::GenericRSResampleImageFilter<XsImageType,XsImageType>  ResamplerType;
 
     typedef otb::SimpleRcsPanSharpeningFusionImageFilter<PanImageType,XsImageType,XsImageType> FusionFilterType;
 
@@ -124,84 +104,55 @@ int main(int argc, char* argv[])
     XsReaderType::Pointer xsreader= XsReaderType::New();
     xsreader->SetFileName(parseResult->GetParameterString("--InputXS"));
     xsreader->GenerateOutputInformation();
-
-    // Generate RS Transform
-    RSTransformType::Pointer rsTransform = RSTransformType::New();
-    rsTransform->SetOutputKeywordList(xsreader->GetOutput()->GetImageKeywordlist());
-    rsTransform->SetInputKeywordList(preader->GetOutput()->GetImageKeywordlist());
-
-//    if(parseResult->IsOptionPresent("--DEMDirectory"))
-//      {
-//      rsTransform->SetDEMDirectory(parseResult->GetParameterString("--DEMDirectory",0));
-//      gcp2sensor->SetUseDEM(true);
-//      otb::DEMHandler::Pointer demHandler = otb::DEMHandler::New();
-//      demHandler->OpenDEMDirectory(parseResult->GetParameterString("--DEMDirectory",0).c_str());
-//      gcp2sensor->SetDEMHandler(demHandler);
-//      }
-
-    rsTransform->InstanciateTransform();
-
-
+    
+    // Resample filter 
+    ResamplerType::Pointer    resampler = ResamplerType::New();
+    
+    // Add DEM if any
+    if(parseResult->IsOptionPresent("--DEMDirectory"))
+      {
+      resampler->SetDEMDirectory(parseResult->GetParameterString("--DEMDirectory",0));
+      }
+    
     // Set up output image informations
     XsImageType::SpacingType spacing = preader->GetOutput()->GetSpacing();
     XsImageType::IndexType start = preader->GetOutput()->GetLargestPossibleRegion().GetIndex();
     XsImageType::SizeType size = preader->GetOutput()->GetLargestPossibleRegion().GetSize();
     XsImageType::PointType origin = preader->GetOutput()->GetOrigin();
 
-    DeformationFieldGeneratorType::Pointer dfGenerator = DeformationFieldGeneratorType::New();
-    dfGenerator->SetOutputOrigin(origin);
-    dfGenerator->SetOutputIndex(start);
     if(parseResult->IsOptionPresent("--LocMapSpacing"))
       {
       double defScalarSpacing = parseResult->GetParameterFloat("--LocMapSpacing");
       std::cout<<"Generating coarse deformation field (spacing="<<defScalarSpacing<<")"<<std::endl;
       XsImageType::SpacingType defSpacing;
-      XsImageType::SizeType defSize;
 
       defSpacing[0] = defScalarSpacing;
       defSpacing[1] = defScalarSpacing;
-
-      for(unsigned int dim = 0; dim < XsImageType::ImageDimension;++dim)
-        {
-        defSize[dim] = static_cast<unsigned long>(size[dim]*vcl_abs(spacing[dim]/defSpacing[dim]));
-        }
-      std::cout<<"Deformation field spacing: "<<defSpacing<<std::endl;
-      std::cout<<"Deformation field size: "<<defSize<<std::endl;
-      dfGenerator->SetOutputSpacing(defSpacing);
-      dfGenerator->SetOutputSize(defSize);
+      
+      resampler->SetDeformationFieldSpacing(defSpacing);
       }
     else
       {
-      dfGenerator->SetOutputSpacing(spacing);
-      dfGenerator->SetOutputSize(size);
+      resampler->SetDeformationFieldSpacing(spacing);
       }
-    dfGenerator->SetTransform(rsTransform);
-
-    WarpFilterType::Pointer warper = WarpFilterType::New();
-
-    InterpolatorType::Pointer interpolator = InterpolatorType::New();
-    interpolator->SetInputImage(xsreader->GetOutput());
-    warper->SetInterpolator(interpolator);
-
+    
     XsImageType::PixelType defaultValue;
-    itk::PixelBuilder<XsImageType::PixelType>::Zero(defaultValue,xsreader->GetOutput()->GetNumberOfComponentsPerPixel());
+    itk::PixelBuilder<XsImageType::PixelType>::Zero(defaultValue,
+                                                    xsreader->GetOutput()->GetNumberOfComponentsPerPixel());
 
-    warper->SetInput(xsreader->GetOutput());
-    warper->SetDeformationField(dfGenerator->GetOutput());
-    warper->SetOutputOrigin(origin);
-    warper->SetOutputSpacing(spacing);
-    warper->SetOutputSize(size);
-    warper->SetOutputStartIndex(start);
-    warper->SetEdgePaddingValue(defaultValue);
+    resampler->SetInput(xsreader->GetOutput());
+    resampler->SetOutputOrigin(origin);
+    resampler->SetOutputSpacing(spacing);
+    resampler->SetOutputSize(size);
+    resampler->SetOutputStartIndex(start);
+    resampler->SetOutputKeywordList(preader->GetOutput()->GetImageKeywordlist());
+    resampler->SetEdgePaddingValue(defaultValue);
 
     FusionFilterType::Pointer  fusionFilter = FusionFilterType::New();
     fusionFilter->SetPanInput(preader->GetOutput());
-    fusionFilter->SetXsInput(warper->GetOutput());
+    fusionFilter->SetXsInput(resampler->GetOutput());
     fusionFilter->GetOutput()->UpdateOutputInformation();
-
-    std::cout<<"Output size: "<<fusionFilter->GetOutput()->GetLargestPossibleRegion().GetSize()<<std::endl;
-
-
+    
     WriterType::Pointer writer = WriterType::New();
     writer->SetFileName(parseResult->GetOutputImage());
     writer->SetInput(fusionFilter->GetOutput());
-- 
GitLab