diff --git a/Code/BasicFilters/otbStreamingResampleImageFilter.h b/Code/BasicFilters/otbStreamingResampleImageFilter.h index 40954530629e3cf42f34f1790745214e4e4f3e33..1d0127919a61da9255790b28663b695d37e3a0e4 100644 --- a/Code/BasicFilters/otbStreamingResampleImageFilter.h +++ b/Code/BasicFilters/otbStreamingResampleImageFilter.h @@ -65,12 +65,11 @@ public: typedef typename Superclass::InterpolatorType InterpolatorType; typedef typename InterpolatorType::PointType PointType; - /** Set size of neighborhood needed to interpolate points */ itkSetMacro(InterpolatorNeighborhoodRadius,unsigned int); - - /** Get size of neighborhood needed to interpolate points */ itkGetMacro(InterpolatorNeighborhoodRadius,unsigned int); + itkSetMacro(AddedRadius,unsigned int); + itkGetMacro(AddedRadius,unsigned int); /** ResampleImageFilter needs a different input requested region than * the output requested region. As such, ResampleImageFilter needs @@ -93,8 +92,8 @@ private: // Determine size of pad needed for interpolators neighborhood unsigned int m_InterpolatorNeighborhoodRadius; - // Determine if interpolator radius is determined by class user -// bool m_RadiusIsDeterminedByUser; + // Used to be sure that each final region will be contiguous + unsigned int m_AddedRadius; }; diff --git a/Code/BasicFilters/otbStreamingResampleImageFilter.txx b/Code/BasicFilters/otbStreamingResampleImageFilter.txx index 637353f4ee7eadc9b06b3634ec89c416c9316be6..6e27f3d5d149d131fa862f0713f4704851eab572 100644 --- a/Code/BasicFilters/otbStreamingResampleImageFilter.txx +++ b/Code/BasicFilters/otbStreamingResampleImageFilter.txx @@ -40,7 +40,7 @@ StreamingResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType // Default neighborhood interpolation radius is one pixel m_InterpolatorNeighborhoodRadius = 1 ; -// m_RadiusIsDeterminedByUser = false; + m_AddedRadius = 2; } @@ -53,7 +53,7 @@ StreamingResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType if ( this->GetInput() ) { - otbDebugMacro(<< "-------------- GenerateInputRequestedRegion ---------------" << std::endl); + otbMsgDebugMacro(<< "-------------- GenerateInputRequestedRegion ---------------" << std::endl); InputImagePointer inputImage = const_cast< typename Superclass::InputImageType *>( this->GetInput() ); OutputImagePointer outputImage = const_cast< typename Superclass::OutputImageType *>( this->GetOutput() ); @@ -66,27 +66,27 @@ StreamingResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType std::vector<IndexType> vPoints; typename std::vector<IndexType>::iterator it; - otbDebugMacro(<< "Size : " << size[0] << " " << size[1]); + otbMsgDebugMacro(<< "Size : " << size[0] << " " << size[1]); indexTmp[0]=index[0]; indexTmp[1]=index[1]; vPoints.push_back(indexTmp); - otbDebugMacro(<< "indexUL : (" << indexTmp[0] << "," << indexTmp[1] << ")"); + //otbGenericMsgDebugMacro(<< "indexUL : (" << indexTmp[0] << "," << indexTmp[1] << ")"); indexTmp[0]=index[0]+size[0]; indexTmp[1]=index[1]; vPoints.push_back(indexTmp); - otbDebugMacro(<< "indexUR : (" << indexTmp[0] << "," << indexTmp[1] << ")"); + //otbGenericMsgDebugMacro(<< "indexUR : (" << indexTmp[0] << "," << indexTmp[1] << ")"); indexTmp[0]=index[0]+size[0]; indexTmp[1]=index[1]+size[1]; vPoints.push_back(indexTmp); - otbDebugMacro(<< "indexLR : (" << indexTmp[0] << "," << indexTmp[1] << ")"); + //otbGenericMsgDebugMacro(<< "indexLR : (" << indexTmp[0] << "," << indexTmp[1] << ")"); indexTmp[0]=index[0]; indexTmp[1]=index[1]+size[1]; vPoints.push_back(indexTmp); - otbDebugMacro(<< "indexLL : (" << indexTmp[0] << "," << indexTmp[1] << ")"); + //otbGenericMsgDebugMacro(<< "indexLL : (" << indexTmp[0] << "," << indexTmp[1] << ")"); typedef itk::ContinuousIndex<TInterpolatorPrecisionType, 2> ContinuousIndexType; typename ContinuousIndexType::ValueType minX = itk::NumericTraits<typename ContinuousIndexType::ValueType>::max(); @@ -106,24 +106,32 @@ StreamingResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType // Calculate transformed points needed for previous filter in the pipeline outputImage->TransformIndexToPhysicalPoint( *it, outputPoint ); + otbMsgDebugMacro(<< "Pour l'Index Ncurrent:(" << (*it)[0]<<","<< (*it)[1] << ")"<< std::endl + << "Le point physique correspondant est: ("<< outputPoint[0]<< ","<< outputPoint[1]<< ")"); + // Compute corresponding input pixel continuous index inputPoint = this->GetTransform()->TransformPoint(outputPoint); - inputImage->TransformPhysicalPointToContinuousIndex(inputPoint, indexTmpTr); - + inputImage->TransformPhysicalPointToContinuousIndex(inputPoint, indexTmpTr); + + otbMsgDebugMacro(<< "L'index correspondant a ce point est:" << std::endl + << indexTmpTr[0] << ","<< indexTmpTr[1] ); + if (indexTmpTr[0]>maxX) maxX = indexTmpTr[0]; - else if (indexTmpTr[0]<minX) + + if (indexTmpTr[0]<minX) minX = indexTmpTr[0]; if (indexTmpTr[1]>maxY) maxY = indexTmpTr[1]; - else if (indexTmpTr[1]<minY) + + if (indexTmpTr[1]<minY) minY = indexTmpTr[1]; - otbDebugMacro(<< "indexTr : (" << indexTmpTr[0] << "," << indexTmpTr[1] << ")"); + //otbGenericMsgDebugMacro(<< "indexTr : (" << indexTmpTr[0] << "," << indexTmpTr[1] << ")"); } - otbDebugMacro(<< "MinX : " << minX << " MinY : " << minY << " MaxX : " << maxX << " MaxY " << maxY); + otbMsgDebugMacro(<< "MinX : " << minX << " MinY : " << minY << " MaxX : " << maxX << " MaxY " << maxY); // Create region needed in previous filter in the pipeline, which is the bounding box of previous transformed points InputImageRegionType region; @@ -132,7 +140,7 @@ StreamingResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType size[0] = static_cast<long unsigned int>(maxX - minX); size[1] = static_cast<long unsigned int>(maxY - minY); - otbDebugMacro(<< "Index : (" << index[0] << "," << index[1] << ") Size : (" << size[0] << "," << size[1] << ")"); + otbMsgDebugMacro(<< "Index : (" << index[0] << "," << index[1] << ") Size : (" << size[0] << "," << size[1] << ")"); region.SetSize(size); region.SetIndex(index); @@ -146,13 +154,13 @@ StreamingResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType neededRadius = m_InterpolatorNeighborhoodRadius; } - otbDebugMacro(<< "Interpolation needed radius : " << neededRadius); - region.PadByRadius(neededRadius); + otbMsgDebugMacro(<< "Interpolation needed radius : " << neededRadius); + region.PadByRadius(neededRadius+m_AddedRadius); - otbDebugMacro(<< "Initial Region : Index(" << inputImage->GetLargestPossibleRegion().GetIndex()[0] << "," << inputImage->GetLargestPossibleRegion().GetIndex()[1] << ") Size(" << inputImage->GetLargestPossibleRegion().GetSize()[0] << "," << inputImage->GetLargestPossibleRegion().GetSize()[1] << ")"); + otbMsgDebugMacro(<< "Initial Region : Index(" << inputImage->GetLargestPossibleRegion().GetIndex()[0] << "," << inputImage->GetLargestPossibleRegion().GetIndex()[1] << ") Size(" << inputImage->GetLargestPossibleRegion().GetSize()[0] << "," << inputImage->GetLargestPossibleRegion().GetSize()[1] << ")"); // To be sure that requested region in pipeline is not largest than real input image - otbDebugMacro(<< "Final Region (Before Crop) : Index(" << region.GetIndex()[0] << "," << region.GetIndex()[1] << ") Size(" << region.GetSize()[0] << "," << region.GetSize()[1] << ")"); + otbMsgDebugMacro(<< "Final Region (Before Crop) : Index(" << region.GetIndex()[0] << "," << region.GetIndex()[1] << ") Size(" << region.GetSize()[0] << "," << region.GetSize()[1] << ")"); // If requested region is not contained in input image, then result region is null if (!region.Crop(inputImage->GetLargestPossibleRegion())) @@ -167,7 +175,7 @@ StreamingResampleImageFilter<TInputImage,TOutputImage,TInterpolatorPrecisionType inputImage->SetRequestedRegion(region); - otbDebugMacro(<< "Final Region (After Crop) : Index(" << region.GetIndex()[0] << "," << region.GetIndex()[1] << ") Size(" << region.GetSize()[0] << "," << region.GetSize()[1] << ")"); + otbMsgDebugMacro(<< "Final Region (After Crop) : Index(" << region.GetIndex()[0] << "," << region.GetIndex()[1] << ") Size(" << region.GetSize()[0] << "," << region.GetSize()[1] << ")"); } } diff --git a/Code/IO/otbDEMReader.cxx b/Code/IO/otbDEMReader.cxx new file mode 100644 index 0000000000000000000000000000000000000000..6500a1fe4226d8848bbd46e059497d618081cea6 --- /dev/null +++ b/Code/IO/otbDEMReader.cxx @@ -0,0 +1,86 @@ +/*========================================================================= + + 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 __otbDEMReader_cxx +#define __otbDEMReader_cxx + +#include "otbDEMReader.h" +#include "otbMacro.h" + +namespace otb +{ + + + +DEMReader +::DEMReader() +{ + m_ElevManager=ossimElevManager::instance(); +} + + + +DEMReader +::~DEMReader() +{ + // not needed, m_ElevManager created with instance() method + // delete m_ElevManager; + +} + + + +bool +DEMReader +::OpenDEMDirectory(const char* DEMDirectory) +{ + ossimFilename ossimDEMDir; + ossimDEMDir=ossimFilename(DEMDirectory); + bool result= false; + if (m_ElevManager->openDirectory(ossimDEMDir)) + { + result= true; + } + return result; +} + + +double +DEMReader +::GetHeightAboveMSL(const PointType& geoPoint) +{ + float height; + ossimGpt ossimWorldPoint; + ossimWorldPoint.lat=geoPoint[0]; + ossimWorldPoint.lon=geoPoint[1]; + height=m_ElevManager->getHeightAboveMSL(ossimWorldPoint); + return height; +} + + +void +DEMReader +::PrintSelf(std::ostream& os, Indent indent) const +{ + Superclass::PrintSelf(os,indent); + + os << indent << "DEMReader" << std::endl; +} + +} // namespace otb + +#endif diff --git a/Code/IO/otbDEMReader.h b/Code/IO/otbDEMReader.h index e0f2fa3ffadf448e678ab23caaf97d9fbfd21c06..2b6326a97d4b4a8d44988ad53613bab81d237e36 100755 --- a/Code/IO/otbDEMReader.h +++ b/Code/IO/otbDEMReader.h @@ -39,77 +39,43 @@ namespace otb * \ingroup Images * */ -template <class TDEMImage> -class ITK_EXPORT DEMReader: -public itk::ImageSource<Image<typename TDEMImage::PixelType,2, 0> > + +class ITK_EXPORT DEMReader: public itk::Object { public : /** Standard class typedefs. */ typedef itk::Indent Indent; - typedef TDEMImage DEMImageType; - typedef typename DEMImageType::Pointer DEMImagePointerType; - typedef typename DEMImageType::PixelType PixelType; - - typedef DEMReader Self; - typedef itk::ImageSource<Image<typename DEMImageType::PixelType,2, 0> > Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; - typedef Image<PixelType,2> OutputImageType; - - typedef typename Superclass::Pointer OutputImagePointer; - typedef typename OutputImageType::SpacingType SpacingType; - typedef typename OutputImageType::SizeType SizeType; - typedef typename OutputImageType::PointType PointType; - typedef typename OutputImageType::IndexType IndexType; - typedef typename Superclass::OutputImageRegionType OutputImageRegionType; - typedef itk::ImageRegionIteratorWithIndex< Image<PixelType,2, 0> > ImageIteratorType; - + typedef DEMReader Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + typedef itk::Point<double, 2> PointType; + /** Method for creation through the object factory. */ itkNewMacro(Self); /** Run-time type information (and related methods). */ - itkTypeMacro(DEMReader,ImageSource); + itkTypeMacro(DEMReader,Object); - /** Set the spacing. */ - itkSetMacro(Spacing,SpacingType); - itkGetConstReferenceMacro(Spacing,SpacingType); - - /** Set the Upper Left coordinates. */ - itkSetMacro(Ul,PointType); - itkGetConstReferenceMacro(Ul,PointType); - - /** Set the Lower Right coordinates. */ - itkSetMacro(Lr,PointType); - itkGetConstReferenceMacro(Lr,PointType); - - /** Set the spacing. */ - void SetSpacing(const double* spacing); - /** Try to open the DEM directory. */ - bool OpenDEMDirectory(char* &DEMDirectory); + bool OpenDEMDirectory(const char* DEMDirectory); - /** Compute the height above MSL(Mean Sea Level) of the point. */ - virtual double GetHeightAboveMSL(const PointType& worldPoint); + /** Compute the height above MSL(Mean Sea Level) of a geographic point. */ + virtual double GetHeightAboveMSL(const PointType& geoPoint); protected: DEMReader(); ~DEMReader(); void PrintSelf(std::ostream& os, Indent indent) const; - void GenerateData(); - virtual void GenerateOutputInformation(); ossimElevManager* m_ElevManager; - //DEMImagePointerType m_DEMImage; - SpacingType m_Spacing; - PointType m_Ul; - PointType m_Lr; + }; } // namespace otb -#ifndef OTB_MANUAL_INSTANTIATION -#include "otbDEMReader.txx" -#endif + #endif diff --git a/Code/IO/otbDEMReader.txx b/Code/IO/otbDEMReader.txx deleted file mode 100755 index 7c7bfd99a978eb4bbe13de88a7eb124685984a57..0000000000000000000000000000000000000000 --- a/Code/IO/otbDEMReader.txx +++ /dev/null @@ -1,152 +0,0 @@ -/*========================================================================= - - 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 __otbDEMReader_txx -#define __otbDEMReader_txx - -#include "otbDEMReader.h" -#include "otbMacro.h" - -namespace otb -{ - - -template<class TDEMImage> DEMReader<TDEMImage>::DEMReader() -{ - m_ElevManager=ossimElevManager::instance(); -//m_DEMImage = DEMImageType::New(); - m_Spacing[0]=0; - m_Spacing[1]=0; - m_Ul[0]=0; - m_Ul[1]=0; - m_Lr[0]=0; - m_Lr[1]=0; -} - -template<class TDEMImage> -DEMReader<TDEMImage>::~DEMReader() -{ - delete m_ElevManager; -} - -///Methode pour specifier un dossier contenant des DEM -template<class TDEMImage> bool DEMReader<TDEMImage>::OpenDEMDirectory(char* &DEMDirectory) -{ - ossimFilename ossimDEMDir; - ossimDEMDir=ossimFilename(DEMDirectory); - bool result= false; - if (m_ElevManager->openDirectory(ossimDEMDir)) - { - result= true; - } - return result; -} - -///Methode pour calculer l'altitude d'un point geographique -template<class TDEMImage> double DEMReader<TDEMImage>::GetHeightAboveMSL(const PointType& worldPoint) -{ - float height; - ossimGpt ossimWorldPoint; - ossimWorldPoint.lat=worldPoint[0]; - ossimWorldPoint.lon=worldPoint[1]; - height=m_ElevManager->getHeightAboveMSL(ossimWorldPoint); - return height; -} - -///Methode SetSpacing -template<class TDEMImage> void DEMReader<TDEMImage>::SetSpacing( const double* spacing) -{ - SpacingType s(spacing); - this->SetSpacing(s); -} - -///Methode GenerateOutputInformation -template <class TDEMImage> void DEMReader<TDEMImage>::GenerateOutputInformation() -{ - DEMImageType *output; - IndexType start; - start[0]=0; - start[1]=0; - - PointType origin; - origin[0]=m_Ul[0]; //latitude de l'origine. - origin[1]=m_Ul[1]; //longitude de l'origine. - output = this->GetOutput(0); - - // Calcul de la taille de l'image: - SizeType size; - size[0]=int (abs(((m_Lr[0]-m_Ul[0])/m_Spacing[0]))+1.5); - size[1]=int (abs(((m_Lr[1]-m_Ul[1])/m_Spacing[1]))+1.5); - - // On specifie les parametres de la region - OutputImageRegionType largestPossibleRegion; - largestPossibleRegion.SetSize( size ); - largestPossibleRegion.SetIndex( start ); - output->SetLargestPossibleRegion( largestPossibleRegion ); - - output->SetSpacing(m_Spacing); - output->SetOrigin(m_Ul); -} - -///Methode GenerateData -template <class TDEMImage> void DEMReader<TDEMImage>::GenerateData() -{ - DEMImagePointerType m_DEMImage = this->GetOutput(); - - // allocate the output buffer - m_DEMImage->SetBufferedRegion( m_DEMImage->GetRequestedRegion() ); - m_DEMImage->Allocate(); - - // Create an iterator that will walk the output region - ImageIteratorType outIt = ImageIteratorType(m_DEMImage,m_DEMImage->GetRequestedRegion()); - - // Walk the output image, evaluating the height at each pixel - IndexType currentindex; - PointType phyPoint; - double height; - for (outIt.GoToBegin(); !outIt.IsAtEnd(); ++outIt) - { - currentindex=outIt.GetIndex(); - m_DEMImage->TransformIndexToPhysicalPoint(currentindex, phyPoint); - ossimGpt ossimWorldPoint; - ossimWorldPoint.lat=phyPoint[0]; - ossimWorldPoint.lon=phyPoint[1]; - height=m_ElevManager->getHeightAboveMSL(ossimWorldPoint); //Calcul de l'altitude - otbMsgDebugMacro(<<" HeightAboveMSL: "<<height); - if (height>-static_cast<double>(32768)) //On teste si les fichiers MNT recouvre la zone g�ographique demand�e (-32768 = theNullHeightValue) - { - m_DEMImage->SetPixel(currentindex, static_cast<PixelType>(height) ); - } //On remplit l'image - else - { - m_DEMImage->SetPixel(currentindex, static_cast<PixelType>(0) ); - } - } -} - -template <class TDEMImage> void DEMReader<TDEMImage>::PrintSelf(std::ostream& os, Indent indent) const -{ - Superclass::PrintSelf(os,indent); - - os << indent << "Spacing:"<< m_Spacing[0] << ","<< m_Spacing[1] << std::endl; - os << indent << "Lr:"<< m_Lr[0] << ","<< m_Lr[1] << std::endl; - os << indent << "Ul:"<< m_Ul[0] << ","<< m_Ul[1] << std::endl; -} - -} // namespace otb - -#endif diff --git a/Code/IO/otbGDALImageIO.cxx b/Code/IO/otbGDALImageIO.cxx index 9bd835e25ce9b3340b085af29fcb806898b9f76f..5b27bf944d0761ba5be98abcef7de9442950adf3 100755 --- a/Code/IO/otbGDALImageIO.cxx +++ b/Code/IO/otbGDALImageIO.cxx @@ -593,8 +593,8 @@ void GDALImageIO::InternalReadImageInformation() } - m_Origin[0] = minGCPX; - m_Origin[1] = minGCPY; + m_Origin[0] = minGCPX; + m_Origin[1] = minGCPY; } @@ -996,27 +996,41 @@ void GDALImageIO::InternalWriteImageInformation() /* Set the GCPs */ /* -------------------------------------------------------------------- */ - - if(ImageBase::GetGCPCount(dico)>0) - { - unsigned int gcpCount = ImageBase::GetGCPCount(dico); - GDAL_GCP * gdalGcps = new GDAL_GCP[gcpCount]; - - for(unsigned int gcpIndex = 0; gcpIndex < gcpCount;++gcpIndex) - { - gdalGcps[gcpIndex].pszId = const_cast<char *>(ImageBase::GetGCPId(dico,gcpIndex).c_str()); - gdalGcps[gcpIndex].pszInfo = const_cast<char *>(ImageBase::GetGCPInfo(dico,gcpIndex).c_str()); - gdalGcps[gcpIndex].dfGCPPixel = ImageBase::GetGCPCol(dico,gcpIndex); - gdalGcps[gcpIndex].dfGCPLine = ImageBase::GetGCPRow(dico,gcpIndex); - gdalGcps[gcpIndex].dfGCPX = ImageBase::GetGCPX(dico,gcpIndex); - gdalGcps[gcpIndex].dfGCPY = ImageBase::GetGCPY(dico,gcpIndex); - gdalGcps[gcpIndex].dfGCPZ = ImageBase::GetGCPZ(dico,gcpIndex); - } - + unsigned int gcpCount = ImageBase::GetGCPCount(dico); + + GDAL_GCP * gdalGcps = new GDAL_GCP[gcpCount+1]; + + bool gcpHasOrigin = false; + for(unsigned int gcpIndex = 0; gcpIndex < gcpCount;++gcpIndex) + { + gdalGcps[gcpIndex].pszId = const_cast<char *>(ImageBase::GetGCPId(dico,gcpIndex).c_str()); + gdalGcps[gcpIndex].pszInfo = const_cast<char *>(ImageBase::GetGCPInfo(dico,gcpIndex).c_str()); + gdalGcps[gcpIndex].dfGCPPixel = ImageBase::GetGCPCol(dico,gcpIndex); + gdalGcps[gcpIndex].dfGCPLine = ImageBase::GetGCPRow(dico,gcpIndex); + gdalGcps[gcpIndex].dfGCPX = ImageBase::GetGCPX(dico,gcpIndex); + gdalGcps[gcpIndex].dfGCPY = ImageBase::GetGCPY(dico,gcpIndex); + gdalGcps[gcpIndex].dfGCPZ = ImageBase::GetGCPZ(dico,gcpIndex); + gcpHasOrigin = ImageBase::GetGCPCol(dico,gcpIndex)==0 && ImageBase::GetGCPRow(dico,gcpIndex)==0; + } + gdalGcps[gcpCount].pszId = "Origin"; + gdalGcps[gcpCount].pszInfo = "Origin gcp added by OTB"; + gdalGcps[gcpCount].dfGCPPixel = 0; + gdalGcps[gcpCount].dfGCPLine = 0; + gdalGcps[gcpCount].dfGCPX = m_Origin[0]; + gdalGcps[gcpCount].dfGCPY = m_Origin[1]; + gdalGcps[gcpCount].dfGCPZ = 0; + + if(gcpHasOrigin) + { m_poDataset->SetGCPs(gcpCount,gdalGcps,ImageBase::GetGCPProjection(dico).c_str()); - delete [] gdalGcps; } + else + { + std::cout<<"GCPs do not contain origin."<<std::endl; + m_poDataset->SetGCPs(gcpCount+1,gdalGcps,ImageBase::GetGCPProjection(dico).c_str()); + } + delete [] gdalGcps; /* -------------------------------------------------------------------- */ /* Set the six coefficients of affine geoTtransform */ diff --git a/Code/IO/otbStreamingImageFileWriter.txx b/Code/IO/otbStreamingImageFileWriter.txx index 2fef00b833adbfc3a6389d6a8f8fd467274da3d4..c8677cb34a1c61e342ae7ac5f4d057e8e1698f12 100755 --- a/Code/IO/otbStreamingImageFileWriter.txx +++ b/Code/IO/otbStreamingImageFileWriter.txx @@ -467,6 +467,8 @@ StreamingImageFileWriter<TInputImage> * Loop over the number of pieces, execute the upstream pipeline on each * piece, and copy the results into the output image. */ + otbGenericMsgDebugMacro(<< "Number Of Stream Divisions : " << numDivisionsFromSplitter); + unsigned int piece; for (piece = 0; piece < numDivisionsFromSplitter && !this->GetAbortGenerateData(); diff --git a/Code/Projections/otbCompositeTransform.h b/Code/Projections/otbCompositeTransform.h index f75aa79ef7051ade64d99a91aa45c96f6f0a550f..bd13645b73894ac133dde9f0e4c87e94f7699a5b 100644 --- a/Code/Projections/otbCompositeTransform.h +++ b/Code/Projections/otbCompositeTransform.h @@ -62,20 +62,24 @@ public : typedef typename TSecondTransform::Pointer SecondTransformPointerType; /** Standard vector type for this class. */ - typedef typename Superclass::InputVectorType InputVectorType; - typedef typename Superclass::OutputVectorType OutputVectorType; +// typedef typename TFirstTransform::InputVectorType FirstTransformInputVectorType; +// typedef typename TFirstTransform::OutputVectorType FirstTransformOutputVectorType; /** Standard covariant vector type for this class */ - typedef typename Superclass::InputCovariantVectorType InputCovariantVectorType; - typedef typename Superclass::OutputCovariantVectorType OutputCovariantVectorType; +// typedef typename TFirstTransform::InputCovariantVectorType FirstTransformInputCovariantVectorType; +// typedef typename TFirstTransform::OutputCovariantVectorType FirstTransformOutputCovariantVectorType; /** Standard vnl_vector type for this class. */ - typedef typename Superclass::InputVnlVectorType InputVnlVectorType; - typedef typename Superclass::OutputVnlVectorType OutputVnlVectorType; +// typedef typename TFirstTransform::InputVnlVectorType FirstTransformInputVnlVectorType; +// typedef typename TFirstTransform::OutputVnlVectorType FirstTransformOutputVnlVectorType; /** Standard coordinate point type for this class */ - typedef typename Superclass::InputPointType InputPointType; - typedef typename Superclass::OutputPointType OutputPointType; + typedef typename Superclass::InputPointType FirstTransformInputPointType; + typedef typename TFirstTransform::OutputPointType FirstTransformOutputPointType; + /** TSecondTransform::InputPointType is purposely not defined + * It contrains user to choose First Transform Output compatible + * with Second Transform Input */ + typedef typename Superclass::OutputPointType SecondTransformOutputPointType; /** Method for creation through the object factory. */ itkNewMacro( Self ); @@ -96,22 +100,16 @@ public : /** Method to transform a point. */ - virtual OutputPointType TransformPoint(const InputPointType & ) const; + virtual SecondTransformOutputPointType TransformPoint(const FirstTransformInputPointType & ) const; /** Method to transform a vector. */ - virtual OutputVectorType TransformVector(const InputVectorType &) const; +// virtual OutputVectorType TransformVector(const InputVectorType &) const; /** Method to transform a vnl_vector. */ - virtual OutputVnlVectorType TransformVector(const InputVnlVectorType &) const; +// virtual OutputVnlVectorType TransformVector(const InputVnlVectorType &) const; /** Method to transform a CovariantVector. */ - virtual OutputCovariantVectorType TransformCovariantVector(const InputCovariantVectorType &) const; - - /** Compute MapProjection1 coordinates to MapProjection2 coordinates. */ - OutputPointType ComputeProjection1ToProjection2(const InputPointType &point1); - - /** Compute MapProjection1 coordinates to MapProjection2 coordinates. */ - InputPointType ComputeProjection2ToProjection1(const OutputPointType &point2); +// virtual OutputCovariantVectorType TransformCovariantVector(const InputCovariantVectorType &) const; protected: CompositeTransform(); diff --git a/Code/Projections/otbCompositeTransform.txx b/Code/Projections/otbCompositeTransform.txx index 5ea6d7dd1f08cca78cdf4cca58460e114f8d0a85..4566dec01139e107009d3839e249f98fd8b50ccd 100644 --- a/Code/Projections/otbCompositeTransform.txx +++ b/Code/Projections/otbCompositeTransform.txx @@ -19,6 +19,7 @@ #define __otbCompositeTransform_txx #include "otbCompositeTransform.h" +#include "otbMapProjections.h" namespace otb { @@ -39,20 +40,25 @@ CompositeTransform<TFirstTransform, TSecondTransform, TScalarType, NInputDimensi template<class TFirstTransform, class TSecondTransform, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -typename CompositeTransform<TFirstTransform, TSecondTransform, TScalarType, NInputDimensions, NOutputDimensions>::OutputPointType +typename CompositeTransform<TFirstTransform, TSecondTransform, TScalarType,NInputDimensions, NOutputDimensions>::SecondTransformOutputPointType CompositeTransform<TFirstTransform, TSecondTransform, TScalarType, NInputDimensions, NOutputDimensions> -::TransformPoint(const InputPointType &point1) const +::TransformPoint(const FirstTransformInputPointType &point1) const { - InputPointType pointTmp; - OutputPointType point2; + FirstTransformOutputPointType pointTmp; + SecondTransformOutputPointType point2; pointTmp=m_FirstTransform->TransformPoint(point1); + //pointTmp=utmprojection->TransformPoint(point1); + otbMsgDevMacro(<< "Le point g�ographique correspondant est: ("<< pointTmp[0]<< ","<< pointTmp[1]<< ")"); + point2=m_SecondTransform->TransformPoint(pointTmp); + otbMsgDevMacro(<< "Les coordonn�es en pixel sur l'image capteur correspondant � ce point sont:" << std::endl + << point2[0] << ","<< point2[1] ); return point2; } -template<class TFirstTransform, class TSecondTransform, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> +/*template<class TFirstTransform, class TSecondTransform, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> typename CompositeTransform<TFirstTransform, TSecondTransform, TScalarType, NInputDimensions, NOutputDimensions>::OutputVectorType CompositeTransform<TFirstTransform, TSecondTransform, TScalarType,NInputDimensions, NOutputDimensions> ::TransformVector(const InputVectorType &vector1) const @@ -92,7 +98,7 @@ CompositeTransform<TFirstTransform, TSecondTransform, TScalarType,NInputDimensio covariantVector2=m_SecondTransform->TransformCovariantVector(covariantVectorTmp); return covariantVector2; -} +}*/ diff --git a/Code/Projections/otbInverseSensorModel.h b/Code/Projections/otbInverseSensorModel.h index d54820ffdc9a4fd468ceff7a13076a1afc8fe43a..9f8426f32e476787da57480a06134bd4d12cacaa 100644 --- a/Code/Projections/otbInverseSensorModel.h +++ b/Code/Projections/otbInverseSensorModel.h @@ -19,6 +19,8 @@ #define __otbInverseSensorModel_h #include "otbSensorModelBase.h" +#include "otbDEMReader.h" + #include "itkMacro.h" #include "itkSmartPointer.h" #include "itkObject.h" @@ -51,16 +53,21 @@ class ITK_EXPORT InverseSensorModel : public SensorModelBase<TScalarType, public : /** Standard class typedefs. */ - typedef InverseSensorModel Self; + typedef InverseSensorModel Self; typedef SensorModelBase< TScalarType, NInputDimensions, NOutputDimensions, - NParametersDimensions > Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; + NParametersDimensions > Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; typedef typename Superclass::InputPointType InputPointType; - typedef typename Superclass::OutputPointType OutputPointType; +// typedef itk::Point<TScalarType, 3> InputPointType; + typedef typename Superclass::OutputPointType OutputPointType; + + //typedef otb::Image<double, NInputDimensions> ImageType; + typedef DEMReader DEMHandlerType; + typedef typename DEMHandlerType::Pointer DEMHandlerPointerType; /** Method for creation through the object factory. */ itkNewMacro( Self ); @@ -77,6 +84,30 @@ public : // Pour projeter un point g�o connaissant son altitude. OutputPointType TransformPoint(const InputPointType &point, double height) const; + itkGetMacro(UseDEM, bool); + itkSetMacro(UseDEM, bool); + + itkGetObjectMacro(DEMHandler, DEMHandlerType); + + virtual void SetDEMHandler(DEMHandlerType* _arg) + { + if (this->m_DEMHandler != _arg) + { + this->m_DEMHandler = _arg; + this->Modified(); + m_UseDEM = true; + } + } + + virtual void SetDEMDirectory(const std::string& directory) + { + m_DEMHandler->OpenDEMDirectory(directory.c_str()); + m_UseDEM = true; + } + + void ActiveDEM() { m_UseDEM = true; this->Modified(); } + void DesactiveDEM() { m_UseDEM = false; this->Modified();} + protected: InverseSensorModel(); virtual ~InverseSensorModel(); @@ -88,6 +119,13 @@ private : InverseSensorModel(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented + + /** Object that read and use DEM */ + DEMHandlerPointerType m_DEMHandler; + + /** Specify if DEM is used in Point Transformation */ + bool m_UseDEM ; + }; diff --git a/Code/Projections/otbInverseSensorModel.txx b/Code/Projections/otbInverseSensorModel.txx index f24a6511db388209f0b858dd46bc0b30d3c9f659..b2c8f02e4cdf6d6ed1037f84ad92956ef734bc7e 100644 --- a/Code/Projections/otbInverseSensorModel.txx +++ b/Code/Projections/otbInverseSensorModel.txx @@ -33,6 +33,8 @@ template < class TScalarType, InverseSensorModel< TScalarType,NInputDimensions,NOutputDimensions,NParametersDimensions> ::InverseSensorModel() { + m_DEMHandler = DEMHandlerType::New(); + m_UseDEM = false; } template < class TScalarType, @@ -52,12 +54,13 @@ template < class TScalarType, unsigned int NParametersDimensions > typename InverseSensorModel< TScalarType,NInputDimensions,NOutputDimensions,NParametersDimensions>::OutputPointType InverseSensorModel< TScalarType,NInputDimensions,NOutputDimensions,NParametersDimensions> -::TransformPoint(const InputPointType &point) const +::TransformPoint(const InputPointType &point, double height) const { - +// std::cout << "INVERSESENSORMODEL MIARY" << std::endl; // std::cout << "InverseSensorModel Before : (" << point[0] << "," << point[1] << ")" << std::endl; + //On transforme le type "itk::point" en type "ossim::ossimGpt" - ossimGpt ossimGPoint(point[0], point[1]); + ossimGpt ossimGPoint(point[0], point[1], height); //On calcule ossimDpt ossimDPoint; @@ -85,10 +88,20 @@ template < class TScalarType, unsigned int NParametersDimensions > typename InverseSensorModel< TScalarType,NInputDimensions,NOutputDimensions,NParametersDimensions>::OutputPointType InverseSensorModel< TScalarType,NInputDimensions,NOutputDimensions,NParametersDimensions> -::TransformPoint(const InputPointType &point, double height) const +::TransformPoint(const InputPointType &point) const { //On transforme le type "itk::point" en type "ossim::ossimGpt" - ossimGpt ossimGPoint(point[0], point[1], height); + ossimGpt ossimGPoint(point[0], point[1]); + + if (m_UseDEM) + { + otbMsgDebugMacro(<< "USING DEM ! ") ; + double height = m_DEMHandler->GetHeightAboveMSL(point); +// ossimGpt ossimGPoint(point[0], point[1], height); + otbMsgDebugMacro(<< "height : " << height) ; + ossimGPoint.height(height); +// ossimGPoint.height(height); + } //On calcule ossimDpt ossimDPoint; @@ -103,7 +116,7 @@ InverseSensorModel< TScalarType,NInputDimensions,NOutputDimensions,NParametersDi OutputPointType outputPoint; outputPoint[0]=ossimDPoint.x; outputPoint[1]=ossimDPoint.y; - + return outputPoint; } diff --git a/Code/Projections/otbMapProjection.h b/Code/Projections/otbMapProjection.h index afbe6450ee42057930223bf47dc167b6ed928a35..788fbe1d9a45bb45da7d429da9958c31263ec9f9 100644 --- a/Code/Projections/otbMapProjection.h +++ b/Code/Projections/otbMapProjection.h @@ -43,10 +43,13 @@ namespace otb * Map projection a sa propre classe. **/ +typedef enum {FORWARD=0, INVERSE=1} InverseOrForwardTransformationEnum; + template <class TOssimMapProjection, class TScalarType = double, unsigned int NInputDimensions=2, - unsigned int NOutputDimensions=2> + unsigned int NOutputDimensions=2, + InverseOrForwardTransformationEnum transform=INVERSE> class ITK_EXPORT MapProjection: public itk::Transform<TScalarType, // Data type for scalars NInputDimensions, // Number of dimensions in the input space NOutputDimensions> // Number of dimensions in the output space @@ -89,9 +92,11 @@ public : void SetEllipsoid(const double &major_axis, const double &minor_axis); - OutputPointType Forward(const InputPointType &point) const; +// OutputPointType Forward(const InputPointType &point) const; - InputPointType Inverse(const OutputPointType &point) const; +// InputPointType Inverse(const OutputPointType &point) const; + + OutputPointType TransformPoint(const InputPointType &point) const; virtual InputPointType Origin(); diff --git a/Code/Projections/otbMapProjection.txx b/Code/Projections/otbMapProjection.txx index a5cfd607f3ec08aec411a3f6c8af46a92a8272b3..140a84e570a6c5ebf66fe1032ce205d561102e7f 100644 --- a/Code/Projections/otbMapProjection.txx +++ b/Code/Projections/otbMapProjection.txx @@ -23,16 +23,16 @@ namespace otb { -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::MapProjection() : Superclass(SpaceDimension,ParametersDimension) { m_MapProjection = new OssimMapProjectionType(); } -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::~MapProjection() { delete m_MapProjection; @@ -41,8 +41,8 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi ///M�thode pour attribuer une ellipse � la projection(4 m�thodes) //Par defaut: -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::SetEllipsoid() { ossimEllipsoid ellipsoid=ossimEllipsoid(); @@ -50,16 +50,16 @@ void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDi } //Par copie: -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::SetEllipsoid(const ossimEllipsoid &ellipsoid) { m_MapProjection->setEllipsoid(ellipsoid); } //En connaissant son "code" -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::SetEllipsoid(std::string code) { const ossimEllipsoid *ellipsoid = ossimEllipsoidFactory::instance()->create(ossimString(code)); @@ -76,8 +76,8 @@ void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDi } //En connaissant ses axes -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::SetEllipsoid(const double &major_axis, const double &minor_axis) { ossimEllipsoid ellipse(major_axis,minor_axis); @@ -88,9 +88,9 @@ void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDi ///M�thode qui transforme un point g�ographique en un point cartographique: //Encapsulation de la m�thode "forward" -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -typename MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions>::OutputPointType -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +/*template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +typename MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform>::OutputPointType +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::Forward(const InputPointType & point) const { @@ -106,16 +106,16 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi outputPoint[1]=ossimDPoint.y; return outputPoint; -} +}*/ ///M�thode qui transforme un point cartographique en un point g�ographique: //Encapsulation de la m�thode "inverse" -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -typename MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions>::InputPointType -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +/*template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +typename MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform>::InputPointType +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::Inverse(const OutputPointType & point) const { - //On transforme le type "itk::point" en type "ossim::ossimGpt" + //On transforme le type "itk::point" en type "ossim::ossimDpt" ossimDpt ossimDPoint(point[0], point[1]); //On le proj�te sur la carte @@ -127,12 +127,53 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi inputPoint[1]=ossimGPoint.lon; return inputPoint; +}*/ + +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +typename MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform>::OutputPointType +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> +::TransformPoint(const InputPointType & point) const +{ + OutputPointType outputPoint; + + if (Transform == INVERSE) + { + //On transforme le type "itk::point" en type "ossim::ossimDpt" + ossimDpt ossimDPoint(point[0], point[1]); + + //On le proj�te sur la carte + ossimGpt ossimGPoint; + ossimGPoint=m_MapProjection->inverse(ossimDPoint); + + outputPoint[0]=ossimGPoint.lat; + outputPoint[1]=ossimGPoint.lon; + } + else if (Transform == FORWARD) + { + //On transforme le type "itk::point" en type "ossim::ossimGpt" + ossimGpt ossimGPoint(point[0], point[1]); + + //On le proj�te sur la carte + ossimDpt ossimDPoint; + ossimDPoint=m_MapProjection->forward(ossimGPoint); + + outputPoint[0]=ossimDPoint.x; + outputPoint[1]=ossimDPoint.y; + } + else + { + itkExceptionMacro(<<"Model is INVERSE or FORWARD only !!"); + } + + return outputPoint; } + + ///M�thode qui retourne le point g�ographique correspondant � l'index (0;0) -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -typename MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions>::InputPointType -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +typename MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform>::InputPointType +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::Origin() { ossimGpt ossimOrigin=m_MapProjection->origin(); @@ -145,9 +186,9 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi } ///M�thode pour r�cup�rer le FalseNorthing(pour �viter les coordonn�es n�gatives) -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> double -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::GetFalseNorthing() const { double falseNorthing=m_MapProjection->getFalseNorthing(); @@ -156,9 +197,9 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi } ///M�thode pour r�cup�rer le FalseEasting(pour �viter les coordonn�es n�gatives) -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> double -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::GetFalseEasting() const { double falseEasting=m_MapProjection->getFalseEasting(); @@ -167,9 +208,9 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi } ///M�thode pour r�cup�rer le StandardParallel1(# selon le type de projection) -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> double -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::GetStandardParallel1() const { double standardParallel1=m_MapProjection->getStandardParallel1(); @@ -178,9 +219,9 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi } ///M�thode pour r�cup�rer le StandardParallel2(# selon le type de projection) -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> double -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::GetStandardParallel2() const { double standardParallel2=m_MapProjection->getStandardParallel2(); @@ -189,9 +230,9 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi } ///M�thode pour r�cup�rer le nom de la projection -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> std::string -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::GetProjectionName() const { std::string projectionName; @@ -201,18 +242,18 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi } ///M�thode pour v�rifier si... -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> bool -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::IsGeographic() const { return (m_MapProjection->isGeographic()); } ///M�thode pour r�cup�rer le "major axis" d'une ellipse -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> double -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::GetA() const { double majorAxis=m_MapProjection->getA(); @@ -221,9 +262,9 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi } ///M�thode pour r�cup�rer le "minor axis" d'une ellipse -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> double -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::GetB() const { double minorAxis=m_MapProjection->getB(); @@ -232,9 +273,9 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi } ///M�thode pour r�cup�rer le "Flattening" d'une ellipse -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> double -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::GetF() const { double flattening=m_MapProjection->getF(); @@ -243,9 +284,9 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi } ///M�thode qui retourne la r�solution (m�tres/pixel) -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -typename MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions>::OutputPointType -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +typename MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform>::OutputPointType +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::GetMetersPerPixel() const { ossimDpt ossimMetersPerPixels=m_MapProjection->getMetersPerPixel(); @@ -258,9 +299,9 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi } ///M�thode qui retourne la r�solution (degr�s/pixel) -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -typename MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions>::OutputPointType -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +typename MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform>::OutputPointType +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::GetDecimalDegreesPerPixel() const { ossimDpt ossimDecimalDegreesPerPixels=m_MapProjection->getDecimalDegreesPerPixel(); @@ -273,16 +314,16 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi } ///M�thode pour attribuer les axes � l'ellipse -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::SetAB(double a, double b) { m_MapProjection->setAB(a,b); } ///M�thode pour attribuer l'origine -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::SetOrigin(const InputPointType &origin) { ossimGpt ossimOrigin(origin[0], origin[1]); @@ -290,16 +331,16 @@ void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDi } ///M�thode pour fixer une r�solution de la carte -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::SetMetersPerPixel(const OutputPointType &point) {ossimDpt ossimDPoint(point[0], point[1]); m_MapProjection->setMetersPerPixel(ossimDPoint); } ///M�thode pour fixer une r�solution de la carte -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::SetDecimalDegreesPerPixel(const OutputPointType &point) { ossimDpt ossimDPoint(point[0], point[1]); @@ -307,8 +348,8 @@ void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDi } ///M�thode pour calculer approximativement la r�solution(degr�/pixel) -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> -void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> +void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::ComputeDegreesPerPixel(const InputPointType &ground, const OutputPointType &metersPerPixel, double &deltaLat, double &deltaLon) { ossimDpt ossimMetersPerPixel(metersPerPixel[0], metersPerPixel[1]); @@ -318,9 +359,9 @@ void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDi } ///M�thode pour calculer approximativement la r�solution(m�tre/pixel) -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> void -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::ComputeMetersPerPixel(const InputPointType ¢er, double deltaDegreesPerPixelLat, double deltaDegreesPerPixelLon, OutputPointType &metersPerPixel) { //Correction @@ -332,9 +373,9 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi } ///M�thode pour calculer approximativement la r�solution(m�tre/pixel) -template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> +template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> void -MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions> +MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, Transform> ::ComputeMetersPerPixel(double deltaDegreesPerPixelLat, double deltaDegreesPerPixelLon, OutputPointType &metersPerPixel) { ossimDpt ossimMetersPerPixel; @@ -344,7 +385,7 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi } ///Non impl�ment� dans ossim. -// template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions>::SetMatrix(double rotation, const OutputPointType &scale, const OutputPointType &translation) +// template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, InverseOrForwardTransformationEnum>::SetMatrix(double rotation, const OutputPointType &scale, const OutputPointType &translation) // { // ossimDpt ossimScale(scale[0], scale[1]); // ossimDpt ossimTranslation(translation[0], translation[1]); @@ -352,7 +393,7 @@ MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensi // } ///M�thode pour attribuer un falseEasting -// template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions> void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions>::SetFalseEasting(double falseEasting) +// template<class TOssimMapProjection, class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions, InverseOrForwardTransformationEnum Transform> void MapProjection<TOssimMapProjection, TScalarType, NInputDimensions, NOutputDimensions, InverseOrForwardTransformationEnum>::SetFalseEasting(double falseEasting) // { std::string projectionName=m_MapProjection->getProjectionName(); // if (projectionName=="ossimUtmProjection") // {m_MapProjection->setFalseEasting(falseEasting); diff --git a/Code/Projections/otbOrthoRectificationFilter.h b/Code/Projections/otbOrthoRectificationFilter.h index b52b7441b323d2ca15c75ce0a0df62cbbf1fac49..3f322bca8778b23457c8b76be4498de9872ae8be 100644 --- a/Code/Projections/otbOrthoRectificationFilter.h +++ b/Code/Projections/otbOrthoRectificationFilter.h @@ -25,20 +25,12 @@ #include "otbMapProjection.h" #include "otbStreamingResampleImageFilter.h" #include "otbCompositeTransform.h" +#include "otbInverseSensorModel.h" namespace otb { -/** Macro used to access Resampler class */ -/*#define otbSetResampleMacro(name, type)\ - virtual void SetResample##name( const type & _arg)\ - {\ - if (m_Resampler->Get##name() != _arg)\ - {\ - m_Resampler->Set##name(_arg);\ - this->Modified();\ - }\ - }\*/ + /**Template otbOrthoRectificationFilter.txx * Cette classe permet de passer d'une MapProjection � une autre en passant par un point g�ographique. @@ -59,9 +51,6 @@ public : typedef itk::SmartPointer<Self> Pointer; typedef itk::SmartPointer<const Self> ConstPointer; - -// typedef otb::StreamingResampleImageFilter< TInputImage, TOutputImage > ResamplerType; -// typedef typename ResamplerType::Pointer ResamplerPointerType; typedef typename TInputImage::IndexType IndexType; typedef typename TInputImage::SizeType SizeType; typedef typename TInputImage::SpacingType SpacingType; @@ -76,12 +65,9 @@ public : typedef InverseSensorModel<double> SensorModelType; typedef typename SensorModelType::Pointer SensorModelPointerType; - typedef CompositeTransform<SensorModelType,MapProjectionType> CompositeTransformType; + typedef CompositeTransform< MapProjectionType,SensorModelType> CompositeTransformType; typedef typename CompositeTransformType::Pointer CompositeTransformPointerType; - typedef typename Superclass::InterpolatorType InterpolatorType; -// typedef typename InterpolatorType::Pointer InterpolatorPointerType; - /** Method for creation through the object factory. */ itkNewMacro( Self ); @@ -89,8 +75,35 @@ public : itkTypeMacro( OrthoRectificationFilter, StreamingResampleImageFilter ); /** Accessors */ - itkSetMacro(MapProjection, MapProjectionPointerType); - itkGetMacro(MapProjection, MapProjectionPointerType); + virtual void SetMapProjection (MapProjectionType* _arg) + { + if (this->m_MapProjection != _arg) + { + this->m_MapProjection = _arg; + m_CompositeTransform->SetFirstTransform(_arg); + m_IsComputed = false; + this->Modified(); + } + } + + itkGetObjectMacro(MapProjection, MapProjectionType); + + virtual void SetDEMDirectory(const std::string& directory) + { + m_SensorModel->SetDEMDirectory(directory); + this->Modified(); + } + + virtual void ActiveDEM() + { + m_SensorModel->ActiveDEM(); + } + + virtual void DesactiveDEM() + { + m_SensorModel->DesactiveDEM(); + } + protected: OrthoRectificationFilter(); @@ -108,10 +121,10 @@ private: /** Calculate transformation model from sensor model & map projection composition */ void ComputeResampleTransformationModel(); - void Modified(); /** Boolean used to know if transformation model computation is needed */ bool m_IsComputed; + SensorModelPointerType m_SensorModel; MapProjectionPointerType m_MapProjection; CompositeTransformPointerType m_CompositeTransform; diff --git a/Code/Projections/otbOrthoRectificationFilter.txx b/Code/Projections/otbOrthoRectificationFilter.txx index 2ff358482c8bcdcf82f26f16dfec82bc4569583a..29a71bc5a5b3a130ab73c740884c2c9497ebfd38 100644 --- a/Code/Projections/otbOrthoRectificationFilter.txx +++ b/Code/Projections/otbOrthoRectificationFilter.txx @@ -90,13 +90,14 @@ OrthoRectificationFilter<TInputImage, TOutputImage, TMapProjection, TInterpolato { if (m_IsComputed == false) { + otbMsgDevMacro(<< "COMPUTE RESAMPLE TRANSFORMATION MODEL"); typename TOutputImage::Pointer output = this->GetOutput(); // Get OSSIM sensor model from image keywordlist m_SensorModel->SetImageGeometry(output->GetImageKeywordlist()); - m_CompositeTransform->SetFirstTransform(m_SensorModel); - m_CompositeTransform->SetSecondTransform(m_MapProjection); + m_CompositeTransform->SetFirstTransform(m_MapProjection); + m_CompositeTransform->SetSecondTransform(m_SensorModel); this->SetTransform(m_CompositeTransform); @@ -104,15 +105,6 @@ OrthoRectificationFilter<TInputImage, TOutputImage, TMapProjection, TInterpolato } } -template <class TInputImage, class TOutputImage, class TMapProjection, class TInterpolatorPrecision> -void -OrthoRectificationFilter<TInputImage, TOutputImage, TMapProjection, TInterpolatorPrecision> -::Modified() -{ - m_IsComputed = false; - this->Modified(); -} - } //namespace otb #endif diff --git a/Code/Projections/otbSensorModelBase.h b/Code/Projections/otbSensorModelBase.h index eef7b278bd1ed365aeaad44a1426cf832736344d..6b6d24a23dfb45573e18adafdacb181acf3ba023 100755 --- a/Code/Projections/otbSensorModelBase.h +++ b/Code/Projections/otbSensorModelBase.h @@ -31,7 +31,7 @@ namespace otb { template <class TScalarType, - unsigned int NInputDimensions=2, + unsigned int NInputDimensions=3, unsigned int NOutputDimensions=2, unsigned int NParametersDimensions=3> class ITK_EXPORT SensorModelBase : public itk::Transform<TScalarType, diff --git a/Testing/Code/IO/CMakeLists.txt b/Testing/Code/IO/CMakeLists.txt index 1258a79af0998c9cafc36a6fc03d9e0c58ed583d..0677bbede969161294a14d67bf3d62f5137a3031 100755 --- a/Testing/Code/IO/CMakeLists.txt +++ b/Testing/Code/IO/CMakeLists.txt @@ -1110,7 +1110,7 @@ ADD_TEST(ioTuDEMReaderNew ${IO_TESTS} otbDEMReaderNew ) ADD_TEST(ioTvDEMReader ${IO_TESTS} --compare-ascii ${TOL} ${BASELINE_FILES}/ioDEMGetHeightAboveMSL.txt ${TEMP}/ioDEMGetHeightAboveMSL.txt - otbDEMReader + otbDEMReaderTest ${INPUTDATA}/DEM/srtm_directory ${TEMP}/ioDEMGetHeightAboveMSL.txt ) @@ -1183,7 +1183,7 @@ otbSpatialObjectDXFReaderNew.cxx otbSpatialObjectDXFReader.cxx otbOSSIMImageMetaDataReaderTest.cxx otbDEMReaderNew.cxx -otbDEMReader.cxx +otbDEMReaderTest.cxx otbImageGeometryHandlerNew.cxx ) diff --git a/Testing/Code/IO/otbDEMReaderNew.cxx b/Testing/Code/IO/otbDEMReaderNew.cxx index 52962cd779b9c6f80a1c234aae11438a4aabb8be..6a99e94c5e3bd058feb48dce07013c89def3e63d 100644 --- a/Testing/Code/IO/otbDEMReaderNew.cxx +++ b/Testing/Code/IO/otbDEMReaderNew.cxx @@ -26,7 +26,7 @@ int otbDEMReaderNew(int argc, char * argv[]) { const unsigned int Dimension = 2; typedef otb::Image<unsigned char,Dimension> ImageType; - typedef otb::DEMReader<ImageType> DEMReaderType; + typedef otb::DEMReader DEMReaderType; // Instantiating object DEMReaderType::Pointer object = DEMReaderType::New(); diff --git a/Testing/Code/IO/otbDEMReader.cxx b/Testing/Code/IO/otbDEMReaderTest.cxx similarity index 95% rename from Testing/Code/IO/otbDEMReader.cxx rename to Testing/Code/IO/otbDEMReaderTest.cxx index 725a99d75d0e617fbf264092c34e36d3317ad1f9..de93d021be20b3a22a9a9ff04c49595ca9d295aa 100644 --- a/Testing/Code/IO/otbDEMReader.cxx +++ b/Testing/Code/IO/otbDEMReaderTest.cxx @@ -23,7 +23,7 @@ #include "otbImageFileWriter.h" #include "otbMapProjections.h" -int otbDEMReader(int argc, char * argv[]) +int otbDEMReaderTest(int argc, char * argv[]) { try { @@ -34,7 +34,7 @@ int otbDEMReader(int argc, char * argv[]) bool bOpenDirectory(false); typedef otb::Image<float,Dimension> ImageType; - typedef otb::DEMReader<ImageType> DEMReaderType; + typedef otb::DEMReader DEMReaderType; // Instantiating object DEMReaderType::Pointer demReader = DEMReaderType::New(); diff --git a/Testing/Code/IO/otbIOTests.cxx b/Testing/Code/IO/otbIOTests.cxx index 1490fe7b0fe2b3d965fdbfb635aa58cdc16a4046..b4c97868ba63abaa887155150589e5d8cfa58f3a 100755 --- a/Testing/Code/IO/otbIOTests.cxx +++ b/Testing/Code/IO/otbIOTests.cxx @@ -75,6 +75,6 @@ REGISTER_TEST(otbSpatialObjectDXFReader); REGISTER_TEST(otbOSSIMImageMetaDataReaderTest); REGISTER_TEST(otbDEMReaderNew); -REGISTER_TEST(otbDEMReader); +REGISTER_TEST(otbDEMReaderTest); REGISTER_TEST(otbImageGeometryHandlerNew); } diff --git a/Testing/Code/Projections/CMakeLists.txt b/Testing/Code/Projections/CMakeLists.txt index 4a60abd88888065633ea002039922da193efa6ad..b90146aecf794863a509f31173a56344f5c0a955 100755 --- a/Testing/Code/Projections/CMakeLists.txt +++ b/Testing/Code/Projections/CMakeLists.txt @@ -15,135 +15,186 @@ SET(EPSILON 0.0000000001) SET(PROJECTIONS_TESTS ${CXX_TEST_PATH}/otbProjectionsTests) - +# TEST 1 ADD_TEST(prTuProjectionBaseNew ${PROJECTIONS_TESTS} otbProjectionBaseNew ) +# TEST 2 ADD_TEST(prTuMapProjectionsNew ${PROJECTIONS_TESTS} otbMapProjectionsNew ) +# TEST 3 ADD_TEST(prTuSensorModelsNew ${PROJECTIONS_TESTS} otbSensorModelsNew ) +# TEST 4 ADD_TEST(prTvTestCreateProjectionWithOSSIM_Toulouse ${PROJECTIONS_TESTS} otbCreateProjectionWithOSSIM ${IMAGEDATA}/TOULOUSE/QuickBird/000000128955_01_P001_PAN/02APR01105228-P1BS-000000128955_01_P001.TIF ) +# TEST 5 ADD_TEST(prTvTestCreateProjectionWithOSSIM_Cevennes ${PROJECTIONS_TESTS} otbCreateProjectionWithOSSIM ${IMAGEDATA}/CEVENNES/06FEB12104912-P1BS-005533998070_01_P001.TIF ) - +# TEST 6 ADD_TEST(prTvTestCreateProjectionWithOTB_Toulouse ${PROJECTIONS_TESTS} otbCreateProjectionWithOTB ${IMAGEDATA}/TOULOUSE/QuickBird/000000128955_01_P001_PAN/02APR01105228-P1BS-000000128955_01_P001.TIF ) +# TEST 7 ADD_TEST(prTvTestCreateProjectionWithOTB_Cevennes ${PROJECTIONS_TESTS} otbCreateProjectionWithOTB ${IMAGEDATA}/CEVENNES/06FEB12104912-P1BS-005533998070_01_P001.TIF ) +# TEST 8 ADD_TEST(prTvTestCreateInverseForwardSensorModel_Toulouse ${PROJECTIONS_TESTS} otbCreateInverseForwardSensorModel ${IMAGEDATA}/TOULOUSE/QuickBird/000000128955_01_P001_PAN/02APR01105228-P1BS-000000128955_01_P001.TIF ) + +# TEST 9 ADD_TEST(prTvTestCreateInverseForwardSensorModel_Cevennes ${PROJECTIONS_TESTS} otbCreateInverseForwardSensorModel ${IMAGEDATA}/CEVENNES/06FEB12104912-P1BS-005533998070_01_P001.TIF ) - -ADD_TEST(prTvRegionProjectionToulouse ${PROJECTIONS_TESTS} +# TEST 10 +ADD_TEST(prTvRegionProjectionToulouse ${PROJECTIONS_TESTS} + --compare-image ${TOL} ${BASELINE}/prTvRegionProjectionToulouse.tif + ${TEMP}/prTvRegionProjectionToulouse.tif otbRegionProjection ${IMAGEDATA}/TOULOUSE/QuickBird/000000128955_01_P001_PAN/02APR01105228-P1BS-000000128955_01_P001.TIF - ${TEMP}/RegionProjectionToulouse.png - 43.60 1.44 # Coordonn�es Lat/Long du Capitole Toulouse + ${TEMP}/prTvRegionProjectionToulouse.tif + 43.60 1.44 # Coordonnees Lat/Long du Capitole Toulouse 500 500 10 -# 0.1 -# 0.1 ) +# TEST 11 ADD_TEST(prTvRegionProjectionResamplerToulouse ${PROJECTIONS_TESTS} + --compare-image ${TOL} ${BASELINE}/prTvRegionProjectionResamplerToulouse.tif + ${TEMP}/prTvRegionProjectionResamplerToulouse.tif otbRegionProjectionResampler ${IMAGEDATA}/TOULOUSE/QuickBird/000000128955_01_P001_PAN/02APR01105228-P1BS-000000128955_01_P001.TIF - ${TEMP}/prTvRegionProjectionResamplerToulouse.hdr - 43.60 1.44 # Coordonn�es Lat/Long du Capitole Toulouse - 2000 - 2000 - 200 + ${TEMP}/prTvRegionProjectionResamplerToulouse.tif + 43.60 1.44 # Coordonnees Lat/Long du Capitole Toulouse + 500 + 500 + 100 0.00001 0.00001 ) -ADD_TEST(prTvRegionProjectionResamplerCevennes ${PROJECTIONS_TESTS} - otbRegionProjectionResampler +# TEST 12 +ADD_TEST(prTvRegionProjectionCevennes ${PROJECTIONS_TESTS} + --compare-image ${TOL} ${BASELINE}/prTvRegionProjectionCevennes.tif + ${TEMP}/prTvRegionProjectionCevennes.tif + otbRegionProjection ${IMAGEDATA}/CEVENNES/06FEB12104912-P1BS-005533998070_01_P001.TIF - ${TEMP}/prTvRegionProjectionResamplerCevennes.tif + ${TEMP}/prTvRegionProjectionCevennes.tif 44.08 3.7 500 500 10 - 0.00001 - 0.00001 ) -ADD_TEST(prTvRegionProjectionCevennes ${PROJECTIONS_TESTS} - otbRegionProjection +# TEST 13 +ADD_TEST(prTvRegionProjectionResamplerCevennes ${PROJECTIONS_TESTS} + --compare-image ${TOL} ${BASELINE}/prTvRegionProjectionResamplerCevennes.tif + ${TEMP}/prTvRegionProjectionResamplerCevennes.tif + otbRegionProjectionResampler ${IMAGEDATA}/CEVENNES/06FEB12104912-P1BS-005533998070_01_P001.TIF - ${TEMP}/prTvRegionProjectionCevennes.png + ${TEMP}/prTvRegionProjectionResamplerCevennes.tif 44.08 3.7 500 500 10 -# 0.00001 -# 0.00001 + 0.00001 + 0.00001 ) + -ADD_TEST(prTvOrthoRectificationCevennes ${PROJECTIONS_TESTS} +#ADD_TEST(prTvOrthoRectificationToulouse ${PROJECTIONS_TESTS} +# otbOrthoRectificationFilter +# ${IMAGEDATA}/TOULOUSE/QuickBird/000000128955_01_P001_PAN/02APR01105228-P1BS-000000128955_01_P001.TIF +# ${TEMP}/prTvOrthoRectificationToulouse.tif +# 43.60 +# 1.44 # Coordonnees Lat/Long du Capitole Toulouse +# 500 +# 500 +# 10 +# 0.00001 +# 0.00001 +# ) + +# TEST 14 +ADD_TEST(prTvSensorImageToCartoCevennes ${PROJECTIONS_TESTS} + --compare-image ${TOL} ${BASELINE}/prTvSensorImageToCartoCevennes_UTM.tif + ${TEMP}/prTvSensorImageToCartoCevennes_UTM.tif + otbSensorImageToCarto + ${IMAGEDATA}/CEVENNES/06FEB12104912-P1BS-005533998070_01_P001.TIF + ${TEMP}/prTvSensorImageToCartoCevennes_UTM.tif + 556046 + 4.881e+06 + 500 + 500 + 220 + ) + +# TEST 15 +ADD_TEST(prTvOrthoRectificationCevennes ${PROJECTIONS_TESTS} + --compare-image ${TOL} ${BASELINE}/prTvOrthoRectificationCevennes_UTM.tif + ${TEMP}/prTvOrthoRectificationCevennes_UTM.tif otbOrthoRectificationFilter ${IMAGEDATA}/CEVENNES/06FEB12104912-P1BS-005533998070_01_P001.TIF - ${TEMP}/prTvOrthoRectificationCevennes.tif - 44.08 - 3.7 + ${TEMP}/prTvOrthoRectificationCevennes_UTM.tif + 556046 + 4.881e+06 500 - 500 - 10 - 0.00001 - 0.00001 + 500 + 220 + 0.5 + 0.5 ) +# TEST 16 +ADD_TEST(prTvSensorImageDEMToCartoCevennes ${PROJECTIONS_TESTS} + --compare-image ${TOL} ${BASELINE}/prTvSensorImageDEMToCartoCevennes_UTM.tif + ${TEMP}/prTvSensorImageDEMToCartoCevennes_UTM.tif + otbSensorImageDEMToCarto + ${IMAGEDATA}/CEVENNES/06FEB12104912-P1BS-005533998070_01_P001.TIF + ${TEMP}/prTvSensorImageDEMToCartoCevennes_UTM.tif + 556046 + 4.881e+06 + 500 + 500 + 200 + ${INPUTDATA}/DEM/srtm_directory/ + ${TEMP}/SensorImageDEMToCartoCevennes.hdr + ) -# ADD_TEST(prTvSensorImageToCartoCevennes ${PROJECTIONS_TESTS} -# otbSensorImageToCarto -# ${IMAGEDATA}/CEVENNES/06FEB12104912-P1BS-005533998070_01_P001.TIF -# ${TEMP}/SensorImageToCartoCevennes_UTM.png -# 556046 -# 4.881e+06 -# 2000 -# 2000 -# 220 -# ) - - -# ADD_TEST(prTvSensorImageDEMToCartoCevennes ${PROJECTIONS_TESTS} -# otbSensorImageDEMToCarto -# ${IMAGEDATA}/CEVENNES/06FEB12104912-P1BS-005533998070_01_P001.TIF -# ${TEMP}/SensorImageDEMToCartoCevennes_UTM_DEM.png -# 556046 -# 4.881e+06 -# 2000 -# 2000 -# 200 -# ${INPUTDATA}/DEM/srtm_directory/ -# ${TEMP}/SensorImageDEMToCartoCevennes.hdr -# ) - - - +# TEST 17 +ADD_TEST(prTvOrthoRectificationCevennesWithDEM ${PROJECTIONS_TESTS} + --compare-image ${TOL} ${BASELINE}/prTvOrthoRectificationCevennesWithDEM_UTM.tif + ${TEMP}/prTvOrthoRectificationCevennesWithDEM_UTM.tif + otbOrthoRectificationFilterWithDEM + ${IMAGEDATA}/CEVENNES/06FEB12104912-P1BS-005533998070_01_P001.TIF + ${TEMP}/prTvOrthoRectificationCevennesWithDEM_UTM.tif + 556046 + 4.881e+06 + 500 + 500 + 200 + 0.7 + 0.7 + ${INPUTDATA}/DEM/srtm_directory/ + ${TEMP}/prTvOrthoRectificationCevennesWithDEM_UTM.tif + ) SET(Projections_SRCS otbProjectionBaseNew.cxx @@ -155,6 +206,7 @@ otbCreateInverseForwardSensorModel.cxx otbRegionProjection.cxx otbRegionProjectionResampler.cxx otbOrthoRectificationFilter.cxx +otbOrthoRectificationFilterWithDEM.cxx otbSensorImageToCarto.cxx otbSensorImageDEMToCarto.cxx ) diff --git a/Testing/Code/Projections/otbOrthoRectificationFilter.cxx b/Testing/Code/Projections/otbOrthoRectificationFilter.cxx index abfb02cba70a2ac546186e359866f43241632168..9a30f0cd04154d8303cfc0048cf11cae56b572ec 100644 --- a/Testing/Code/Projections/otbOrthoRectificationFilter.cxx +++ b/Testing/Code/Projections/otbOrthoRectificationFilter.cxx @@ -76,14 +76,14 @@ int otbOrthoRectificationFilter( int argc, char* argv[] ) typedef otb::ImageFileReader<ImageType> ReaderType; // typedef otb::ImageFileWriter<ImageType> WriterType; typedef otb::StreamingImageFileWriter<ImageType> WriterType; - + typedef itk::Transform<double,2,2> TransformType; // typedef otb::InverseSensorModel<double> ModelType; typedef itk::LinearInterpolateImageFunction<ImageType, double > InterpolatorType; typedef otb::StreamingResampleImageFilter<ImageType,ImageType> ResamplerType; - typedef otb::UtmProjection utmMapProjectionType ; - typedef otb::OrthoRectificationFilter<ImageType, ImageType, utmMapProjectionType> OrthoRectifFilterType ; + typedef otb::UtmProjection UtmMapProjectionType ; + typedef otb::OrthoRectificationFilter<ImageType, ImageType, UtmMapProjectionType> OrthoRectifFilterType ; typedef otb::InverseSensorModel<double> ModelType; //Allocate pointer @@ -95,6 +95,7 @@ int otbOrthoRectificationFilter( int argc, char* argv[] ) // ResamplerType::Pointer orthoRectifFilter=ResamplerType::New(); OrthoRectifFilterType::Pointer orthoRectifFilter=OrthoRectifFilterType::New(); + UtmMapProjectionType::Pointer utmMapProjection = UtmMapProjectionType::New(); // Set parameters ... reader->SetFileName(argv[1]); @@ -127,7 +128,12 @@ int otbOrthoRectificationFilter( int argc, char* argv[] ) origin[0]=strtod(argv[3], NULL); //latitude de l'origine. origin[1]=strtod(argv[4], NULL); //longitude de l'origine. orthoRectifFilter->SetOutputOrigin(origin); -// orthoRectifFilter->SetTransform( model ); + + utmMapProjection->SetZone(31); + utmMapProjection->SetHemisphere('N'); + orthoRectifFilter->SetMapProjection(utmMapProjection); + + // orthoRectifFilter->SetTransform( model ); // model->SetImageGeometry(reader->GetOutput()->GetImageKeywordlist()); // resampler->SetTransform(model); // resampler->SetInterpolator(interpolator); @@ -139,7 +145,7 @@ int otbOrthoRectificationFilter( int argc, char* argv[] ) writer->SetInput(orthoRectifFilter->GetOutput()); - writer->SetTilingStreamDivisions(20); + writer->SetTilingStreamDivisions(); // writer->SetNumberOfStreamDivisions(1000); otbGenericMsgDebugMacro(<< "Update writer ..." ); writer->Update(); diff --git a/Testing/Code/Projections/otbOrthoRectificationFilterWithDEM.cxx b/Testing/Code/Projections/otbOrthoRectificationFilterWithDEM.cxx new file mode 100644 index 0000000000000000000000000000000000000000..ee497bc22df1431ed1d9e76178db7ae3a2b03e2d --- /dev/null +++ b/Testing/Code/Projections/otbOrthoRectificationFilterWithDEM.cxx @@ -0,0 +1,164 @@ +/*========================================================================= + + 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. + +=========================================================================*/ +#if defined(_MSC_VER) +#pragma warning ( disable : 4786 ) +#endif + +/*! + * + * PURPOSE: + * + * Application pour projeter une r�gion d'une image en coordonn�es g�ographiques + * en utilisant un Interpolator+regionextractor et un Iterator. + * + */ + +// iostream is used for general output +#include <iostream> +#include <iterator> +#include <stdlib.h> + +#include "otbMacro.h" +#include "otbImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "otbStreamingImageFileWriter.h" +#include "otbInverseSensorModel.h" +#include "otbStreamingResampleImageFilter.h" + +#include "itkExceptionObject.h" +#include "itkExtractImageFilter.h" +#include "itkResampleImageFilter.h" +#include "itkRescaleIntensityImageFilter.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkLinearInterpolateImageFunction.h" + +#include "otbOrthoRectificationFilter.h" +#include "otbMapProjections.h" + +#include "init/ossimInit.h" + + + +int otbOrthoRectificationFilterWithDEM( int argc, char* argv[] ) +{ + try + { + + ossimInit::instance()->initialize(argc, argv); + + if(argc!=12) + { + std::cout << argv[0] <<" <input filename> <output filename> <latitude de l'origine> <longitude de l'origine> <taille_x> <taille_y>"; + std::cout << " <NumberOfstreamDivisions> <x_spacing> <y_spacing> <srtm directory> <DEM Image Filename>" << std::endl; + + return EXIT_FAILURE; + } + + + typedef otb::Image<unsigned int, 2> ImageType; + typedef otb::ImageFileReader<ImageType> ReaderType; +// typedef otb::ImageFileWriter<ImageType> WriterType; + typedef otb::StreamingImageFileWriter<ImageType> WriterType; + typedef itk::Transform<double,2,2> TransformType; + +// typedef otb::InverseSensorModel<double> ModelType; + typedef itk::LinearInterpolateImageFunction<ImageType, double > InterpolatorType; + + typedef otb::StreamingResampleImageFilter<ImageType,ImageType> ResamplerType; + typedef otb::UtmProjection UtmMapProjectionType ; + typedef otb::OrthoRectificationFilter<ImageType, ImageType, UtmMapProjectionType> OrthoRectifFilterType ; + typedef otb::InverseSensorModel<double> ModelType; + + //Allocate pointer + ReaderType::Pointer reader=ReaderType::New(); + WriterType::Pointer writer=WriterType::New(); +// ModelType::Pointer model=ModelType::New(); + InterpolatorType::Pointer interpolator=InterpolatorType::New(); + ModelType::Pointer model= ModelType::New(); + + OrthoRectifFilterType::Pointer orthoRectifFilter=OrthoRectifFilterType::New(); + UtmMapProjectionType::Pointer utmMapProjection = UtmMapProjectionType::New(); + + // Set parameters ... + reader->SetFileName(argv[1]); + writer->SetFileName(argv[2]); + + // Read meta data (ossimKeywordlist) + reader->GenerateOutputInformation(); + model->SetImageGeometry(reader->GetOutput()->GetImageKeywordlist()); + +// std::cout << "Model1: " << model << std::endl; + + orthoRectifFilter->SetInput(reader->GetOutput()); + + ImageType::IndexType start; + start[0]=0; + start[1]=0; + orthoRectifFilter->SetOutputStartIndex(start); + + ImageType::SizeType size; + size[0]=atoi(argv[5]); //Taille en X. + size[1]=atoi(argv[6]); //Taille en Y. + orthoRectifFilter->SetSize(size); + + ImageType::SpacingType spacing; + spacing[0]=atof(argv[8]); + spacing[1]=atof(argv[9]); + orthoRectifFilter->SetOutputSpacing(spacing); + + ImageType::PointType origin; + origin[0]=strtod(argv[3], NULL); //latitude de l'origine. + origin[1]=strtod(argv[4], NULL); //longitude de l'origine. + orthoRectifFilter->SetOutputOrigin(origin); + + utmMapProjection->SetZone(31); + utmMapProjection->SetHemisphere('N'); + orthoRectifFilter->SetMapProjection(utmMapProjection); + + std::string srtmDirectory(argv[10]); + orthoRectifFilter->SetDEMDirectory(srtmDirectory); + + writer->SetInput(orthoRectifFilter->GetOutput()); + + writer->SetTilingStreamDivisions(); + // writer->SetNumberOfStreamDivisions(1000); + otbGenericMsgDebugMacro(<< "Update writer ..." ); + writer->Update(); + + } + catch( itk::ExceptionObject & err ) + { + std::cout << "Exception itk::ExceptionObject levee !" << std::endl; + std::cout << err << std::endl; + return EXIT_FAILURE; + } + catch( std::bad_alloc & err ) + { + std::cout << "Exception bad_alloc : "<<(char*)err.what()<< std::endl; + return EXIT_FAILURE; + } + catch( ... ) + { + std::cout << "Exception levee inconnue !" << std::endl; + return EXIT_FAILURE; + } + return EXIT_SUCCESS; + + }//Fin main() + diff --git a/Testing/Code/Projections/otbProjectionsTests.cxx b/Testing/Code/Projections/otbProjectionsTests.cxx index edb451afa5b52fabae395037ad01c220e37b247a..8ead8d00bf190f94b0536b90a543f472c464330c 100755 --- a/Testing/Code/Projections/otbProjectionsTests.cxx +++ b/Testing/Code/Projections/otbProjectionsTests.cxx @@ -36,6 +36,7 @@ REGISTER_TEST(otbCreateInverseForwardSensorModel); REGISTER_TEST(otbRegionProjection); REGISTER_TEST(otbRegionProjectionResampler); REGISTER_TEST(otbOrthoRectificationFilter); +REGISTER_TEST(otbOrthoRectificationFilterWithDEM); REGISTER_TEST(otbSensorImageToCarto); REGISTER_TEST(otbSensorImageDEMToCarto); } diff --git a/Testing/Code/Projections/otbRegionProjection.cxx b/Testing/Code/Projections/otbRegionProjection.cxx index 0337cc05570fd415e2c24ad21894032167a2f7e2..a8227b3dc1da7e968edf50f736e1f9d2359ffb1e 100755 --- a/Testing/Code/Projections/otbRegionProjection.cxx +++ b/Testing/Code/Projections/otbRegionProjection.cxx @@ -157,6 +157,8 @@ otbGenericMsgDebugMacro(<< "ossimKeywordlist: "<<geom_kwl ); otbGenericMsgDebugMacro(<< "InverseSensorModel created " ); ModelType::OutputPointType inputpoint; + + otbGenericMsgDebugMacro(<< "Miary Sensor Model :" << model); //reader1->SetFileName(argv[2]); @@ -197,7 +199,8 @@ otbGenericMsgDebugMacro(<< "Interpolator created " ); /********************************************************/ /* Cr�ation de notre writer */ /********************************************************/ -typedef otb::ImageFileWriter<CharImageType> CharWriterType; +//typedef otb::ImageFileWriter<CharImageType> CharWriterType; +typedef otb::ImageFileWriter<ImageType> CharWriterType; typedef otb::ImageFileWriter<ImageType> WriterType; WriterType::Pointer extractorwriter=WriterType::New(); CharWriterType::Pointer writer=CharWriterType::New(); @@ -380,10 +383,13 @@ otbGenericMsgDebugMacro(<< "currentIndexArray deleted" ); writer->SetFileName(argv[2]); otbGenericMsgDebugMacro(<< "FilenameSet" ); -rescaler->SetInput(outputimage); +//rescaler->SetInput(outputimage); CharImageType::Pointer charoutputimage=CharImageType::New(); -charoutputimage=rescaler->GetOutput(); -writer->SetInput(charoutputimage); +//charoutputimage=rescaler->GetOutput(); +//writer->SetInput(charoutputimage); + +std::cout << "Image before Writer" << outputimage->GetOrigin() << std::endl; +writer->SetInput(outputimage); writer->Update(); otbGenericMsgDebugMacro(<< "Outputimage created" ); diff --git a/Testing/Code/Projections/otbRegionProjectionResampler.cxx b/Testing/Code/Projections/otbRegionProjectionResampler.cxx index d33f662cec974391b267a4443b397ea6222813a0..29d35ac479aa05a0959b5bcaece8e1474d3a6fcd 100755 --- a/Testing/Code/Projections/otbRegionProjectionResampler.cxx +++ b/Testing/Code/Projections/otbRegionProjectionResampler.cxx @@ -46,8 +46,10 @@ #include "itkRescaleIntensityImageFilter.h" #include "itkImageRegionIteratorWithIndex.h" #include "itkLinearInterpolateImageFunction.h" +#include "itkTranslationTransform.h" #include "otbInverseSensorModel.h" +#include "otbCompositeTransform.h" #include "init/ossimInit.h" @@ -70,7 +72,7 @@ int otbRegionProjectionResampler( int argc, char* argv[] ) typedef otb::Image<unsigned char, 2> CharImageType; - typedef otb::Image<unsigned short, 2> ImageType; + typedef otb::Image<unsigned int, 2> ImageType; typedef otb::ImageFileReader<ImageType> ReaderType; // typedef otb::ImageFileWriter<ImageType> WriterType; typedef otb::StreamingImageFileWriter<ImageType> WriterType; @@ -83,6 +85,10 @@ int otbRegionProjectionResampler( int argc, char* argv[] ) // ResamplerType; typedef otb::StreamingResampleImageFilter< ImageType, ImageType > ResamplerType; + typedef itk::TranslationTransform<double,2> TransformType; + + typedef otb::CompositeTransform<ModelType,TransformType> CompositeType; + ImageType::IndexType start; ImageType::SizeType size; ImageType::SpacingType spacing; @@ -96,7 +102,9 @@ int otbRegionProjectionResampler( int argc, char* argv[] ) InterpolatorType::Pointer interpolator=InterpolatorType::New(); RescalerType::Pointer rescaler=RescalerType::New(); ResamplerType::Pointer resampler = ResamplerType::New(); - +// TransformType::Pointer translationTransform = TransformType::New(); +// CompositeType::Pointer compositeTransform = CompositeType::New(); + // Set parameters ... reader->SetFileName(argv[1]); writer->SetFileName(argv[2]); @@ -110,6 +118,14 @@ int otbRegionProjectionResampler( int argc, char* argv[] ) // otb::ImageKeywordlist ImageKeyworlist = reader->GetOutput()->GetImageKeywordlist(); model->SetImageGeometry(reader->GetOutput()->GetImageKeywordlist()); +/* compositeTransform->SetFirstTransform(model); + + itk::Vector<double,2> offset; + offset[0]=atof(argv[10]); + offset[1]=atof(argv[11]); + translationTransform->SetOffset(offset); + + compositeTransform->SetSecondTransform(translationTransform);*/ start[0]=0; start[1]=0; @@ -131,6 +147,8 @@ int otbRegionProjectionResampler( int argc, char* argv[] ) // region = inputImage->GetBufferedRegion(); */ + std::cout << "Origin " << origin << std::endl; + resampler->SetOutputSpacing( spacing ); resampler->SetOutputOrigin( origin ); resampler->SetSize( region.GetSize() ); @@ -141,6 +159,8 @@ int otbRegionProjectionResampler( int argc, char* argv[] ) resampler->SetTransform( model ); resampler->SetInterpolator( interpolator ); + otbGenericMsgDebugMacro(<< "Romain Sensor Model :" << model); + // otbGenericMsgDebugMacro(<< "Resampler initialized !! " ); // resampler->Update(); @@ -148,10 +168,10 @@ int otbRegionProjectionResampler( int argc, char* argv[] ) writer->SetInput(resampler->GetOutput()); // writer->SetNumberOfStreamDivisions(1000); - writer->SetTilingStreamDivisions(); +// writer->SetTilingStreamDivisions(atoi(argv[7])); + writer->SetTilingStreamDivisions(10); otbGenericMsgDebugMacro(<< "Update writer ..." ); writer->Update(); - } catch( itk::ExceptionObject & err ) { diff --git a/Testing/Code/Projections/otbSensorImageDEMToCarto.cxx b/Testing/Code/Projections/otbSensorImageDEMToCarto.cxx index 871170494c02f1043172ca7ba343df4d8eda3f20..37373db5f7a5f80faf9ab17dc83ed1ae79d5e512 100755 --- a/Testing/Code/Projections/otbSensorImageDEMToCarto.cxx +++ b/Testing/Code/Projections/otbSensorImageDEMToCarto.cxx @@ -124,7 +124,7 @@ ossimKeywordlist geom_kwl; /* Cr�ation d'un DEMReader */ /********************************************************/ -typedef otb::DEMReader<DEMImageType> DEMReaderType; +typedef otb::DEMReader DEMReaderType; DEMReaderType::Pointer otbElevManager=DEMReaderType::New(); double height; @@ -261,7 +261,7 @@ outputimage->TransformIndexToPhysicalPoint(currentindex, outputpoint); << "Le point physique correspondant est: ("<< outputpoint[0]<< ","<< outputpoint[1]<< ")"); //On applique la projection: -geoPoint= utmprojection->Inverse(outputpoint); +geoPoint= utmprojection->TransformPoint(outputpoint); otbMsgDevMacro(<< "Le point g�ographique correspondant est: ("<< geoPoint[0]<< ","<< geoPoint[1]<< ")"); //on calcule son altitude @@ -308,7 +308,7 @@ min_y=pixelIndexArray[1]; if(j%2!=0 && pixelIndexArray[j]<min_y){min_y=pixelIndexArray[j];} }//Fin while - otbGenericMsgDebugMacro(<< "max_x=" << max_x<< std::endl + otbMsgDevMacro(<< "max_x=" << max_x<< std::endl << "max_y=" << max_y<< std::endl << "min_x=" << min_x<< std::endl << "min_y=" << min_y); @@ -365,11 +365,11 @@ DEMimage->SetPixel(currentindexbis,DEMpixelvalue); otbMsgDevMacro(<< "Altitude stock�e:" << heightArray[k] ); } delete pixelIndexArray; -otbGenericMsgDebugMacro(<< "pixelIndexArray deleted" ); +otbMsgDevMacro(<< "pixelIndexArray deleted" ); delete currentIndexArray; -otbGenericMsgDebugMacro(<< "currentIndexArray deleted" ); +otbMsgDevMacro(<< "currentIndexArray deleted" ); delete heightArray; -otbGenericMsgDebugMacro(<< "heightArray deleted" ); +otbMsgDevMacro(<< "heightArray deleted" ); }//Fin boucle principale //Cr�ation de l'image de sortie diff --git a/Testing/Code/Projections/otbSensorImageToCarto.cxx b/Testing/Code/Projections/otbSensorImageToCarto.cxx index b2447971e01ff88999d269606e6d79e5d2c8d8b5..c3d9506dd85e4c81c25d5759f9af9ff1ec2a9a9e 100755 --- a/Testing/Code/Projections/otbSensorImageToCarto.cxx +++ b/Testing/Code/Projections/otbSensorImageToCarto.cxx @@ -295,7 +295,7 @@ outputimage->TransformIndexToPhysicalPoint(currentindex, outputpoint); << "Le point physique correspondant est: ("<< outputpoint[0]<< ","<< outputpoint[1]<< ")"); //On applique la projection: -geoPoint= utmprojection->Inverse(outputpoint); +geoPoint= utmprojection->TransformPoint(outputpoint); otbMsgDevMacro(<< "Le point g�ographique correspondant est: ("<< geoPoint[0]<< ","<< geoPoint[1]<< ")"); //On calcule les coordonn�es pixeliques sur l'image capteur inputpoint = model->TransformPoint(geoPoint);