Commit 9e6d8f09 authored by remi cresson's avatar remi cresson
Browse files

ENH: Reproject input vector data bounding box

No related merge requests found
Showing with 80 additions and 68 deletions
+80 -68
...@@ -29,6 +29,9 @@ ...@@ -29,6 +29,9 @@
// ogr // ogr
#include "otbOGR.h" #include "otbOGR.h"
// Projection
#include "otbGenericRSTransform.h"
using namespace std; using namespace std;
namespace otb namespace otb
...@@ -41,10 +44,10 @@ class ExtractGeom : public Application ...@@ -41,10 +44,10 @@ class ExtractGeom : public Application
{ {
public: public:
/** Standard class typedefs. */ /** Standard class typedefs. */
typedef ExtractGeom Self; typedef ExtractGeom Self;
typedef Application Superclass; typedef Application Superclass;
typedef itk::SmartPointer<Self> Pointer; typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer; typedef itk::SmartPointer<const Self> ConstPointer;
/** Standard macro */ /** Standard macro */
itkNewMacro(Self); itkNewMacro(Self);
...@@ -52,19 +55,24 @@ public: ...@@ -52,19 +55,24 @@ public:
/** Filters */ /** Filters */
typedef otb::MultiChannelExtractROI<FloatVectorImageType::InternalPixelType, typedef otb::MultiChannelExtractROI<FloatVectorImageType::InternalPixelType,
FloatVectorImageType::InternalPixelType> ExtractFilterType; FloatVectorImageType::InternalPixelType> ExtractFilterType;
/** mask */ /** mask */
typedef bool TMaskPixelValueType; typedef bool TMaskPixelValueType;
typedef otb::Image<TMaskPixelValueType, 2> MaskImageType; typedef otb::Image<TMaskPixelValueType, 2> MaskImageType;
typedef itk::MaskImageFilter<FloatVectorImageType, MaskImageType, typedef itk::MaskImageFilter<FloatVectorImageType, MaskImageType,
FloatVectorImageType> MaskFilterType; FloatVectorImageType> MaskFilterType;
/* vector data filters */ /* vector data filters */
typedef otb::VectorDataIntoImageProjectionFilter<VectorDataType, typedef otb::VectorDataIntoImageProjectionFilter<VectorDataType,
FloatVectorImageType> VectorDataReprojFilterType; FloatVectorImageType> VectorDataReprojFilterType;
typedef otb::VectorDataToLabelImageCustomFilter<VectorDataType, 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() void DoInit()
{ {
...@@ -83,8 +91,8 @@ public: ...@@ -83,8 +91,8 @@ public:
AddParameter(ParameterType_InputImage, "in", "Input Image"); AddParameter(ParameterType_InputImage, "in", "Input Image");
SetParameterDescription("in"," Input image."); SetParameterDescription("in"," Input image.");
AddParameter(ParameterType_InputVectorData, "shp", "Input Shapefile" ); AddParameter(ParameterType_InputVectorData, "vec", "Input vector" );
SetParameterDescription("shp","Input Shapefile" ); SetParameterDescription("vec","Input vector" );
AddParameter(ParameterType_OutputImage, "out", "Output image"); AddParameter(ParameterType_OutputImage, "out", "Output image");
SetParameterDescription("out"," Output image."); SetParameterDescription("out"," Output image.");
...@@ -92,13 +100,10 @@ public: ...@@ -92,13 +100,10 @@ public:
AddRAMParameter(); AddRAMParameter();
// Doc example parameter settings // Doc example parameter settings
SetDocExampleParameterValue("shp", "vecteur.shp"); SetDocExampleParameterValue("vec", "vecteur.shp");
SetDocExampleParameterValue("in", "QB_Toulouse_Ortho_XS.tif"); SetDocExampleParameterValue("in", "QB_Toulouse_Ortho_XS.tif");
SetDocExampleParameterValue("out", "QB_Toulouse_Ortho_XS_decoup.tif uint16"); SetDocExampleParameterValue("out", "QB_Toulouse_Ortho_XS_decoup.tif uint16");
} }
void DoUpdateParameters() void DoUpdateParameters()
...@@ -111,55 +116,66 @@ public: ...@@ -111,55 +116,66 @@ public:
// Get inputs // Get inputs
FloatVectorImageType::Pointer xs = GetParameterImage("in"); FloatVectorImageType::Pointer xs = GetParameterImage("in");
VectorDataType* shp = GetParameterVectorData("shp"); VectorDataType* shp = GetParameterVectorData("vec");
/* Reproject vector data */ /* Reproject vector data */
vdReproj = VectorDataReprojFilterType::New(); m_VectorDataReprojectionFilter = VectorDataReprojFilterType::New();
vdReproj->SetInputVectorData(shp); m_VectorDataReprojectionFilter->SetInputVectorData(shp);
vdReproj->SetInputImage(xs); m_VectorDataReprojectionFilter->SetInputImage(xs);
vdReproj->Update(); m_VectorDataReprojectionFilter->Update();
/* Get VectorData bounding box */ /* Get VectorData bounding box */
OGREnvelope env; OGREnvelope env;
otb::ogr::DataSource::Pointer ogrDS; otb::ogr::DataSource::Pointer ogrDS;
ogrDS = otb::ogr::DataSource::New(GetParameterString("shp") , ogrDS = otb::ogr::DataSource::New(GetParameterString("vec") ,
otb::ogr::DataSource::Modes::Read); otb::ogr::DataSource::Modes::Read);
double ulx, uly, lrx, lry; itk::Point<double, 2> ulp_in, lrp_in;
bool extentAvailable = true; bool extentAvailable = true;
std::string inputProjectionRef = ""; std::string inputProjectionRef = "";
// First try to get the extent in the metadata // First try to get the extent in the metadata
try 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&) catch(const itk::ExceptionObject&)
{ {
extentAvailable = false; extentAvailable = false;
} }
// If no extent available force the computation of the extent // If no extent available force the computation of the extent
if (!extentAvailable) if (!extentAvailable)
{ {
try try
{ {
inputProjectionRef = ogrDS->GetGlobalExtent(ulx,uly,lrx,lry,true); inputProjectionRef = ogrDS->GetGlobalExtent(ulp_in[0],ulp_in[1],lrp_in[0],lrp_in[1],true);
extentAvailable = 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 " // Reproject region
"in force mode. The spatial extent is mandatory when " RSTransformType::Pointer rsTransform = RSTransformType::New();
"orx, ory, spx and spy parameters are not set, consider " rsTransform->SetInputProjectionRef(shp->GetProjectionRef());
"setting them. Error from library: "<<err.GetDescription()); 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; otb::RegionComparator<FloatVectorImageType, FloatVectorImageType>::RSRegion::PointType rsRoiOrigin;
rsRoiOrigin[0] = ulx; rsRoiOrigin[0] = ulp_out[0];
rsRoiOrigin[1] = uly; rsRoiOrigin[1] = ulp_out[1];
otb::RegionComparator<FloatVectorImageType, FloatVectorImageType>::RSRegion::SizeType rsRoiSize; otb::RegionComparator<FloatVectorImageType, FloatVectorImageType>::RSRegion::SizeType rsRoiSize;
rsRoiSize[0] = lrx - ulx; rsRoiSize[0] = lrp_out[ 0 ] - ulp_out[0];
rsRoiSize[1] = lry - uly; rsRoiSize[1] = lrp_out[ 1 ] - ulp_out[1];
otb::RegionComparator<FloatVectorImageType, FloatVectorImageType>::RSRegion rsRoi; otb::RegionComparator<FloatVectorImageType, FloatVectorImageType>::RSRegion rsRoi;
rsRoi.SetOrigin(rsRoiOrigin); rsRoi.SetOrigin(rsRoiOrigin);
rsRoi.SetSize(rsRoiSize); rsRoi.SetSize(rsRoiSize);
...@@ -167,46 +183,42 @@ public: ...@@ -167,46 +183,42 @@ public:
/* Compute intersecting region */ /* Compute intersecting region */
otb::RegionComparator<FloatVectorImageType, FloatVectorImageType> comparator; otb::RegionComparator<FloatVectorImageType, FloatVectorImageType> comparator;
comparator.SetImage1(xs); comparator.SetImage1(xs);
FloatVectorImageType::RegionType roi = FloatVectorImageType::RegionType roi = comparator.RSRegionToImageRegion(rsRoi);
comparator.RSRegionToImageRegion(rsRoi /*vdProperties->GetBoundingRegion()*/);
roi.PadByRadius(1); // avoid extrapolation roi.PadByRadius(1); // avoid extrapolation
if (!roi.Crop(xs->GetLargestPossibleRegion())) if (!roi.Crop(xs->GetLargestPossibleRegion()))
{ {
otbAppLogFATAL( << " Input VectorData is outside image !" ); otbAppLogFATAL( << " Input VectorData is outside image !" );
return; return;
} }
/* Rasterize vector data */ /* Rasterize vector data */
rasterizer = RasteriseFilterType::New(); m_RasterizeFilter = RasteriseFilterType::New();
rasterizer->AddVectorData(vdReproj->GetOutput()); m_RasterizeFilter->AddVectorData(m_VectorDataReprojectionFilter->GetOutput());
rasterizer->SetOutputOrigin(xs->GetOrigin()); m_RasterizeFilter->SetOutputOrigin(xs->GetOrigin());
rasterizer->SetOutputSpacing(xs->GetSpacing()); m_RasterizeFilter->SetOutputSpacing(xs->GetSpacing());
rasterizer->SetOutputSize( m_RasterizeFilter->SetOutputSize(xs->GetLargestPossibleRegion().GetSize());
xs->GetLargestPossibleRegion().GetSize()); m_RasterizeFilter->SetBurnMaxValueMode(true);
rasterizer->SetBurnMaxValueMode(true); m_RasterizeFilter->SetOutputProjectionRef(xs->GetProjectionRef());
rasterizer->SetOutputProjectionRef( m_RasterizeFilter->Update();
xs->GetProjectionRef());
rasterizer->Update();
/* Mask input image */ /* Mask input image */
maskFilter = MaskFilterType::New(); m_MaskFilter = MaskFilterType::New();
maskFilter->SetInput(xs); m_MaskFilter->SetInput(xs);
maskFilter->SetMaskImage(rasterizer->GetOutput()); m_MaskFilter->SetMaskImage(m_RasterizeFilter->GetOutput());
/* Extract ROI */ /* Extract ROI */
m_Filter = ExtractFilterType::New(); m_ExtractFilter = ExtractFilterType::New();
m_Filter->SetInput(maskFilter->GetOutput()); m_ExtractFilter->SetInput(m_MaskFilter->GetOutput());
m_Filter->SetExtractionRegion(roi); m_ExtractFilter->SetExtractionRegion(roi);
SetParameterOutputImage("out", m_Filter->GetOutput());
SetParameterOutputImage("out", m_ExtractFilter->GetOutput());
} }
ExtractFilterType::Pointer m_Filter; ExtractFilterType::Pointer m_ExtractFilter;
VectorDataReprojFilterType::Pointer vdReproj; VectorDataReprojFilterType::Pointer m_VectorDataReprojectionFilter;
RasteriseFilterType::Pointer rasterizer; RasteriseFilterType::Pointer m_RasterizeFilter;
MaskFilterType::Pointer maskFilter; MaskFilterType::Pointer m_MaskFilter;
}; };
} }
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment