From aff4925a9b9c413e95829f407ec2a42064b37614 Mon Sep 17 00:00:00 2001 From: Mickael Savinaud <mickael.savinaud@c-s.fr> Date: Tue, 24 May 2011 18:32:25 +0200 Subject: [PATCH] WIP: add test about vectorDataProjectionFilter GeoToMap but compare fail even if qgis indicate the same shapefile --- Testing/Code/Projections/CMakeLists.txt | 91 +++++++++++--- .../Code/Projections/otbProjectionsTests3.cxx | 1 + ...ectorDataIntoImageProjectionFilterTest.cxx | 114 ++++++++++++------ ...VectorDataProjectionFilterFromGeoToMap.cxx | 66 ++++++++++ 4 files changed, 220 insertions(+), 52 deletions(-) create mode 100644 Testing/Code/Projections/otbVectorDataProjectionFilterFromGeoToMap.cxx diff --git a/Testing/Code/Projections/CMakeLists.txt b/Testing/Code/Projections/CMakeLists.txt index a28d559382..16e069660f 100644 --- a/Testing/Code/Projections/CMakeLists.txt +++ b/Testing/Code/Projections/CMakeLists.txt @@ -932,48 +932,104 @@ ADD_TEST(prTvVectorDataProjectionFilterFromMapToImage ${PROJECTIONS_TESTS3} ${TEMP}/prTvVectorDataProjectionFilterFromMapToImage.kml ) -# extract from image QB and MidiPyrenees vectordata -#ADD_TEST(prTvVectorData_ProjRef_ProjIntoImage_WithProjRef_Extract_Filter0 ${PROJECTIONS_TESTS3} +# With QGIS the two vectordata are equal +ADD_TEST(prTvVectorDataProjectionFilterFromGeoToMap ${PROJECTIONS_TESTS3} + #--compare-ascii ${NOTOL} + #${INPUTDATA}/ToulousePoints-examples.shp + #${TEMP}/prTvVectorDataProjectionFilterFromGeoToMap.shp + otbVectorDataProjectionFilterFromGeoToMap + ${TEMP}/prTvVectorDataProjectionFilterFromMapToGeo.kml + ${INPUTDATA}/QB_Toulouse_Ortho_PAN.tif + ${TEMP}/prTvVectorDataProjectionFilterFromGeoToMap.shp +) +SET_TESTS_PROPERTIES(prTvVectorDataProjectionFilterFromGeoToMap PROPERTIES DEPENDS prTvVectorDataProjectionFilterFromMapToGeo) + + + +################### +# TESTS ABOUT REPROJECTION INTO AN IMAGE OF A VECTORDATA +# QB ortho and extract and extract of MidiPyrenees roads file +#ADD_TEST(prTvVectorData_ProjRef_ProjIntoImage_WithProjRef ${PROJECTIONS_TESTS3} # otbVectorDataIntoImageProjectionFilterTest -# ${INPUTDATA}/QB_extract_512x512.tif +# ${INPUTDATA}/Dempster-Shafer/ROI_QB_TOULOUSE.TIF # ${LARGEINPUT}/VECTOR/MidiPyrenees/roads.shp # ${INPUTDATA}/DEM/srtm_directory -# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithProjRef10.shp -# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithProjRef10.txt +# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithProjRef.shp +# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithProjRef.txt # 0 #) -#ADD_TEST(prTvVectorData_ProjRef_ProjIntoImage_WithProjRef_ExtractLeftUp_Filter0 ${PROJECTIONS_TESTS3} +#ADD_TEST(prTvVectorData_ProjRef_ProjIntoImage_WithProjRef1 ${PROJECTIONS_TESTS3} +# otbVectorDataIntoImageProjectionFilterTest +# ${INPUTDATA}/Dempster-Shafer/ROI_QB_TOULOUSE.TIF +# ${LARGEINPUT}/VECTOR/MidiPyrenees/roads.shp +# ${INPUTDATA}/DEM/srtm_directory +# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithProjRef1.shp +# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithProjRef1.txt +# 1 +#) + +#ADD_TEST(prTvVectorData_ProjRef_ProjIntoImage_WithSensorModel ${PROJECTIONS_TESTS3} # otbVectorDataIntoImageProjectionFilterTest -# ${INPUTDATA}/QB_extract_leftup_256x2048.tif +# ${INPUTDATA}/Dempster-Shafer/ROI_QB_TOULOUSE_EXTRACT.HDR # ${LARGEINPUT}/VECTOR/MidiPyrenees/roads.shp # ${INPUTDATA}/DEM/srtm_directory -# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithProjRef20.shp -# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithProjRef20.txt +# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithWithSensorModel.shp +# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithWithSensorModel.txt # 0 #) -#ADD_TEST(prTvVectorData_ProjRef_ProjIntoImage_WithProjRef_Extract_Filter1 ${PROJECTIONS_TESTS3} +#ADD_TEST(prTvVectorData_ProjRef_ProjIntoImage_WithSensorModel1 ${PROJECTIONS_TESTS3} # otbVectorDataIntoImageProjectionFilterTest -# ${INPUTDATA}/QB_extract_512x512.tif +# ${INPUTDATA}/Dempster-Shafer/ROI_QB_TOULOUSE_EXTRACT.HDR # ${LARGEINPUT}/VECTOR/MidiPyrenees/roads.shp # ${INPUTDATA}/DEM/srtm_directory -# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithProjRef11.shp -# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithProjRef11.txt +# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithWithSensorModel1.shp +# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithWithSensorModel1.txt # 1 #) -#ADD_TEST(prTvVectorData_ProjRef_ProjIntoImage_WithProjRef_ExtractLeftUp_Filter1 ${PROJECTIONS_TESTS3} +#ADD_TEST(prTvVectorData_ProjRef_ProjIntoImage_Index ${PROJECTIONS_TESTS3} +# otbVectorDataIntoImageProjectionFilterTest +# ${INPUTDATA}/Dempster-Shafer/ROI_QB_TOULOUSE_EXTRACT_INDEX.HDR +# ${LARGEINPUT}/VECTOR/MidiPyrenees/roads.shp +# ${INPUTDATA}/DEM/srtm_directory +# ${TEMP}/prTvVectorDataProjRefProjIntoImageIndex.shp +# ${TEMP}/prTvVectorDataProjRefProjIntoImageindex.txt +# 0 +#) + +#ADD_TEST(prTvVectorData_ProjRef_ProjIntoImage_Index1 ${PROJECTIONS_TESTS3} # otbVectorDataIntoImageProjectionFilterTest -# ${INPUTDATA}/QB_extract_leftup_256x2048.tif +# ${INPUTDATA}/Dempster-Shafer/ROI_QB_TOULOUSE_EXTRACT_INDEX.HDR # ${LARGEINPUT}/VECTOR/MidiPyrenees/roads.shp # ${INPUTDATA}/DEM/srtm_directory -# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithProjRef21.shp -# ${TEMP}/prTvVectorDataProjRefProjIntoImageWithProjRef21.txt +# ${TEMP}/prTvVectorDataProjRefProjIntoImageIndex1.shp +# ${TEMP}/prTvVectorDataProjRefProjIntoImageindex1.txt +# 1 +#) + +#ADD_TEST(prTvVectorData_NoProjRef_ProjIntoImage_Index ${PROJECTIONS_TESTS3} +# otbVectorDataIntoImageProjectionFilterTest +# ${INPUTDATA}/QB_Toulouse_extract_leftup_400x128_index.hdr +# ${INPUTDATA}/QB_Toulouse_extract_leftup_400x128_index.shp +# ${INPUTDATA}/DEM/srtm_directory +# ${TEMP}/prTvVectorDataNoProjRefProjIntoImageIndex.shp +# ${TEMP}/prTvVectorDataNoProjRefProjIntoImageindex.txt +# 0 +#) +#ADD_TEST(prTvVectorData_NoProjRef_ProjIntoImage_Index1 ${PROJECTIONS_TESTS3} +# otbVectorDataIntoImageProjectionFilterTest +# ${INPUTDATA}/QB_Toulouse_extract_leftup_400x128_index.hdr +# ${INPUTDATA}/QB_Toulouse_extract_leftup_400x128_index.shp +# ${INPUTDATA}/DEM/srtm_directory +# ${TEMP}/prTvVectorDataNoProjRefProjIntoImageIndex1.shp +# ${TEMP}/prTvVectorDataNoProjRefProjIntoImageindex1.txt # 1 #) + ADD_TEST(prTuGeocentricTransformNew ${PROJECTIONS_TESTS3} otbGeocentricTransformNew ) @@ -1181,6 +1237,7 @@ otbVectorDataProjectionFilter.cxx otbVectorDataProjectionFilterFromMapToSensor.cxx otbVectorDataProjectionFilterFromMapToGeo.cxx otbVectorDataProjectionFilterFromMapToImage.cxx +otbVectorDataProjectionFilterFromGeoToMap.cxx otbVectorDataIntoImageProjectionFilterTest.cxx otbGeocentricTransformNew.cxx otbGeocentricTransform.cxx diff --git a/Testing/Code/Projections/otbProjectionsTests3.cxx b/Testing/Code/Projections/otbProjectionsTests3.cxx index 2b7d83655b..c10325ea16 100644 --- a/Testing/Code/Projections/otbProjectionsTests3.cxx +++ b/Testing/Code/Projections/otbProjectionsTests3.cxx @@ -38,6 +38,7 @@ void RegisterTests() REGISTER_TEST(otbVectorDataProjectionFilterFromMapToSensor); REGISTER_TEST(otbVectorDataProjectionFilterFromMapToGeo); REGISTER_TEST(otbVectorDataProjectionFilterFromMapToImage); + REGISTER_TEST(otbVectorDataProjectionFilterFromGeoToMap); REGISTER_TEST(otbGeocentricTransformNew); REGISTER_TEST(otbGeocentricTransform); REGISTER_TEST(otbTileMapTransform); diff --git a/Testing/Code/Projections/otbVectorDataIntoImageProjectionFilterTest.cxx b/Testing/Code/Projections/otbVectorDataIntoImageProjectionFilterTest.cxx index 90ece04590..d1c8ffad4d 100644 --- a/Testing/Code/Projections/otbVectorDataIntoImageProjectionFilterTest.cxx +++ b/Testing/Code/Projections/otbVectorDataIntoImageProjectionFilterTest.cxx @@ -29,6 +29,9 @@ #include "otbVectorDataIntoImageProjectionFilter.h" +#include "otbVectorDataExtractROI.h" +#include "otbVectorDataProjectionFilter.h" + int otbVectorDataIntoImageProjectionFilterTest(int argc, char * argv[]) { typedef float PixelType; @@ -45,19 +48,24 @@ int otbVectorDataIntoImageProjectionFilterTest(int argc, char * argv[]) typedef otb::VectorDataIntoImageProjectionFilter <VectorDataType, VectorImageType> VectorDataReProjFilter; + typedef otb::VectorDataProjectionFilter + <VectorDataType, VectorDataType> VectorDataProjectionFilterType; + typedef otb::VectorDataExtractROI<VectorDataType> VectorDataExtractROIType; + typedef VectorDataExtractROIType::RegionType RemoteSensingRegionType; + typedef itk::PreOrderTreeIterator<VectorDataType::DataTreeType> TreeIteratorType; std::string imageInputFilename = argv[1]; std::string vectorDataInputFilename = argv[2]; std::string demDirectory = argv[3]; std::string vectorDataOutputFilename = argv[4]; - std::string txtOutputFilename = argv[5]; + std::string vectorDataOutputFilename2 = argv[5]; std::cout << imageInputFilename << "\n" << vectorDataInputFilename << "\n" << demDirectory << "\n" << vectorDataOutputFilename << "\n" - << txtOutputFilename << std::endl; + << vectorDataOutputFilename2 << std::endl; // Read the image ReaderType::Pointer reader = ReaderType::New(); @@ -69,29 +77,8 @@ int otbVectorDataIntoImageProjectionFilterTest(int argc, char * argv[]) vdReader->SetFileName(vectorDataInputFilename); vdReader->Update(); - std::cout<<"Set input data to read data DONE"<<std::endl; - - std::ofstream file; - file.open(txtOutputFilename.c_str()); - - file << "--- IMAGE INPUT ---" << std::endl; - file << "Spacing of the input image: "<< reader->GetOutput()->GetSpacing() << std::endl; - file << "Origin of the input image: "<< reader->GetOutput()->GetOrigin() << std::endl; - file << "Size of the input image: "<< reader->GetOutput()->GetLargestPossibleRegion() << std::endl; - file << "ProjRef of the input image: "<< reader->GetOutput()->GetProjectionRef() << std::endl; - - VectorDataReProjFilter::Pointer vdReProjFilter = VectorDataReProjFilter::New(); - //---------- - // Display input - file << "\n--- VECTOR DATA INPUT ---" << std::endl; - file << "ProjRef of the input vector data: "<< vdReader->GetOutput()->GetProjectionRef() << std::endl; - file << "Origin of the input vector data: "<< vdReader->GetOutput()->GetOrigin() << std::endl; - file << "Spacing of the input vector data: "<< vdReader->GetOutput()->GetSpacing() << std::endl; - - //---------- - // Set input of the filter vdReProjFilter->SetInputImage(reader->GetOutput()); vdReProjFilter->SetInputVectorData(vdReader->GetOutput()); @@ -110,8 +97,6 @@ int otbVectorDataIntoImageProjectionFilterTest(int argc, char * argv[]) vdReProjFilter->SetUseOutputSpacingAndOriginFromImage(false); } - std::cout<<"Set input data for the filter DONE (" << stateOutput << ")" <<std::endl; - //---------- // WRITE //vdReProjFilter->Update(); @@ -121,18 +106,77 @@ int otbVectorDataIntoImageProjectionFilterTest(int argc, char * argv[]) vdwriter->SetInput(vdReProjFilter->GetOutput()); vdwriter->Update(); - std::cout<<"Update"<<std::endl; + //------- + // do the same with old code + + VectorDataProjectionFilterType::Pointer vproj; + VectorDataExtractROIType::Pointer vdextract; + + // Extract The part of the VectorData that actually overlaps with + // the image extent + vdextract = VectorDataExtractROIType::New(); + vdextract->SetInput(vdReader->GetOutput()); + + // Find the geographic region of interest + + // Ge the index of the corner of the image + ImageType::IndexType ul, ur, ll, lr; + ImageType::PointType pul, pur, pll, plr; + ul = reader->GetOutput()->GetLargestPossibleRegion().GetIndex(); + ur = ul; + ll = ul; + lr = ul; + ur[0] += reader->GetOutput()->GetLargestPossibleRegion().GetSize()[0]; + lr[0] += reader->GetOutput()->GetLargestPossibleRegion().GetSize()[0]; + lr[1] += reader->GetOutput()->GetLargestPossibleRegion().GetSize()[1]; + ll[1] += reader->GetOutput()->GetLargestPossibleRegion().GetSize()[1]; + + // Transform to physical point + reader->GetOutput()->TransformIndexToPhysicalPoint(ul, pul); + reader->GetOutput()->TransformIndexToPhysicalPoint(ur, pur); + reader->GetOutput()->TransformIndexToPhysicalPoint(ll, pll); + reader->GetOutput()->TransformIndexToPhysicalPoint(lr, plr); + + // Build the cartographic region + RemoteSensingRegionType rsRegion; + RemoteSensingRegionType::IndexType rsOrigin; + RemoteSensingRegionType::SizeType rsSize; + rsOrigin[0] = std::min(pul[0], plr[0]); + rsOrigin[1] = std::min(pul[1], plr[1]); + rsSize[0] = vcl_abs(pul[0] - plr[0]); + rsSize[1] = vcl_abs(pul[1] - plr[1]); + + rsRegion.SetOrigin(rsOrigin); + rsRegion.SetSize(rsSize); + rsRegion.SetRegionProjection(reader->GetOutput()->GetProjectionRef()); + rsRegion.SetKeywordList(reader->GetOutput()->GetImageKeywordlist()); + + // Set the cartographic region to the extract roi filter + vdextract->SetRegion(rsRegion); + if (!demDirectory.empty()) + { + vdextract->SetDEMDirectory(demDirectory); + } + + // Reproject VectorData in image projection + vproj = VectorDataProjectionFilterType::New(); + vproj->SetInput(vdextract->GetOutput()); + vproj->SetInputProjectionRef(vdReader->GetOutput()->GetProjectionRef()); + vproj->SetOutputKeywordList(reader->GetOutput()->GetImageKeywordlist()); + vproj->SetOutputProjectionRef(reader->GetOutput()->GetProjectionRef()); + vproj->SetOutputOrigin(reader->GetOutput()->GetOrigin()); + vproj->SetOutputSpacing(reader->GetOutput()->GetSpacing()); + if (!demDirectory.empty()) + { + vproj->SetDEMDirectory(demDirectory); + } //---------- - // check output - file << "\n--- VECTRODATA OUTPUT ---" << std::endl; - file << "ProjRef of the output vector data: "<< vdReProjFilter->GetOutput()->GetProjectionRef() << std::endl; - file << "Origin of the output vector data: "<< vdReProjFilter->GetOutput()->GetOrigin() << std::endl; - file << "Spacing of the output vector data: "<< vdReProjFilter->GetOutput()->GetSpacing() << std::endl; - - file.close(); - - std::cout<<"END"<<std::endl; + // WRITE + VectorDataWriterType::Pointer vdwriter2 = VectorDataWriterType::New(); + vdwriter2->SetFileName(vectorDataOutputFilename2); + vdwriter2->SetInput(vproj->GetOutput()); + vdwriter2->Update(); return EXIT_SUCCESS; diff --git a/Testing/Code/Projections/otbVectorDataProjectionFilterFromGeoToMap.cxx b/Testing/Code/Projections/otbVectorDataProjectionFilterFromGeoToMap.cxx new file mode 100644 index 0000000000..5f4535fde2 --- /dev/null +++ b/Testing/Code/Projections/otbVectorDataProjectionFilterFromGeoToMap.cxx @@ -0,0 +1,66 @@ +/*========================================================================= + + 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 "itkExceptionObject.h" +#include "otbVectorDataProjectionFilter.h" +#include "otbVectorData.h" +#include "otbVectorDataFileReader.h" +#include "otbVectorDataFileWriter.h" + +#include "otbImage.h" +#include "otbImageFileReader.h" + +int otbVectorDataProjectionFilterFromGeoToMap(int argc, char * argv[]) +{ + + if (argc < 2) + { + std::cout << argv[0] << " <input vector filename> <output vector filename> " << std::endl; + + return EXIT_FAILURE; + } + + typedef otb::VectorData<double> InputVectorDataType; + typedef otb::VectorData<double> OutputVectorDataType; + + typedef otb::VectorDataFileReader<InputVectorDataType> VectorDataFileReaderType; + VectorDataFileReaderType::Pointer reader = VectorDataFileReaderType::New(); + + reader->SetFileName(argv[1]); + reader->UpdateOutputInformation(); + + typedef otb::Image<unsigned short int, 2> ImageType; + typedef otb::ImageFileReader<ImageType> ImageReaderType; + ImageReaderType::Pointer imageReader = ImageReaderType::New(); + imageReader->SetFileName(argv[2]); + imageReader->UpdateOutputInformation(); + + typedef otb::VectorDataProjectionFilter<InputVectorDataType, OutputVectorDataType> VectorDataFilterType; + VectorDataFilterType::Pointer vectorDataProjection = VectorDataFilterType::New(); + + vectorDataProjection->SetInput(reader->GetOutput()); + vectorDataProjection->SetOutputProjectionRef(imageReader->GetOutput()->GetProjectionRef()); + + typedef otb::VectorDataFileWriter<OutputVectorDataType> VectorDataFileWriterType; + VectorDataFileWriterType::Pointer writer = VectorDataFileWriterType::New(); + + writer->SetFileName(argv[3]); + writer->SetInput(vectorDataProjection->GetOutput()); + writer->Update(); + + return EXIT_SUCCESS; +} -- GitLab