ENH: add SAR sensor model tests

Testing/Code/Projections/CMakeLists.txt

@@ -237,6 +237,35 @@ ADD_TEST(prTvSensorModelErs ${PROJECTIONS_TESTS1}
          ${TEMP}/prTvSensorModelErs.txt
 )
 
+
+SET(SARSENSORMODEL_TYPES 
+  "RADARSAT1/GOMA/SCENE01/"
+  "RADARSAT2/ALTONA/Fine_Quad-Pol_Dataset/PK6621_DK406_FQ9_20080405_124900_HH_VV_HV_VH_SLC_Altona/product.xml"
+  "ENVISAT_ASAR_SCENE01/ASA_APS_1PNDPA20030821_7713.N1"
+  "PALSAR/200801280007/l1data/VOL-ALPSRP037120700-H1.1__A"
+  "TERRASARX/PANGKALANBUUN/IMAGEDATA/IMAGE_HH_SRA_stripFar_008.cos"
+)
+
+FOREACH( file ${SARSENSORMODEL_TYPES})
+  # Get the sensor name 
+  SET(sharp_regexp "([0-9A-Za-z_]*)[ ]*/[ ]*(.*)")
+  STRING(REGEX REPLACE "${sharp_regexp}" "\\1" sensor_name "${file}")
+  # Tests
+	ADD_TEST(prTvSensorModelGrid_${sensor_name} ${PROJECTIONS_TESTS1}
+	  --compare-n-images ${NOTOL} 2
+	         ${BASELINE}/prTvSensorModelGridDeltaX_${sensor_name}.tif
+	         ${TEMP}/prTvSensorModelGridDeltaX_${sensor_name}.tif
+	         ${BASELINE}/prTvSensorModelGridDeltaY_${sensor_name}.tif
+	         ${TEMP}/prTvSensorModelGridDeltaY_${sensor_name}.tif
+    	otbSensorModelGrid
+    		${LARGEINPUT}/${file}
+         	${TEMP}/prTvSensorModelGridDeltaX_${sensor_name}.tif
+         	${TEMP}/prTvSensorModelGridDeltaY_${sensor_name}.tif
+         	10 10
+		)
+ENDFOREACH( file ${SARSENSORMODEL_TYPES})
+
+
 ENDIF(OTB_DATA_USE_LARGEINPUT)
 
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -667,7 +696,9 @@ otbCreateProjectionWithOTB.cxx
 otbCreateInverseForwardSensorModel.cxx
 otbRegionProjectionResampler.cxx
 otbSensorModel.cxx
+otbSensorModelGrid.cxx
 )
+
 SET(Projections_SRCS2
 otbProjectionsTests2.cxx
 otbOrthoRectificationFilterNew.cxx

Testing/Code/Projections/otbProjectionsTests1.cxx

@@ -37,4 +37,5 @@ void RegisterTests()
   REGISTER_TEST(otbCreateInverseForwardSensorModel);
   REGISTER_TEST(otbRegionProjectionResampler);
   REGISTER_TEST(otbSensorModel);
+  REGISTER_TEST(otbSensorModelGrid);
 }

Testing/Code/Projections/otbSensorModelGrid.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.
+
+=========================================================================*/
+
+#if defined(_MSC_VER)
+#pragma warning ( disable : 4786 )
+#endif
+
+#include <iostream>
+
+#include "otbImage.h"
+#include "otbVectorImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "otbForwardSensorModel.h"
+#include "otbInverseSensorModel.h"
+
+int otbSensorModelGrid(int argc, char* argv[])
+{
+  if (argc != 6)
+    {
+    std::cout << argv[0] << " <input filename> <output filename in x> <output filename in y> <size_x> <size_y>" << std::endl;
+
+    return EXIT_FAILURE;
+    }
+
+  char * filename = argv[1];
+  char * outFilenameX = argv[2];
+  char * outFilenameY = argv[3];
+  unsigned int size_x((unsigned int)atoi(argv[4]));
+  unsigned int size_y((unsigned int)atoi(argv[5]));
+
+  typedef otb::VectorImage<double, 2>           VectorImageType;
+  typedef otb::Image<double, 2>                 ImageType;
+  typedef otb::ImageFileReader<VectorImageType> ReaderType;
+  typedef otb::ImageFileWriter<ImageType> WriterType;
+  typedef VectorImageType::SizeType             SizeType;
+
+  ReaderType::Pointer reader = ReaderType::New();
+  reader->SetFileName(filename);
+  reader->UpdateOutputInformation();
+
+  WriterType::Pointer writeErrorX = WriterType::New();
+  writeErrorX->SetFileName(outFilenameX);
+
+  WriterType::Pointer writeErrorY = WriterType::New();
+  writeErrorY->SetFileName(outFilenameY);
+
+  SizeType sizeIn = reader->GetOutput()->GetLargestPossibleRegion().GetSize();
+  SizeType sizeOut;
+  sizeOut[0] = size_x;
+  sizeOut[1] = size_y;
+
+  ImageType::RegionType region;
+  region.SetIndex(0,0);
+  region.SetSize(sizeOut);
+
+  ImageType::Pointer outputErrorX = ImageType::New();
+  outputErrorX->SetRegions(region);
+  outputErrorX->Allocate();
+
+  ImageType::Pointer outputErrorY = ImageType::New();
+  outputErrorY->SetRegions(region);
+  outputErrorY->Allocate();
+
+  double averageElevation = 16.19688987731934;
+
+  typedef otb::ForwardSensorModel<double> ForwardSensorModelType;
+  ForwardSensorModelType::Pointer forwardSensorModel = ForwardSensorModelType::New();
+  forwardSensorModel->SetImageGeometry(reader->GetOutput()->GetImageKeywordlist());
+  forwardSensorModel->SetAverageElevation(averageElevation);
+
+  typedef otb::InverseSensorModel<double> InverseSensorModelType;
+  InverseSensorModelType::Pointer inverseSensorModel = InverseSensorModelType::New();
+  inverseSensorModel->SetImageGeometry(reader->GetOutput()->GetImageKeywordlist());
+  inverseSensorModel->SetAverageElevation(averageElevation);
+
+  double deltaX = static_cast<double>(sizeIn[0]) / static_cast<double>(size_x);
+  double deltaY = static_cast<double>(sizeIn[1]) / static_cast<double>(size_y);
+
+  for(unsigned int i = 0 ; i < size_x; ++i)
+    {
+      for(unsigned int j = 0 ; j < size_y; ++j)
+        {
+          itk::Point<double, 2> imagePoint;
+          imagePoint[0] = deltaX * i;
+          imagePoint[1] = deltaY * j;
+
+          itk::Point<double, 2> geoPoint;
+          geoPoint = forwardSensorModel->TransformPoint(imagePoint);
+
+          itk::Point<double, 2> reversedImagePoint;
+          reversedImagePoint = inverseSensorModel->TransformPoint(geoPoint);
+
+          double errorX = imagePoint[0] - reversedImagePoint[0];
+          double errorY = imagePoint[1] - reversedImagePoint[1];
+
+          ImageType::IndexType index;
+          index[0]= i;
+          index[1]= j;
+
+          outputErrorX->SetPixel(index, errorX);
+          outputErrorY->SetPixel(index, errorY);
+        }
+    }
+
+  writeErrorX->SetInput(outputErrorX);
+  writeErrorX->Update();
+
+  writeErrorY->SetInput(outputErrorY);
+  writeErrorY->Update();
+
+  return EXIT_SUCCESS;
+}