Skip to content
Snippets Groups Projects
Commit bfb10a47 authored by Otmane Lahlou's avatar Otmane Lahlou
Browse files

ENH : use the new orhtorectification framework

parent bbc25874
Branches
Tags
No related merge requests found
......@@ -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());
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment