diff --git a/Code/BasicFilters/otbVectorDataToLabelImageFilter.h b/Code/BasicFilters/otbVectorDataToLabelImageFilter.h new file mode 100644 index 0000000000000000000000000000000000000000..4013a5edbdc366058d5b03f48d92160e60c74866 --- /dev/null +++ b/Code/BasicFilters/otbVectorDataToLabelImageFilter.h @@ -0,0 +1,178 @@ +/*========================================================================= + + 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. + +=========================================================================*/ +#ifndef __otbVectorDataToLabelImageFilter_h +#define __otbVectorDataToLabelImageFilter_h + +#include "itkImageToImageFilter.h" +#include "itkImageSource.h" +#include "otbMacro.h" + +#include "otbVectorData.h" + +#include "gdal.h" +#include "ogr_api.h" +#include <ogrsf_frmts.h> + +namespace otb { + +/** \class VectorDataToLabelImageFilter + * \brief Burn geometries from the specified VectorData into raster + * + * This class handles burning several input VectorDatas into the + * output raster. The burn values are extracted from a field set by + * the user.If no burning field is set, the "FID" is used by default. + * + * Setting the output raster informations can be done in two ways by: + * - Setting the Origin/Size/Spacing of the output image + * - Using an existing image as support via SetOutputParametersFromImage(ImageBase) + * + */ +template <class TVectorData, class TOutputImage > +class ITK_EXPORT VectorDataToLabelImageFilter : + public itk::ImageSource<TOutputImage> +{ +public: + /** Standard class typedefs */ + typedef VectorDataToLabelImageFilter Self; + typedef itk::ImageSource<TOutputImage> Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + /** Run-time type information (and related methods). */ + itkTypeMacro(VectorDataToLabelImageFilter, itk::ImageSource); + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + typedef TOutputImage OutputImageType; + typedef typename OutputImageType::Pointer OutputImagePointer; + typedef typename OutputImageType::SizeType OutputSizeType; + typedef typename OutputImageType::IndexType OutputIndexType; + typedef typename OutputImageType::SpacingType OutputSpacingType; + typedef typename OutputImageType::PointType OutputOriginType; + typedef typename OutputImageType::RegionType OutputImageRegionType; + typedef typename OutputImageType::PixelType OutputImagePixelType; + typedef typename OutputImageType::InternalPixelType OutputImageInternalPixelType; + + /** VectorData typedefs*/ + typedef TVectorData VectorDataType; + typedef typename VectorDataType::DataTreeType DataTreeType; + typedef typename DataTreeType::TreeNodeType InternalTreeNodeType; + typedef typename DataTreeType::Pointer DataTreePointerType; + typedef typename DataTreeType::ConstPointer DataTreeConstPointerType; + + typedef itk::ImageBase<OutputImageType::ImageDimension> ImageBaseType; + + /** Get Nth input VectorData */ + const VectorDataType* GetInput(unsigned int idx); + + /** Method for adding VectorData to rasterize */ + virtual void AddVectorData(const VectorDataType* vd); + + /** Set the size of the output image. */ + itkSetMacro(OutputSize, OutputSizeType); + + /** Get the size of the output image. */ + itkGetConstReferenceMacro(OutputSize, OutputSizeType); + + /** Set the origin of the output image. + * \sa GetOrigin() + */ + itkSetMacro(OutputOrigin, OutputOriginType); + virtual void SetOutputOrigin(const double origin[2]); + virtual void SetOutputOrigin(const float origin[2]); + + itkGetConstReferenceMacro(OutputOrigin, OutputOriginType); + + /** Set the spacing (size of a pixel) of the output image. + * \sa GetSpacing() + */ + virtual void SetOutputSpacing(const OutputSpacingType& spacing); + virtual void SetOutputSpacing(const double spacing[2]); + virtual void SetOutputSpacing(const float spacing[2]); + + /** Set/Get Output Projection Ref */ + itkSetStringMacro(OutputProjectionRef); + itkGetStringMacro(OutputProjectionRef); + + itkSetStringMacro(BurnAttribute); + itkGetStringMacro(BurnAttribute); + + /** Useful to set the output parameters from an existing image*/ + void SetOutputParametersFromImage(const ImageBaseType * image); + +protected: + virtual void GenerateData(); + + VectorDataToLabelImageFilter(); + virtual ~VectorDataToLabelImageFilter() + { + // Destroy the geometries stored + for (unsigned int idx = 0; idx < m_SrcDataSetGeometries.size(); ++idx) + { + OGR_G_DestroyGeometry(m_SrcDataSetGeometries[idx]); + } + + if (m_OGRDataSourcePointer != NULL) + { + OGRDataSource::DestroyDataSource(m_OGRDataSourcePointer); + } + // Clean up the environment + OGRCleanupAll(); + GDALDestroyDriverManager(); + } + + virtual void GenerateOutputInformation(); + + void PrintSelf(std::ostream& os, itk::Indent indent) const; + +private: + VectorDataToLabelImageFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + OGRDataSource* m_OGRDataSourcePointer; + + // Vector Of OGRGeometyH + std::vector< OGRGeometryH > m_SrcDataSetGeometries; + + std::vector<double> m_BurnValues; + std::vector<double> m_FullBurnValues; + std::vector<int> m_BandsToBurn; + + // Field used to extract the burn value + std::string m_BurnAttribute; + + // Default burn value + double m_DefaultBurnValue; + + // Output params + std::string m_OutputProjectionRef; + OutputSpacingType m_OutputSpacing; + OutputOriginType m_OutputOrigin; + OutputSizeType m_OutputSize; + OutputIndexType m_OutputStartIndex; +}; // end of class VectorDataToLabelImageFilter + +} // end of namespace otb + + +#ifndef OTB_MANUAL_INSTANTIATION +#include "otbVectorDataToLabelImageFilter.txx" +#endif + +#endif diff --git a/Code/BasicFilters/otbVectorDataToLabelImageFilter.txx b/Code/BasicFilters/otbVectorDataToLabelImageFilter.txx new file mode 100644 index 0000000000000000000000000000000000000000..0e0bfb7b3a9fe7732b65fe444b09c420bf9e1ef3 --- /dev/null +++ b/Code/BasicFilters/otbVectorDataToLabelImageFilter.txx @@ -0,0 +1,325 @@ +/*========================================================================= + + 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. + +=========================================================================*/ + +#include "otbVectorDataToLabelImageFilter.h" +#include "otbOGRIOHelper.h" +#include "otbGdalDataTypeBridge.h" +#include "otbMacro.h" + +#include "gdal_alg.h" +#include "ogr_srs_api.h" + +namespace otb +{ +template<class TVectorData, class TOutputImage> +VectorDataToLabelImageFilter<TVectorData, TOutputImage> +::VectorDataToLabelImageFilter() : m_BurnAttribute("FID") +{ + this->SetNumberOfRequiredInputs(1); + + // Output parameters initialization + m_OutputSpacing.Fill(1.0); + m_OutputSize.Fill(0); + m_OutputStartIndex.Fill(0); + + // This filter is intended to work with otb::Image + m_BandsToBurn.push_back(1); + + // Default burn value if no burnAttribute availabke + m_DefaultBurnValue = 1.; +} + +template<class TVectorData, class TOutputImage> +void +VectorDataToLabelImageFilter<TVectorData, TOutputImage> +::AddVectorData(const VectorDataType* vd) +{ + this->itk::ProcessObject::PushBackInput( vd ); +} + +template <class TVectorData, class TOutputImage> +const typename VectorDataToLabelImageFilter<TVectorData, TOutputImage>::VectorDataType * +VectorDataToLabelImageFilter<TVectorData, TOutputImage> +::GetInput(unsigned int idx) +{ + return static_cast<const TVectorData *> + (this->itk::ProcessObject::GetInput(idx)); +} + +template <class TVectorData, class TOutputImage> +void +VectorDataToLabelImageFilter<TVectorData, TOutputImage> +::SetOutputSpacing(const OutputSpacingType& spacing) +{ + if (this->m_OutputSpacing != spacing) + { + this->m_OutputSpacing = spacing; + this->Modified(); + } +} + +template <class TVectorData, class TOutputImage> +void +VectorDataToLabelImageFilter<TVectorData, TOutputImage> +::SetOutputSpacing(const double spacing[2]) +{ + OutputSpacingType s(spacing); + this->SetOutputSpacing(s); +} + +template <class TVectorData, class TOutputImage> +void +VectorDataToLabelImageFilter<TVectorData, TOutputImage> +::SetOutputSpacing(const float spacing[2]) +{ + itk::Vector<float, 2> sf(spacing); + OutputSpacingType s; + s.CastFrom(sf); + this->SetOutputSpacing(s); +} + +template <class TVectorData, class TOutputImage> +void +VectorDataToLabelImageFilter<TVectorData, TOutputImage> +::SetOutputOrigin(const double origin[2]) +{ + OutputOriginType p(origin); + this->SetOutputOrigin(p); +} + +template <class TVectorData, class TOutputImage> +void +VectorDataToLabelImageFilter<TVectorData, TOutputImage> +::SetOutputOrigin(const float origin[2]) +{ + itk::Point<float, 2> of(origin); + OutputOriginType p; + p.CastFrom(of); + this->SetOutputOrigin(p); +} + +template <class TVectorData, class TOutputImage> +void +VectorDataToLabelImageFilter<TVectorData, TOutputImage> +::SetOutputParametersFromImage(const ImageBaseType * image) +{ + const OutputImageType * src = dynamic_cast<const OutputImageType*>(image); + + this->SetOutputOrigin ( src->GetOrigin() ); + this->SetOutputSpacing ( src->GetSpacing() ); + //this->SetOutputStartIndex ( src->GetLargestPossibleRegion().GetIndex() ); + this->SetOutputSize ( src->GetLargestPossibleRegion().GetSize() ); + this->SetOutputProjectionRef(src->GetProjectionRef()); + //this->SetOutputKeywordList(src->GetImageKeywordlist()); +} + +template<class TVectorData, class TOutputImage> +void +VectorDataToLabelImageFilter<TVectorData, TOutputImage> +::GenerateOutputInformation() +{ + // get pointer to the output + OutputImagePointer outputPtr = this->GetOutput(); + if (!outputPtr) + { + return; + } + + // Set the size of the output region + typename TOutputImage::RegionType outputLargestPossibleRegion; + outputLargestPossibleRegion.SetSize(m_OutputSize); + //outputLargestPossibleRegion.SetIndex(m_OutputStartIndex); + outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion); + + // Set spacing and origin + outputPtr->SetSpacing(m_OutputSpacing); + outputPtr->SetOrigin(m_OutputOrigin); + + itk::MetaDataDictionary& dict = outputPtr->GetMetaDataDictionary(); + itk::EncapsulateMetaData<std::string> (dict, MetaDataKey::ProjectionRefKey, + static_cast<std::string>(this->GetOutputProjectionRef())); + + // Generate the OGRLayers from the input VectorDatas + // iteration begin from 1 cause the 0th input is a image + for (unsigned int idx = 0; idx < this->GetNumberOfInputs(); ++idx) + { + const VectorDataType* vd = dynamic_cast< const VectorDataType*>(this->itk::ProcessObject::GetInput(idx)); + + // Get the projection ref of the current VectorData + std::string projectionRefWkt = vd->GetProjectionRef(); + bool projectionInformationAvailable = !projectionRefWkt.empty(); + OGRSpatialReference * oSRS = NULL; + + if (projectionInformationAvailable) + { + oSRS = static_cast<OGRSpatialReference *>(OSRNewSpatialReference(projectionRefWkt.c_str())); + } + else + { + otbMsgDevMacro(<< "Projection information unavailable"); + } + + // Retrieving root node + DataTreeConstPointerType tree = vd->GetDataTree(); + + // Get the input tree root + InternalTreeNodeType * inputRoot = const_cast<InternalTreeNodeType *>(tree->GetRoot()); + + // Iterative method to build the layers from a VectorData + OGRLayer * ogrCurrentLayer = NULL; + std::vector<OGRLayer *> ogrLayerVector; + otb::OGRIOHelper::Pointer IOConversion = otb::OGRIOHelper::New(); + + // The method ConvertDataTreeNodeToOGRLayers create the + // OGRDataSource but don t release it. Destruction is done in the + // desctructor + m_OGRDataSourcePointer = NULL; + ogrLayerVector = IOConversion->ConvertDataTreeNodeToOGRLayers(inputRoot, + m_OGRDataSourcePointer, + ogrCurrentLayer, + oSRS); + + // From OGRLayer* to OGRGeometryH vector + for (unsigned int idx = 0; idx < ogrLayerVector.size(); ++idx) + { + // test if the layers contain a field m_BurnField; + int burnField = -1; + + if( !m_BurnAttribute.empty() ) + { + burnField = OGR_FD_GetFieldIndex( OGR_L_GetLayerDefn( (OGRLayerH)(ogrLayerVector[idx]) ), + m_BurnAttribute.c_str() ); + + // Get the geometries of the layer + OGRFeatureH hFeat; + OGR_L_ResetReading( (OGRLayerH)(ogrLayerVector[idx]) ); + while( ( hFeat = OGR_L_GetNextFeature( (OGRLayerH)(ogrLayerVector[idx]) )) != NULL ) + { + OGRGeometryH hGeom; + if( OGR_F_GetGeometryRef( hFeat ) == NULL ) + { + OGR_F_Destroy( hFeat ); + continue; + } + + hGeom = OGR_G_Clone( OGR_F_GetGeometryRef( hFeat ) ); + m_SrcDataSetGeometries.push_back( hGeom ); + + if (burnField == -1 ) + { + // TODO : if no burnAttribute available, warning or raise an exception?? + m_FullBurnValues.push_back(m_DefaultBurnValue++); + itkWarningMacro(<<"Failed to find attribute "<<m_BurnAttribute << " in layer " + << OGR_FD_GetName( OGR_L_GetLayerDefn( (OGRLayerH)(ogrLayerVector[idx]) )) + <<" .Setting burn value to default = " + << m_DefaultBurnValue); + } + else + { + m_FullBurnValues.push_back( OGR_F_GetFieldAsDouble( hFeat, burnField ) ); + } + + OGR_F_Destroy( hFeat ); + } + } + + // Destroy the oSRS + if (oSRS != NULL) + { + OSRRelease(oSRS); + } + } + } +} + +template<class TVectorData, class TOutputImage> +void +VectorDataToLabelImageFilter<TVectorData, TOutputImage>::GenerateData() +{ + // Call Superclass GenerateData + this->AllocateOutputs(); + + // Get the buffered region + OutputImageRegionType bufferedRegion = this->GetOutput()->GetBufferedRegion(); + + // nb bands + unsigned int nbBands = this->GetOutput()->GetNumberOfComponentsPerPixel(); + + // register drivers + GDALAllRegister(); + + std::ostringstream stream; + stream << "MEM:::" + << "DATAPOINTER=" << (unsigned long)(this->GetOutput()->GetBufferPointer()) << "," + << "PIXELS=" << bufferedRegion.GetSize()[0] << "," + << "LINES=" << bufferedRegion.GetSize()[1]<< "," + << "BANDS=" << nbBands << "," + << "DATATYPE=" << GDALGetDataTypeName(GdalDataTypeBridge::GetGDALDataType<OutputImageInternalPixelType>()) << "," + << "PIXELOFFSET=" << sizeof(OutputImageInternalPixelType) * nbBands << "," + << "LINEOFFSET=" << sizeof(OutputImageInternalPixelType)*nbBands*bufferedRegion.GetSize()[0] << "," + << "BANDOFFSET=" << sizeof(OutputImageInternalPixelType); + + GDALDatasetH dataset = GDALOpen(stream.str().c_str(), GA_Update); + + // Add the projection ref to the dataset + GDALSetProjection (dataset, this->GetOutput()->GetProjectionRef().c_str()); + + // add the geoTransform to the dataset + itk::VariableLengthVector<double> geoTransform(6); + + // Reporting origin and spacing of the buffered region + // the spacing is unchanged, the origin is relative to the buffered region + OutputIndexType bufferIndexOrigin = bufferedRegion.GetIndex(); + OutputOriginType bufferOrigin; + this->GetOutput()->TransformIndexToPhysicalPoint(bufferIndexOrigin, bufferOrigin); + geoTransform[0] = bufferOrigin[0]; + geoTransform[3] = bufferOrigin[1]; + geoTransform[1] = this->GetOutput()->GetSpacing()[0]; + geoTransform[5] = this->GetOutput()->GetSpacing()[1]; + + // FIXME: Here component 1 and 4 should be replaced by the orientation parameters + geoTransform[2] = 0.; + geoTransform[4] = 0.; + GDALSetGeoTransform(dataset,const_cast<double*>(geoTransform.GetDataPointer())); + + // Burn the geometries into the dataset + if (dataset != NULL) + { + GDALRasterizeGeometries( dataset, m_BandsToBurn.size(), + &(m_BandsToBurn[0]), + m_SrcDataSetGeometries.size(), + &(m_SrcDataSetGeometries[0]), + NULL, NULL, &(m_FullBurnValues[0]), + NULL, + GDALDummyProgress, NULL ); + + // release the dataset + GDALClose( dataset ); + } +} + +template<class TVectorData, class TOutputImage> +void +VectorDataToLabelImageFilter<TVectorData, TOutputImage> +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + Superclass::PrintSelf(os, indent); +} + +} // end namespace otb +