An error occurred while loading the file. Please try again.
-
remi cresson authored768e7862
/*========================================================================= Copyright (c) Remi Cresson (IRSTEA). All rights reserved. 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 __gtbVectorDataToLabelImageCustomFilter_txx#define __gtbVectorDataToLabelImageCustomFilter_txx#include "otbVectorDataToLabelImageCustomFilter.h"#include "otbOGRIOHelper.h"#include "otbGdalDataTypeBridge.h"#include "otbMacro.h"#include "gdal_alg.h"#include "ogr_srs_api.h"#include "itkImageRegionIterator.h"namespace otb{template<class TVectorData, class TOutputImage>VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>::VectorDataToLabelImageCustomFilter() : m_OGRDataSourcePointer(0), m_BurnAttribute("FID"), m_BurnMaxValueMode(false) { 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 available m_DefaultBurnValue = 1.; }template<class TVectorData, class TOutputImage>voidVectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>::AddVectorData(const VectorDataType* vd) { this->itk::ProcessObject::PushBackInput( vd ); }template <class TVectorData, class TOutputImage>const typename VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>::VectorDataType *VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>::GetInput(unsigned int idx) { return static_cast<const TVectorData *> (this->itk::ProcessObject::GetInput(idx)); }template <class TVectorData, class TOutputImage>voidVectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>::SetOutputSpacing(const OutputSpacingType& spacing) { if (this->m_OutputSpacing != spacing) {7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
this->m_OutputSpacing = spacing;
this->Modified();
}
}
template <class TVectorData, class TOutputImage>
void
VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
::SetOutputSpacing(const double spacing[2])
{
OutputSpacingType s(spacing);
this->SetOutputSpacing(s);
}
template <class TVectorData, class TOutputImage>
void
VectorDataToLabelImageCustomFilter<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
VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
::SetOutputOrigin(const double origin[2])
{
OutputOriginType p(origin);
this->SetOutputOrigin(p);
}
template <class TVectorData, class TOutputImage>
void
VectorDataToLabelImageCustomFilter<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
VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
::SetOutputParametersFromImage(const ImageBaseType * src)
{
this->SetOutputOrigin ( src->GetOrigin() );
this->SetOutputSignedSpacing ( src->GetSignedSpacing() );
this->SetOutputSize ( src->GetLargestPossibleRegion().GetSize() );
otb::ImageMetadataInterfaceBase::Pointer imi = otb::ImageMetadataInterfaceFactory::CreateIMI(src->GetMetaDataDictionary());
this->SetOutputProjectionRef(imi->GetProjectionRef());
}
template<class TVectorData, class TOutputImage>
void
VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
::GenerateOutputInformation()
{
// get pointer to the output
OutputImagePointer outputPtr = this->GetOutput();
if (!outputPtr)
{
return;
}
// Set the size of the output region
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
typename TOutputImage::RegionType outputLargestPossibleRegion;
outputLargestPossibleRegion.SetSize(m_OutputSize);
//outputLargestPossibleRegion.SetIndex(m_OutputStartIndex);
outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
// Set spacing and origin
outputPtr->SetSignedSpacing(m_OutputSpacing);
outputPtr->SetOrigin(m_OutputOrigin);
itk::MetaDataDictionary& dict = outputPtr->GetMetaDataDictionary();
itk::EncapsulateMetaData<std::string> (dict, otb::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 inputIdx = 0; inputIdx < this->GetNumberOfInputs(); ++inputIdx)
{
const VectorDataType* vd = dynamic_cast< const VectorDataType*>(this->itk::ProcessObject::GetInput(inputIdx));
// 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;
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
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 (m_BurnMaxValueMode)
{
m_FullBurnValues.push_back(static_cast<double>(itk::NumericTraits<OutputImageInternalPixelType>::max()));
}
else 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
VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>::GenerateData()
{
// Call Superclass GenerateData
this->AllocateOutputs();
itk::ImageRegionIterator<TOutputImage> it (this->GetOutput(), this->GetOutput()->GetBufferedRegion());
for (it.GoToBegin(); !it.IsAtEnd(); ++it)
{
it.Set(0);
}
// 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(otb::GdalDataTypeBridge::GetGDALDataType<OutputImageInternalPixelType>()) << ","
<< "PIXELOFFSET=" << sizeof(OutputImageInternalPixelType) * nbBands << ","
<< "LINEOFFSET=" << sizeof(OutputImageInternalPixelType)*nbBands*bufferedRegion.GetSize()[0] << ","
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335
<< "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()->GetSignedSpacing()[0];
geoTransform[5] = this->GetOutput()->GetSignedSpacing()[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
VectorDataToLabelImageCustomFilter<TVectorData, TOutputImage>
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
} // end namespace otb
#endif