From 9e6d8f098603f4c782b3a4c98626e8cd2435d03b Mon Sep 17 00:00:00 2001 From: remi cresson <remi.cresson@teledetection.fr> Date: Mon, 4 Sep 2017 09:36:22 +0000 Subject: [PATCH] ENH: Reproject input vector data bounding box --- app/otbExtractGeom.cxx | 148 ++++++++++++++++++++++------------------- 1 file changed, 80 insertions(+), 68 deletions(-) diff --git a/app/otbExtractGeom.cxx b/app/otbExtractGeom.cxx index 4abbbc9..98a3df9 100644 --- a/app/otbExtractGeom.cxx +++ b/app/otbExtractGeom.cxx @@ -29,6 +29,9 @@ // ogr #include "otbOGR.h" +// Projection +#include "otbGenericRSTransform.h" + using namespace std; namespace otb @@ -41,10 +44,10 @@ class ExtractGeom : public Application { public: /** Standard class typedefs. */ - typedef ExtractGeom Self; - typedef Application Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; + typedef ExtractGeom Self; + typedef Application Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; /** Standard macro */ itkNewMacro(Self); @@ -52,19 +55,24 @@ public: /** Filters */ typedef otb::MultiChannelExtractROI<FloatVectorImageType::InternalPixelType, - FloatVectorImageType::InternalPixelType> ExtractFilterType; + FloatVectorImageType::InternalPixelType> ExtractFilterType; /** mask */ - typedef bool TMaskPixelValueType; - typedef otb::Image<TMaskPixelValueType, 2> MaskImageType; + typedef bool TMaskPixelValueType; + typedef otb::Image<TMaskPixelValueType, 2> MaskImageType; typedef itk::MaskImageFilter<FloatVectorImageType, MaskImageType, - FloatVectorImageType> MaskFilterType; + FloatVectorImageType> MaskFilterType; /* vector data filters */ typedef otb::VectorDataIntoImageProjectionFilter<VectorDataType, - FloatVectorImageType> VectorDataReprojFilterType; + FloatVectorImageType> VectorDataReprojFilterType; typedef otb::VectorDataToLabelImageCustomFilter<VectorDataType, - MaskImageType> RasteriseFilterType; + MaskImageType> RasteriseFilterType; + + /** projection */ + typedef otb::GenericRSTransform<> RSTransformType; + typedef otb::GenericMapProjection<otb::TransformDirection::FORWARD> MapProjectionType; + typedef RSTransformType::InputPointType Point3DType; void DoInit() { @@ -83,8 +91,8 @@ public: AddParameter(ParameterType_InputImage, "in", "Input Image"); SetParameterDescription("in"," Input image."); - AddParameter(ParameterType_InputVectorData, "shp", "Input Shapefile" ); - SetParameterDescription("shp","Input Shapefile" ); + AddParameter(ParameterType_InputVectorData, "vec", "Input vector" ); + SetParameterDescription("vec","Input vector" ); AddParameter(ParameterType_OutputImage, "out", "Output image"); SetParameterDescription("out"," Output image."); @@ -92,13 +100,10 @@ public: AddRAMParameter(); // Doc example parameter settings - SetDocExampleParameterValue("shp", "vecteur.shp"); + SetDocExampleParameterValue("vec", "vecteur.shp"); SetDocExampleParameterValue("in", "QB_Toulouse_Ortho_XS.tif"); SetDocExampleParameterValue("out", "QB_Toulouse_Ortho_XS_decoup.tif uint16"); - - - } void DoUpdateParameters() @@ -111,55 +116,66 @@ public: // Get inputs FloatVectorImageType::Pointer xs = GetParameterImage("in"); - VectorDataType* shp = GetParameterVectorData("shp"); + VectorDataType* shp = GetParameterVectorData("vec"); /* Reproject vector data */ - vdReproj = VectorDataReprojFilterType::New(); - vdReproj->SetInputVectorData(shp); - vdReproj->SetInputImage(xs); - vdReproj->Update(); + m_VectorDataReprojectionFilter = VectorDataReprojFilterType::New(); + m_VectorDataReprojectionFilter->SetInputVectorData(shp); + m_VectorDataReprojectionFilter->SetInputImage(xs); + m_VectorDataReprojectionFilter->Update(); /* Get VectorData bounding box */ OGREnvelope env; otb::ogr::DataSource::Pointer ogrDS; - ogrDS = otb::ogr::DataSource::New(GetParameterString("shp") , + ogrDS = otb::ogr::DataSource::New(GetParameterString("vec") , otb::ogr::DataSource::Modes::Read); - double ulx, uly, lrx, lry; + itk::Point<double, 2> ulp_in, lrp_in; bool extentAvailable = true; std::string inputProjectionRef = ""; + // First try to get the extent in the metadata try { - inputProjectionRef = ogrDS->GetGlobalExtent(ulx,uly,lrx,lry); + inputProjectionRef = ogrDS->GetGlobalExtent(ulp_in[0],ulp_in[1],lrp_in[0],lrp_in[1]); } catch(const itk::ExceptionObject&) { - extentAvailable = false; + extentAvailable = false; } // If no extent available force the computation of the extent if (!extentAvailable) { - try - { - inputProjectionRef = ogrDS->GetGlobalExtent(ulx,uly,lrx,lry,true); - extentAvailable = true; + try + { + inputProjectionRef = ogrDS->GetGlobalExtent(ulp_in[0],ulp_in[1],lrp_in[0],lrp_in[1],true); + extentAvailable = true; + } + catch(itk::ExceptionObject & err) + { + extentAvailable = false; + + otbAppLogFATAL(<<"Failed to retrieve the spatial extent of the dataset " + "in force mode. The spatial extent is mandatory when " + "orx, ory, spx and spy parameters are not set, consider " + "setting them. Error from library: "<<err.GetDescription()); + } } - catch(itk::ExceptionObject & err) - { - extentAvailable = false; - otbAppLogFATAL(<<"Failed to retrieve the spatial extent of the dataset " - "in force mode. The spatial extent is mandatory when " - "orx, ory, spx and spy parameters are not set, consider " - "setting them. Error from library: "<<err.GetDescription()); - } - } + // Reproject region + RSTransformType::Pointer rsTransform = RSTransformType::New(); + rsTransform->SetInputProjectionRef(shp->GetProjectionRef()); + rsTransform->SetOutputKeywordList(xs->GetImageKeywordlist()); + rsTransform->SetOutputProjectionRef(xs->GetProjectionRef()); + rsTransform->InstantiateTransform(); + itk::Point<double, 2> ulp_out , lrp_out; + ulp_out = rsTransform->TransformPoint(ulp_in); + lrp_out = rsTransform->TransformPoint(lrp_in); otb::RegionComparator<FloatVectorImageType, FloatVectorImageType>::RSRegion::PointType rsRoiOrigin; - rsRoiOrigin[0] = ulx; - rsRoiOrigin[1] = uly; + rsRoiOrigin[0] = ulp_out[0]; + rsRoiOrigin[1] = ulp_out[1]; otb::RegionComparator<FloatVectorImageType, FloatVectorImageType>::RSRegion::SizeType rsRoiSize; - rsRoiSize[0] = lrx - ulx; - rsRoiSize[1] = lry - uly; + rsRoiSize[0] = lrp_out[ 0 ] - ulp_out[0]; + rsRoiSize[1] = lrp_out[ 1 ] - ulp_out[1]; otb::RegionComparator<FloatVectorImageType, FloatVectorImageType>::RSRegion rsRoi; rsRoi.SetOrigin(rsRoiOrigin); rsRoi.SetSize(rsRoiSize); @@ -167,46 +183,42 @@ public: /* Compute intersecting region */ otb::RegionComparator<FloatVectorImageType, FloatVectorImageType> comparator; comparator.SetImage1(xs); - FloatVectorImageType::RegionType roi = - comparator.RSRegionToImageRegion(rsRoi /*vdProperties->GetBoundingRegion()*/); + FloatVectorImageType::RegionType roi = comparator.RSRegionToImageRegion(rsRoi); roi.PadByRadius(1); // avoid extrapolation if (!roi.Crop(xs->GetLargestPossibleRegion())) { - otbAppLogFATAL( << " Input VectorData is outside image !" ); - return; + otbAppLogFATAL( << " Input VectorData is outside image !" ); + return; } /* Rasterize vector data */ - rasterizer = RasteriseFilterType::New(); - rasterizer->AddVectorData(vdReproj->GetOutput()); - rasterizer->SetOutputOrigin(xs->GetOrigin()); - rasterizer->SetOutputSpacing(xs->GetSpacing()); - rasterizer->SetOutputSize( - xs->GetLargestPossibleRegion().GetSize()); - rasterizer->SetBurnMaxValueMode(true); - rasterizer->SetOutputProjectionRef( - xs->GetProjectionRef()); - rasterizer->Update(); + m_RasterizeFilter = RasteriseFilterType::New(); + m_RasterizeFilter->AddVectorData(m_VectorDataReprojectionFilter->GetOutput()); + m_RasterizeFilter->SetOutputOrigin(xs->GetOrigin()); + m_RasterizeFilter->SetOutputSpacing(xs->GetSpacing()); + m_RasterizeFilter->SetOutputSize(xs->GetLargestPossibleRegion().GetSize()); + m_RasterizeFilter->SetBurnMaxValueMode(true); + m_RasterizeFilter->SetOutputProjectionRef(xs->GetProjectionRef()); + m_RasterizeFilter->Update(); /* Mask input image */ - maskFilter = MaskFilterType::New(); - maskFilter->SetInput(xs); - maskFilter->SetMaskImage(rasterizer->GetOutput()); + m_MaskFilter = MaskFilterType::New(); + m_MaskFilter->SetInput(xs); + m_MaskFilter->SetMaskImage(m_RasterizeFilter->GetOutput()); /* Extract ROI */ - m_Filter = ExtractFilterType::New(); - m_Filter->SetInput(maskFilter->GetOutput()); - m_Filter->SetExtractionRegion(roi); - - SetParameterOutputImage("out", m_Filter->GetOutput()); + m_ExtractFilter = ExtractFilterType::New(); + m_ExtractFilter->SetInput(m_MaskFilter->GetOutput()); + m_ExtractFilter->SetExtractionRegion(roi); + SetParameterOutputImage("out", m_ExtractFilter->GetOutput()); } - ExtractFilterType::Pointer m_Filter; - VectorDataReprojFilterType::Pointer vdReproj; - RasteriseFilterType::Pointer rasterizer; - MaskFilterType::Pointer maskFilter; + ExtractFilterType::Pointer m_ExtractFilter; + VectorDataReprojFilterType::Pointer m_VectorDataReprojectionFilter; + RasteriseFilterType::Pointer m_RasterizeFilter; + MaskFilterType::Pointer m_MaskFilter; }; } -- GitLab