Commit 1267e4a3 authored by Cyrille Valladeau's avatar Cyrille Valladeau
Browse files

MRG

parents fa0d7e9f 40dfadb8
......@@ -73,7 +73,7 @@ public:
/** Size typedef support. A size is used to define region bounds. */
typedef itk::ContinuousIndex<Type> SizeType;
typedef itk::Size<2> StandardSizeType;
// typedef itk::Size<2> StandardSizeType;
/** ImageRegion typedef needed by the GetImageRegion() method */
typedef itk::ImageRegion<2> ImageRegionType;
......@@ -185,11 +185,11 @@ public:
m_Size = size;
}
void SetSize(const StandardSizeType &size)
{
m_Size[0] = size[0];
m_Size[1] = size[1];
}
// void SetSize(const StandardSizeType &size)
// {
// m_Size[0] = size[0];
// m_Size[1] = size[1];
// }
/** Get the size of the region. */
const SizeType& GetSize() const
......@@ -254,11 +254,13 @@ public:
{
for(unsigned int i=0; i<IndexType::IndexDimension; i++)
{
if( index[i] < m_Index[i] )
if(( index[i] < m_Index[i] )
&& (index[i] < m_Index[i] + m_Size[i]))
{
return false;
}
if( index[i] >= m_Index[i] + m_Size[i] )
if(( index[i] >= m_Index[i] )
&& (index[i] >= m_Index[i] + m_Size[i] ))
{
return false;
}
......
......@@ -283,6 +283,8 @@ VectorDataExtractROI<TVectorData>
std::string regionProjection = m_ROI.GetRegionProjection();
std::string inputVectorProjection = this->GetInput()->GetProjectionRef();
//FIXME: the string comparison is not sufficient to say that two
//projections are different
if(regionProjection == inputVectorProjection)
m_ProjectionNeeded = false;
else
......
This diff is collapsed.
......@@ -26,8 +26,9 @@
namespace otb
{
/** \class VectorDataStyle
* \brief Load style into a mapnik::Map
*
* \brief Load style into a mapnik::Map
*
* \sa VectorDataToImageFilter
*/
......@@ -48,9 +49,20 @@ class ITK_EXPORT VectorDataStyle : public itk::LightObject
void LoadOSMStyle(mapnik::Map& mapnikMap) const;
/** Get/Set the size of the scale factor. */
double GeScaleFactor() const
{
return m_ScaleFactor;
}
void SetScaleFactor( double scale )
{
m_ScaleFactor = scale;
}
protected:
/** Constructor */
VectorDataStyle() {};
VectorDataStyle():
m_ScaleFactor(1.0){};
/** Destructor */
virtual ~VectorDataStyle() {};
/**PrintSelf method */
......@@ -60,6 +72,10 @@ class ITK_EXPORT VectorDataStyle : public itk::LightObject
VectorDataStyle(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
//this parameter is used only in the case of sensor geometry
//to adjust the scale
double m_ScaleFactor;
}; // end class
} // end namespace otb
......
......@@ -133,6 +133,10 @@ template <class TVectorData, class TImage>
mapnik::Map m_Map;
//This factor is used to flip the data on the Y axis when using a
//sensor model geometry (where the Y coordinate increases top-down)
int m_SensorModelFlip;
}; // end class
} // end namespace otb
......
......@@ -21,6 +21,7 @@
#include "otbVectorDataToImageFilter.h"
#include "itkImageRegionIterator.h"
#include "itkImageIteratorWithIndex.h"
#include "otbVectorDataStyle.h"
#include <mapnik/datasource_cache.hpp>
......@@ -53,6 +54,7 @@ namespace otb
m_Size.Fill( 0 );
m_StartIndex.Fill( 0 );
m_Map = mapnik::Map();
m_SensorModelFlip = 1;
}
......@@ -218,12 +220,15 @@ namespace otb
itk::ExposeMetaData<std::string>(input->GetMetaDataDictionary(), MetaDataKey::ProjectionRefKey, vectorDataProjectionWKT);
std::cout << "WKT -> " << vectorDataProjectionWKT << std::endl;
std::string vectorDataProjectionProj4;
m_SensorModelFlip = 1;
if (vectorDataProjectionWKT == "")
{
//We assume that it is an image in sensor model geometry
//and tell mapnik that this is utm
//(with a resolution of 1m per unit)
vectorDataProjectionProj4 = "+proj=utm +zone=31 +ellps=WGS84";
m_SensorModelFlip = -1;
}
else
{
......@@ -232,6 +237,7 @@ namespace otb
oSRS.exportToProj4(&pszProj4);
vectorDataProjectionProj4 = pszProj4;
CPLFree(pszProj4);
m_SensorModelFlip = 1;
}
std::cout << "Proj.4 -> " << vectorDataProjectionProj4 << std::endl;
......@@ -254,10 +260,12 @@ namespace otb
m_Map.addLayer(lyr);
assert( (m_SensorModelFlip == 1)||(m_SensorModelFlip == -1) );
mapnik::Envelope<double> envelope(m_Origin[0],
m_Origin[1]+m_Spacing[1]*m_Size[1],
m_Origin[0]+m_Spacing[0]*m_Size[0],
m_Origin[1]);
m_SensorModelFlip*(m_Origin[1]+m_Spacing[1]*m_Size[1]),
m_Origin[0]+m_Spacing[0]*m_Size[0],
m_SensorModelFlip*m_Origin[1]);
// m_Map.zoomToBox(lyr.envelope());//FIXME: use the Origin/Spacing to calculate this
m_Map.zoomToBox(envelope);
// std::cout << "Envelope: " << lyr.envelope() << std::endl;
......@@ -276,22 +284,23 @@ namespace otb
// unsigned int size = 4*m_Size[0]*m_Size[1];//FIXME: lot of assumption here: RGBA, 2 dimensions
// memcpy(dst, src, size);
itk::ImageRegionIterator<ImageType> it(output, output->GetLargestPossibleRegion());
for (it.GoToBegin(); !it.IsAtEnd(); ++it)
{
itk::RGBAPixel<unsigned char> pix;
pix[0] = *src;
++src;
pix[1] = *src;
++src;
pix[2] = *src;
++src;
pix[3] = *src;
++src;
it.Set(pix);
}
// if (m_SensorModelFlip == false)
// {
itk::ImageRegionIterator<ImageType> it(output, output->GetLargestPossibleRegion());
for (it.GoToBegin(); !it.IsAtEnd(); ++it)
{
itk::RGBAPixel<unsigned char> pix;
pix[0] = *src;
++src;
pix[1] = *src;
++src;
pix[2] = *src;
++src;
pix[3] = *src;
++src;
it.Set(pix);
}
this->AfterThreadedGenerateData();
}
......@@ -333,14 +342,10 @@ namespace otb
ProcessNode((*it),mDatasource);
break;
}
// case FEATURE_POINT:
// {
// newDataNode->SetPoint(this->ReprojectPoint(dataNode->GetPoint()));
// newContainer = OutputInternalTreeNodeType::New();
// newContainer->Set(newDataNode);
// destination->AddChild(newContainer);
// break;
// }
case FEATURE_POINT:
{
itkExceptionMacro(<<"This type (FEATURE_POINT) is not handle (yet) by VectorDataToImageFilter(), please request for it");
}
case otb::FEATURE_LINE:
{
std::cout << std::setprecision(15);
......@@ -355,7 +360,7 @@ namespace otb
while (itVertex != dataNode->GetLine()->GetVertexList()->End())
{
// std::cout << itVertex.Value()[0] << ", " << itVertex.Value()[1] << std::endl;
line->line_to(itVertex.Value()[0],itVertex.Value()[1]);
line->line_to(itVertex.Value()[0],m_SensorModelFlip*itVertex.Value()[1]);
++itVertex;
}
......@@ -381,50 +386,28 @@ namespace otb
mDatasource->push(mfeature);
break;
}
// case FEATURE_POLYGON:
// {
// newDataNode->SetPolygonExteriorRing(this->ReprojectPolygon(dataNode->GetPolygonExteriorRing()));
// newDataNode->SetPolygonInteriorRings(this->ReprojectPolygonList(dataNode->GetPolygonInteriorRings()));
// newContainer = OutputInternalTreeNodeType::New();
// newContainer->Set(newDataNode);
// destination->AddChild(newContainer);
// break;
// }
// case FEATURE_MULTIPOINT:
// {
// newContainer = OutputInternalTreeNodeType::New();
// newContainer->Set(newDataNode);
// destination->AddChild(newContainer);
// ProcessNode((*it),newContainer);
// break;
// }
// case FEATURE_MULTILINE:
// {
// newContainer = OutputInternalTreeNodeType::New();
// newContainer->Set(newDataNode);
// destination->AddChild(newContainer);
// ProcessNode((*it),newContainer);
// break;
// }
// case FEATURE_MULTIPOLYGON:
// {
// newContainer = OutputInternalTreeNodeType::New();
// newContainer->Set(newDataNode);
// destination->AddChild(newContainer);
// ProcessNode((*it),newContainer);
// break;
// }
// case FEATURE_COLLECTION:
// {
// newContainer = OutputInternalTreeNodeType::New();
// newContainer->Set(newDataNode);
// destination->AddChild(newContainer);
// ProcessNode((*it),newContainer);
// break;
// }
case FEATURE_POLYGON:
{
itkExceptionMacro(<<"This type (FEATURE_POLYGON) is not handle (yet) by VectorDataToImageFilter(), please request for it");
}
case FEATURE_MULTIPOINT:
{
itkExceptionMacro(<<"This type (FEATURE_MULTIPOINT) is not handle (yet) by VectorDataToImageFilter(), please request for it");
}
case FEATURE_MULTILINE:
{
itkExceptionMacro(<<"This type (FEATURE_MULTILINE) is not handle (yet) by VectorDataToImageFilter(), please request for it");
}
case FEATURE_MULTIPOLYGON:
{
itkExceptionMacro(<<"This type (FEATURE_MULTIPOLYGON) is not handle (yet) by VectorDataToImageFilter(), please request for it");
}
case FEATURE_COLLECTION:
{
itkExceptionMacro(<<"This type (FEATURE_COLLECTION) is not handle (yet) by VectorDataToImageFilter(), please request for it");
}
}
}
}
......
......@@ -62,7 +62,7 @@ PersistentLineSegmentDetector<TInputImage, TPrecision>
}
}
template<class TInputImage, class TPrecision>
......@@ -113,7 +113,7 @@ PersistentLineSegmentDetector<TInputImage, TPrecision>
{
ImagePointer inputPtr = const_cast< TInputImage * >( this->GetInput() );
ImagePointer outputPtr = const_cast< TInputImage * >( this->GetOutput() );
inputPtr->SetRequestedRegion(outputRegionForThread);
inputPtr->PropagateRequestedRegion();
inputPtr->UpdateOutputData();
......@@ -173,12 +173,12 @@ ImageToLineSegmentVectorData<TInputImage, TPrecision>
{
itkExceptionMacro(<<"No LineSpatialObject to vectorize");
}
VectorDataPointerType vlines = VectorDataType::New();
DataNodePointerType document = DataNodeType::New();
DataNodePointerType folder = DataNodeType::New();
DataNodePointerType root = vlines->GetDataTree()->GetRoot()->Get();
document->SetNodeType(otb::DOCUMENT);
folder->SetNodeType(otb::FOLDER);
vlines->GetDataTree()->Add(document,root);
......@@ -193,7 +193,7 @@ ImageToLineSegmentVectorData<TInputImage, TPrecision>
std::vector<DoubleVertexType> vertexList, vertexNewList;
for(unsigned int i = 0; i<listSize; i++ )
{
{
typename LineSpatialObjectType::const_iterator it = this->GetLines()->GetNthElement(i)->begin();
typename LineSpatialObjectType::const_iterator itEnd = this->GetLines()->GetNthElement(i)->end();
......@@ -208,7 +208,7 @@ ImageToLineSegmentVectorData<TInputImage, TPrecision>
itPoints++;
p2[0] =(*itPoints).GetPosition()[0]; //Second Vertex
p2[1] =(*itPoints).GetPosition()[1];
bool go = false;
// horizontal line
if( vcl_abs(p1[1]-whereAmI[1])<m_ThreadDistanceThreshold && vcl_abs(p2[1]-whereAmI[1])<m_ThreadDistanceThreshold )
......@@ -287,23 +287,23 @@ ImageToLineSegmentVectorData<TInputImage, TPrecision>
{
l->AddVertex(p1);
l->AddVertex(p2);
DataNodePointerType node = DataNodeType::New();
node->SetNodeType(otb::FEATURE_LINE);
node->SetLine(l);
vlines->GetDataTree()->Add(node,folder);
}
++it;
} // End while(it != itEnd) loop
// write the false alarm in vertexList (those that don't have a continuation in the next thread)
for(unsigned int k=0; k<vertexList.size(); k++)
{
typename LineType::Pointer l = LineType::New();
l->AddVertex(vertexList[k][0]);
l->AddVertex(vertexList[k][1]);
DataNodePointerType node = DataNodeType::New();
node->SetNodeType(otb::FEATURE_LINE);
node->SetLine(l);
......
......@@ -103,14 +103,14 @@ DEMToOrthoImageGenerator<TDEMImage, TMapProjection>
DEMImage->TransformIndexToPhysicalPoint(currentindex, cartoPoint);
otbMsgDevMacro(<< "CartoPoint : (" << cartoPoint[0] << "," << cartoPoint[1] << ")");
// otbMsgDevMacro(<< "CartoPoint: (" << cartoPoint[0] << "," << cartoPoint[1] << ")");
geoPoint = m_MapProjection->TransformPoint(cartoPoint);
otbMsgDevMacro(<< "CartoPoint : (" << geoPoint[0] << "," << geoPoint[1] << ")");
// otbMsgDevMacro(<< "CartoPoint: (" << geoPoint[0] << "," << geoPoint[1] << ")");
height=m_DEMHandler->GetHeightAboveMSL(geoPoint); // Altitude calculation
otbMsgDevMacro(<< "height" << height);
// otbMsgDevMacro(<< "height: " << height);
// MNT sets a default value (-32768) at point where it doesn't have altitude information.
// OSSIM has chosen to change this default value in OSSIM_DBL_NAN (-4.5036e15).
if (!ossim::isnan(height))
......
......@@ -141,8 +141,23 @@ ADD_TEST(coTvVectorDataExtractROI ${COMMON_TESTS2}
${BASELINE_FILES}/coVectorDataExtractROIOutput.shp
${TEMP}/coVectorDataExtractROIOutput.shp
otbVectorDataExtractROI
${INPUTDATA}/ToulousePoints-examples.shp #${LARGEINPUT}/TOULOUSE/QuickBird/GIS_FILES/000000128955_01_ORDER_SHAPE.shp
${INPUTDATA}/ToulousePoints-examples.shp
${TEMP}/coVectorDataExtractROIOutput.shp
374369.48850211215904 4853951.786124914 # Origin of the RemoteSensingRegion
1000.25 -25000.2 # Size of the Cartoregion
)
# This test is almost the same as the previous one and should return the
# same output.
# This is to check that region can also be passed with origin
# at the south-west corner
ADD_TEST(coTvVectorDataExtractROI2 ${COMMON_TESTS2}
--compare-ogr ${NOTOL}
${BASELINE_FILES}/coVectorDataExtractROIOutput2.shp
${TEMP}/coVectorDataExtractROIOutput2.shp
otbVectorDataExtractROI
${INPUTDATA}/ToulousePoints-examples.shp
${TEMP}/coVectorDataExtractROIOutput2.shp
374369.48850211215904 4828951.58612491376698 # Origin of the RemoteSensingRegion
1000.25 25000.2 # Size of the Cartoregion
)
......
......@@ -27,25 +27,25 @@ int otbVectorDataExtractROI( int argc, char * argv[] )
const char * infname = argv[1];
const char * outfname = argv[2];
const unsigned int startX = atoi(argv[3]);
const unsigned int startY = atoi(argv[4]);
const unsigned int sizeX = atoi(argv[5]);
const unsigned int sizeY = atoi(argv[6]);
const long int startX = atoi(argv[3]);
const long int startY = atoi(argv[4]);
const long int sizeX = atoi(argv[5]);
const long int sizeY = atoi(argv[6]);
typedef double Type;
typedef otb::VectorData<> VectorDataType;
typedef otb::VectorDataExtractROI< VectorDataType > FilterType;
typedef otb::VectorDataFileReader<VectorDataType> VectorDataFileReaderType;
typedef otb::VectorDataFileWriter<VectorDataType> VectorDataWriterType;
typedef otb::RemoteSensingRegion<Type> TypedRegion;
/** Instanciation of pointer objects*/
FilterType::Pointer filter = FilterType::New();
VectorDataFileReaderType::Pointer reader = VectorDataFileReaderType::New();
VectorDataWriterType::Pointer writer = VectorDataWriterType::New();
/** Edit The cartoRegion*/
TypedRegion region;
TypedRegion::SizeType size;
......@@ -55,26 +55,26 @@ int otbVectorDataExtractROI( int argc, char * argv[] )
size[1] = sizeY;
index[0] = startX;
index[1] = startY;
/*Set the projection */
std::string projectionRefWkt ="PROJCS[\"UTM Zone 31, Northern Hemisphere\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9108\"]],AXIS[\"Lat\",NORTH],AXIS[\"Long\",EAST],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",3],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]]";
region.SetSize(size);
region.SetOrigin(index);
region.SetRegionProjection(projectionRefWkt);
/** */
reader->SetFileName(infname);
filter->SetInput(reader->GetOutput());
filter->SetRegion(region);
/** Write the shapefile*/
writer->SetFileName(outfname);
writer->SetInput(filter->GetOutput());
writer->Update();
return EXIT_SUCCESS;
}
......
......@@ -73,16 +73,20 @@ int otbVectorDataToImageFilterSensorModel(int argc, char * argv[])
size[1] = 500;
ImageType::PointType origin;
origin[0] = imageReader->GetOutput()->GetOrigin()[0];
origin[1] = imageReader->GetOutput()->GetOrigin()[1];
// origin[0] = imageReader->GetOutput()->GetOrigin()[0];
// origin[1] = imageReader->GetOutput()->GetOrigin()[1];
origin[0] = 11052;
origin[1] = 7585;
ImageType::SpacingType spacing;
// spacing[0] = imageReader->GetOutput()->GetSpacing()[0];
// spacing[1] = imageReader->GetOutput()->GetSpacing()[1];
spacing[0] = imageReader->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/static_cast<double>(size[0]);
// /imageReader->GetOutput()->GetSpacing()[0];
spacing[1] = imageReader->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/static_cast<double>(size[1]);
// /imageReader->GetOutput()->GetSpacing()[1];
// spacing[0] = imageReader->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/static_cast<double>(size[0]);
// spacing[1] = imageReader->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/static_cast<double>(size[1]);
spacing[0] = 500/static_cast<double>(size[0]);
spacing[1] = 500/static_cast<double>(size[1]);
std::cout << "Size: " << size <<std::endl;
std::cout << "Origin: " << origin <<std::endl;
std::cout << "Spacing: " << spacing <<std::endl;
typedef otb::VectorDataToImageFilter<VectorDataType, ImageType> VectorDataToImageFilterType;
VectorDataToImageFilterType::Pointer vectorDataRendering = VectorDataToImageFilterType::New();
......
......@@ -1403,13 +1403,13 @@ ADD_TEST(ioTvDEMToOrthoImageGeneratorTest ${IO_TESTS12}
otbDEMToOrthoImageGeneratorTest
${INPUTDATA}/DEM/srtm_directory
${TEMP}/ioTvDEMToOrthoImageGenerator.tif
773441 #6.5
5044457 #44.5
778245.3501 #6.5
4933369.633 #44.5
500
500
200 #0.002
-200 #0.002
48 #31 # utm zone
30 #resolution in m
-30 #resolution in m
31 # utm zone
N # hemisphere
)
......@@ -1624,10 +1624,10 @@ ADD_TEST(ioTvVectorDataFileWriterTwice ${IO_TESTS15}
ADD_TEST(ioTvVectorDataFileWriterPolygons ${IO_TESTS15}
--compare-ogr ${TOL}
${BASELINE_FILES}/ioTvVectorDataFileWriterOutput.shp
${TEMP}/ioTvVectorDataFileWriterOutput.shp
${BASELINE_FILES}/ioTvVectorDataFileWriterPolygonsOutput.shp
${TEMP}/ioTvVectorDataFileWriterPolygonsOutput.shp
otbVectorDataFileWriterPolygons
${TEMP}/ioTvVectorDataFileWriterOutput.shp
${TEMP}/ioTvVectorDataFileWriterPolygonsOutput.shp
)
......
......@@ -35,8 +35,8 @@ int otbDEMToOrthoImageGeneratorTest(int argc, char * argv[])
const unsigned int Dimension = 2;
typedef otb::Image<double, Dimension> ImageType;
typedef otb::UtmForwardProjection MapProjectionType;
typedef otb::DEMToOrthoImageGenerator<ImageType, otb::UtmForwardProjection> DEMToImageGeneratorType;
typedef otb::UtmInverseProjection MapProjectionType;
typedef otb::DEMToOrthoImageGenerator<ImageType, MapProjectionType> DEMToImageGeneratorType;
typedef DEMToImageGeneratorType::DEMHandlerType DEMHandlerType;
typedef DEMHandlerType::PointType PointType;
typedef DEMToImageGeneratorType::SizeType SizeType;
......
Markdown is supported
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