diff --git a/Code/BasicFilters/otbBinaryImageDensityFunction.h b/Code/BasicFilters/otbBinaryImageDensityFunction.h index bff63bbaa8e1596e3484e7361457e34d2e6b9961..9dd890d564dfaa550db0b128bb52b6bbabe766f8 100644 --- a/Code/BasicFilters/otbBinaryImageDensityFunction.h +++ b/Code/BasicFilters/otbBinaryImageDensityFunction.h @@ -45,38 +45,29 @@ template <class TInputImage, class TCoordRep = float > class ITK_EXPORT BinaryImageDensityFunction : public itk::ImageFunction< TInputImage, typename itk::NumericTraits<typename TInputImage::PixelType>::RealType,TCoordRep > { -public: + public: /** Standard class typedefs. */ typedef BinaryImageDensityFunction Self; typedef itk::ImageFunction<TInputImage,typename itk::NumericTraits<typename TInputImage::PixelType>::RealType, - TCoordRep > Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; + TCoordRep > Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; /** Run-time type information (and related methods). */ itkTypeMacro(BinaryImageDensityFunction, itk::ImageFunction); - + /** Method for creation through the object factory. */ itkNewMacro(Self); - - /** InputImageType typedef support. */ - typedef TInputImage InputImageType; - - /** OutputType typdef support. */ - typedef typename Superclass::OutputType OutputType; - - /** Index typedef support. */ - typedef typename Superclass::IndexType IndexType; - /** ContinuousIndex typedef support. */ + /** InputImageType typedef support. */ + typedef TInputImage InputImageType; + typedef typename InputImageType::SizeType RadiusType; + typedef typename Superclass::OutputType OutputType; + typedef typename Superclass::IndexType IndexType; typedef typename Superclass::ContinuousIndexType ContinuousIndexType; - - /** Point typedef support. */ - typedef typename Superclass::PointType PointType; - - /** Dimension of the underlying image. */ -itkStaticConstMacro(ImageDimension, unsigned int, - InputImageType::ImageDimension); + typedef typename Superclass::PointType PointType; + + itkStaticConstMacro(ImageDimension, unsigned int, InputImageType::ImageDimension); /** Datatype used for the density */ typedef typename itk::NumericTraits<typename InputImageType::PixelType>::RealType @@ -102,8 +93,13 @@ itkStaticConstMacro(ImageDimension, unsigned int, /** Get/Set the radius of the neighborhood over which the statistics are evaluated */ - itkSetMacro( NeighborhoodRadius, unsigned int ); - itkGetConstReferenceMacro( NeighborhoodRadius, unsigned int ); + itkSetMacro( NeighborhoodRadius, RadiusType ); + itkGetConstReferenceMacro( NeighborhoodRadius, RadiusType ); + void SetNeighborhoodRadius(unsigned int rad) + { + m_NeighborhoodRadius.Fill( rad ); + this->Modified(); + } protected: @@ -115,7 +111,7 @@ private: BinaryImageDensityFunction( const Self& ); //purposely not implemented void operator=( const Self& ); //purposely not implemented - unsigned int m_NeighborhoodRadius; + RadiusType m_NeighborhoodRadius; }; diff --git a/Code/BasicFilters/otbBinaryImageDensityFunction.txx b/Code/BasicFilters/otbBinaryImageDensityFunction.txx index faabc4abc3d91a0402bb6f024f3e2ede492051ae..6a992cc5d308d6f83609bce024a6a2036da0c0fe 100644 --- a/Code/BasicFilters/otbBinaryImageDensityFunction.txx +++ b/Code/BasicFilters/otbBinaryImageDensityFunction.txx @@ -35,7 +35,7 @@ template <class TInputImage, class TCoordRep> BinaryImageDensityFunction<TInputImage,TCoordRep> ::BinaryImageDensityFunction() { - m_NeighborhoodRadius = 1; + m_NeighborhoodRadius.Fill( 1 ); } @@ -77,11 +77,10 @@ BinaryImageDensityFunction<TInputImage,TCoordRep> } // Create an N-d neighborhood kernel, using a zeroflux boundary condition - typename InputImageType::SizeType kernelSize; - kernelSize.Fill( m_NeighborhoodRadius ); - + typename InputImageType::SizeType kernelSize = m_NeighborhoodRadius; + itk::ConstNeighborhoodIterator<InputImageType> - it(kernelSize, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion()); + it(kernelSize, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion()); // Set the iterator at the desired location it.SetLocation(index); diff --git a/Code/BasicFilters/otbBinaryImageToDensityImageFilter.h b/Code/BasicFilters/otbBinaryImageToDensityImageFilter.h index b8abde15112aef7f3cb19d7c0e36a24e94e9ac2c..6e9564f4cefb5018639806d3d83508693061ea4b 100644 --- a/Code/BasicFilters/otbBinaryImageToDensityImageFilter.h +++ b/Code/BasicFilters/otbBinaryImageToDensityImageFilter.h @@ -20,6 +20,8 @@ #include "itkImageToImageFilter.h" #include "itkDataObject.h" +#include "itkConstNeighborhoodIterator.h" + namespace otb { @@ -48,6 +50,7 @@ public: typedef TInputImage InputImageType; typedef typename InputImageType::RegionType InputImageRegionType; typedef typename InputImageType::Pointer InputImagePointerType; + typedef typename InputImageType::SizeType InputImageSizeType; typedef TOutputImage OutputImageType; typedef typename OutputImageType::Pointer OutputImagePointerType; @@ -55,16 +58,24 @@ public: typedef TCountFunction CountFunctionType; typedef typename CountFunctionType::Pointer CountFunctionPointerType; + typedef itk::ConstNeighborhoodIterator<TInputImage> NeighborhoodIteratorType; + typedef typename NeighborhoodIteratorType::RadiusType RadiusType; + /** Shrink factor accessor */ - itkSetMacro(NeighborhoodRadius,unsigned int); - itkGetMacro(NeighborhoodRadius, unsigned int); + itkSetMacro(NeighborhoodRadius, RadiusType); + itkGetMacro(NeighborhoodRadius, RadiusType); + void SetNeighborhoodRadius(unsigned int rad) + { + m_NeighborhoodRadius.Fill(rad); + this->Modified(); + } /** Main computation method */ - virtual void ThreadedGenerateData( const InputImageRegionType &outputRegionForThread, int threadId ) ; - + virtual void BeforeThreadedGenerateData(); + virtual void GenerateInputRequestedRegion(); protected: /** Constructor */ @@ -81,7 +92,7 @@ private: CountFunctionPointerType m_CountFunction; /** The shrink factor */ - unsigned int m_NeighborhoodRadius; + RadiusType m_NeighborhoodRadius; }; } // End namespace otb diff --git a/Code/BasicFilters/otbBinaryImageToDensityImageFilter.txx b/Code/BasicFilters/otbBinaryImageToDensityImageFilter.txx index 2a95dd2579a77402a9dc662c3eb61f11f77598cf..75691fe081a8faafe0b2ce9d0d2f767068868d0e 100644 --- a/Code/BasicFilters/otbBinaryImageToDensityImageFilter.txx +++ b/Code/BasicFilters/otbBinaryImageToDensityImageFilter.txx @@ -20,7 +20,12 @@ #include "otbBinaryImageToDensityImageFilter.h" #include "itkImageRegionIterator.h" -#include "itkImageRegionConstIterator.h" +//#include "itkImageRegionConstIterator.h" +#include "itkProgressReporter.h" +#include "otbMirrorBoundaryCondition.h" +#include "itkZeroFluxNeumannBoundaryCondition.h" +#include "itkNeighborhoodAlgorithm.h" + #include "otbMacro.h" namespace otb @@ -30,11 +35,10 @@ template <class TInputImage, class TOutputImage, class TCountFunction> BinaryImageToDensityImageFilter<TInputImage, TOutputImage, TCountFunction> ::BinaryImageToDensityImageFilter() { - m_NeighborhoodRadius = 1; - + m_NeighborhoodRadius.Fill( 1 ); m_CountFunction = CountFunctionType::New(); - } + /** Destructor */ template <class TInputImage, class TOutputImage, class TCountFunction> BinaryImageToDensityImageFilter<TInputImage, TOutputImage, TCountFunction> @@ -42,39 +46,119 @@ BinaryImageToDensityImageFilter<TInputImage, TOutputImage, TCountFunction> {} +template <class TInputImage, class TOutputImage, class TCountFunction> +void +BinaryImageToDensityImageFilter<TInputImage, TOutputImage, TCountFunction> +::GenerateInputRequestedRegion() +{ + // call the superclass' implementation of this method + Superclass::GenerateInputRequestedRegion(); + + // get pointers to the input and output + InputImagePointerType inputPtr = const_cast< TInputImage * >( this->GetInput()); + OutputImagePointerType outputPtr = this->GetOutput(); + + if ( !inputPtr || !outputPtr ) + { + return; + } + // get a copy of the input requested region (should equal the output + // requested region) + InputImageRegionType inputRequestedRegion = inputPtr->GetRequestedRegion(); + + // pad the input requested region by the operator radius + inputRequestedRegion.PadByRadius( m_NeighborhoodRadius ); + + // crop the input requested region at the input's largest possible region + if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) ) + { + inputPtr->SetRequestedRegion( inputRequestedRegion ); + return; + } + else + { + // Couldn't crop the region (requested region is outside the largest + // possible region). Throw an exception. + + // store what we tried to request (prior to trying to crop) + inputPtr->SetRequestedRegion( inputRequestedRegion ); + + // build an exception + itk::InvalidRequestedRegionError e(__FILE__, __LINE__); + itk::OStringStream msg; + msg << this->GetNameOfClass() + << "::GenerateInputRequestedRegion()"; + e.SetLocation(msg.str().c_str()); + e.SetDescription("Requested region is (at least partially) outside the largest possible region."); + e.SetDataObject(inputPtr); + throw e; + } +} + +template <class TInputImage, class TOutputImage, class TCountFunction> +void +BinaryImageToDensityImageFilter<TInputImage, TOutputImage, TCountFunction> +::BeforeThreadedGenerateData() +{ + Superclass::BeforeThreadedGenerateData(); + + m_CountFunction->SetInputImage(this->GetInput()); + m_CountFunction->SetNeighborhoodRadius(m_NeighborhoodRadius); +} + + /** Main computation method */ template <class TInputImage, class TOutputImage, class TCountFunction> void BinaryImageToDensityImageFilter<TInputImage, TOutputImage, TCountFunction> ::ThreadedGenerateData( const InputImageRegionType &outputRegionForThread, int threadId ) { - - InputImagePointerType input = const_cast<InputImageType * > (this->GetInput()); - OutputImagePointerType output = this->GetOutput(); + InputImagePointerType inputPtr = const_cast<InputImageType * > (this->GetInput()); + OutputImagePointerType outputPtr = this->GetOutput(); - m_CountFunction->SetInputImage(input); + itk::ZeroFluxNeumannBoundaryCondition<TInputImage> nbc; + RadiusType r; + r[0]=m_NeighborhoodRadius[0]; + r[1]=m_NeighborhoodRadius[1]; - itk::ImageRegionConstIterator<InputImageType> it(input,outputRegionForThread ); - itk::ImageRegionIterator<OutputImageType> itOut(output,outputRegionForThread ); - - it.GoToBegin(); - itOut.GoToBegin(); - - while(!it.IsAtEnd()) - { - m_CountFunction->SetNeighborhoodRadius(m_NeighborhoodRadius); - typename InputImageType::IndexType index = it.GetIndex(); + NeighborhoodIteratorType it; + typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<TInputImage>::FaceListType faceList; + typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<TInputImage> bC; + faceList = bC(inputPtr, outputRegionForThread, r); + typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<TInputImage>::FaceListType::iterator fit; + + itk::ImageRegionIterator<OutputImageType> itOut(outputPtr,outputRegionForThread ); + + itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); + + + typename InputImageType::IndexType index; + + for (fit=faceList.begin(); fit != faceList.end(); ++fit) + { + it = itk::ConstNeighborhoodIterator<TInputImage>(r, inputPtr, *fit); + + itOut = itk::ImageRegionIterator<TOutputImage>(outputPtr, *fit); + it.OverrideBoundaryCondition(&nbc); + it.GoToBegin(); - if(outputRegionForThread.IsInside(index)) + while(!itOut.IsAtEnd()) { - itOut.Set(m_CountFunction->EvaluateAtIndex(index)); + index = it.GetIndex(); + + if(outputRegionForThread.IsInside(index)) + { + itOut.Set(m_CountFunction->EvaluateAtIndex(index)); + } + + ++itOut; + ++it; + progress.CompletedPixel(); // potential exception thrown here } - - ++itOut; - ++it; } } + /** PrintSelf method */ template <class TInputImage, class TOutputImage, class TCountFunction> void @@ -82,8 +166,7 @@ BinaryImageToDensityImageFilter<TInputImage, TOutputImage, TCountFunction> ::PrintSelf(std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf(os,indent); - os << indent << "Neighborhood Radius : " << m_NeighborhoodRadius - << std::endl; + os << indent << "Neighborhood Radius : " << m_NeighborhoodRadius << std::endl; } } // End namespace otb #endif diff --git a/Code/BasicFilters/otbEdgeDensityImageFilter.h b/Code/BasicFilters/otbEdgeDensityImageFilter.h index d56fa7945a2118c50b0de6557babaab5884f0f25..0f8e34e9c0e8602b833d1797cb935d9ace0c6f68 100644 --- a/Code/BasicFilters/otbEdgeDensityImageFilter.h +++ b/Code/BasicFilters/otbEdgeDensityImageFilter.h @@ -85,10 +85,11 @@ public: itkSetMacro( NeighborhoodRadius, unsigned int ); itkGetConstReferenceMacro( NeighborhoodRadius, unsigned int ); + /**Set/Get detector */ + itkSetObjectMacro(Detector, DetectorType); + itkGetObjectMacro(Detector, DetectorType); + itkGetObjectMacro(DensityImageFilter, DensityImageType); - /**Set/Get Descriptor from the otbCountmageFunction*/ - virtual void SetDetector(DetectorType* detector); - virtual DetectorType* GetDetector(); protected: @@ -107,7 +108,6 @@ protected: /** * Main computation method. */ - //virtual void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId ); virtual void GenerateData( ); private: @@ -115,9 +115,9 @@ private: EdgeDensityImageFilter(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented - DetectorPointerType m_Detector; - DensityImagePointerType m_DensityImageFilter; - unsigned int m_NeighborhoodRadius; + DetectorPointerType m_Detector; + DensityImagePointerType m_DensityImageFilter; + unsigned int m_NeighborhoodRadius; }; } #ifndef OTB_MANUAL_INSTANTIATION diff --git a/Code/BasicFilters/otbEdgeDensityImageFilter.txx b/Code/BasicFilters/otbEdgeDensityImageFilter.txx index bef73b4e3410cfac91448832e4cb695eeccef976..e32dcb85261c664d2c971a4d5b51c20aa88fbe0f 100644 --- a/Code/BasicFilters/otbEdgeDensityImageFilter.txx +++ b/Code/BasicFilters/otbEdgeDensityImageFilter.txx @@ -29,7 +29,7 @@ namespace otb template <class TInputImage , class TOutputImage, class TEdgeDetector, class TDensityCount> EdgeDensityImageFilter<TInputImage, TOutputImage, TEdgeDetector, TDensityCount> ::EdgeDensityImageFilter() -{ +{ this->SetNumberOfRequiredInputs( 1 ); m_NeighborhoodRadius = 1; @@ -46,65 +46,35 @@ EdgeDensityImageFilter<TInputImage, TOutputImage, TEdgeDetector, TDensityCount> ::~EdgeDensityImageFilter() {} + /** * threaded Generate Data */ - -/** -* ThreadedGenerateData Performs the pixel-wise addition -*/ template <class TInputImage , class TOutputImage, class TEdgeDetector, class TDensityCount> void EdgeDensityImageFilter<TInputImage, TOutputImage, TEdgeDetector, TDensityCount> ::GenerateData() - //::ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int threadId ) { - typename Superclass::OutputImagePointer outputImage = this->GetOutput(); - InputImagePointerType ptr = const_cast<InputImageType *>(this->GetInput()); - if (!ptr) - return ; + m_Detector->SetInput( this->GetInput() ); - /** Apply Canny Detector*/ - m_Detector->SetInput(ptr); - + std::cout<<"###"<<this->GetInput()->GetLargestPossibleRegion()<<std::endl; + m_Detector->Update(); + m_Detector->UpdateOutputInformation(); + //m_Detector->SetRequestedRegionToLargestPossibleRegion(); + std::cout<<"~~~"<<m_Detector->GetOutput()->GetLargestPossibleRegion()<<std::endl; - /** Compute density on the binaruzed Image */ - m_DensityImageFilter->SetInput(m_Detector->GetOutput()); m_DensityImageFilter->SetNeighborhoodRadius(m_NeighborhoodRadius); + m_DensityImageFilter->SetInput(m_Detector->GetOutput()); + + m_DensityImageFilter->UpdateOutputInformation(); + std::cout<<"***"<<m_DensityImageFilter->GetOutput()->GetLargestPossibleRegion()<<std::endl; - /** updating the output*/ m_DensityImageFilter->GraftOutput(this->GetOutput()); m_DensityImageFilter->Update(); this->GraftOutput(m_DensityImageFilter->GetOutput()); } -/** - * Set Detector - */ -template <class TInputImage , class TOutputImage, class TEdgeDetector, class TDensityCount> -void -EdgeDensityImageFilter<TInputImage, TOutputImage, TEdgeDetector, TDensityCount> -::SetDetector(DetectorType* detector) -{ - m_Detector = detector; -} - - -/** - * Get Detector - */ -template <class TInputImage , class TOutputImage, class TEdgeDetector, class TDensityCount> -typename EdgeDensityImageFilter<TInputImage, TOutputImage, TEdgeDetector, TDensityCount> -::DetectorType * -EdgeDensityImageFilter<TInputImage, TOutputImage, TEdgeDetector, TDensityCount> -::GetDetector() -{ - return m_Detector; -} - - - /*---------------------------------------------------------------- PrintSelf -----------------------------------------------------------------*/ diff --git a/Code/Common/otbRectangle.h b/Code/Common/otbRectangle.h new file mode 100644 index 0000000000000000000000000000000000000000..fa4a161b140d6da58222eb73335b92e83853a35c --- /dev/null +++ b/Code/Common/otbRectangle.h @@ -0,0 +1,128 @@ +/*========================================================================= + +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 __otbRectangle_h +#define __otbRectangle_h + +#include "itkObject.h" +#include "itkImageRegion.h" + +#include "itkContinuousIndex.h" +#include "itkVectorContainer.h" + + + +namespace otb +{ + /** \class Rectangle + * \brief This class represent a Rectangle. + * + * A rectangle is defined by the median of the width, an orientation and a width + * + */ + template<class TValue=double> + class ITK_EXPORT Rectangle + : public itk::Object + { + public: + /** Standard typedefs */ + typedef Rectangle Self; + typedef itk::Object Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + typedef TValue ValueType; + + /** Type macro */ + itkNewMacro(Self); + + /** Creation through object factory macro */ + itkTypeMacro(Rectangle, itk::Object); + + /** Derived typedefs */ + typedef itk::ContinuousIndex<ValueType,2> ContinuousIndexType ; + typedef ContinuousIndexType VertexType; + typedef itk::VectorContainer<unsigned, VertexType> VertexListType; + typedef typename VertexListType::Pointer VertexListPointerType; + typedef typename VertexListType::ConstIterator VertexListConstIteratorType; + + + /* typedef typename Superclass::VertexType VertexType; */ + /* typedef typename Superclass::VertexListType VertexListType; */ + /* typedef typename VertexListType::Pointer VertexListPointerType; */ + /* typedef typename Superclass::ContinuousIndexType ContinuousIndexType; */ + /* typedef typename Superclass::VertexListConstIteratorType VertexListConstIteratorType; */ + + typedef itk::ImageRegion<2> RegionType; + + /** + * Set/Get the parameters the a rectangle + * - width and orientation + */ + + itkGetMacro(Width,ValueType); + itkSetMacro(Width,ValueType); + + itkSetMacro(Orientation,ValueType); + itkGetMacro(Orientation,ValueType); + + /** + * Check wether point is strictly inside the rectangle. + * \param point The point to check. + * \return True if the point is inside the polygon. + */ + bool IsInside(VertexType point) const; + + virtual void AddVertex (const ContinuousIndexType &vertex); + + /** GetBounding region*/ + virtual RegionType GetBoundingRegion() const; + + protected: + /** Constructor */ + Rectangle() + { + m_VertexList = VertexListType::New(); + }; + + /** Destructor */ + virtual ~Rectangle() {}; + + /**PrintSelf method */ + virtual void PrintSelf(std::ostream& os, itk::Indent indent) const; + + /** */ + virtual double ComputeEuclideanDistanceToSegment(VertexType q1, VertexType q2, VertexType p) const; + + /** */ + + + private: + Rectangle(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + ValueType m_Orientation; + ValueType m_Width; + + VertexListPointerType m_VertexList; + }; +}// End namespace otb + +#ifndef OTB_MANUAL_INSTANTIATION +#include "otbRectangle.txx" +#endif +#endif diff --git a/Code/Common/otbRectangle.txx b/Code/Common/otbRectangle.txx new file mode 100644 index 0000000000000000000000000000000000000000..a938b1567008bc8c30e9e2693e5ca2914c53676c --- /dev/null +++ b/Code/Common/otbRectangle.txx @@ -0,0 +1,223 @@ +/*========================================================================= + +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 __otbRectangle_txx +#define __otbRectangle_txx + +#include "otbRectangle.h" + +namespace otb +{ + + +/** + * Check wether point is strictly inside the polygon. + * \param point The point to check. + * \return True if the point is inside the polygon. + */ + +template<class TValue> +bool +Rectangle<TValue> +::IsInside(VertexType point) const +{ + /*Iterate through the inputVertexList*/ + if(m_VertexList->Size()<2) + itkGenericExceptionMacro(<<"Rectangle needs TWO vertex, up-to-date the start and the end of the segments with AdDVertex Method "); + + VertexListConstIteratorType it = m_VertexList->Begin(); + + VertexType p1 = it.Value(); + ++it; + VertexType p2 = it.Value(); + + + /** Compute Length of the rectangle*/ + double lengthSeg = vcl_sqrt( (p1[0]-p2[0])* (p1[0]-p2[0]) + (p1[1]-p2[1])* (p1[1]-p2[1]) ); + + /** Orthogonal segment containing the middle of the segment */ + VertexType middleP; + middleP[0] = (p1[0] + p2[0])/2.; + middleP[1] = (p1[1] + p2[1])/2.; + + + VertexType corner; + corner[0] = middleP[0] + (m_Width/2)*vcl_sin(m_Orientation); + corner[1] = middleP[1] - (m_Width/2)*vcl_cos(m_Orientation); + + + /** Compute the distance to the orthogonal median of the rectangle*/ + if( this->ComputeEuclideanDistanceToSegment( p1, p2 , point) - (m_Width/2.) < 1e-10 + && this->ComputeEuclideanDistanceToSegment(middleP , corner, point) - (lengthSeg/2.) < 1e-10) + return true; + else + return false; +} + + + /** + * Method to compute the distance between a point and a segment + */ + template<class TValue> + double + Rectangle<TValue> + ::ComputeEuclideanDistanceToSegment(VertexType q1, VertexType q2, VertexType p) const + { + double Xq1 = q1[0]; + double Yq1 = q1[1]; + double Xq2 = q2[0]; + double Yq2 = q2[1]; + double xp = p[0]; + double yp = p[1]; + + double SegmentLength = vcl_sqrt((Xq1-Xq2)* (Xq1-Xq2) + (Yq1-Yq2) *(Yq1-Yq2)); + double CrossProduct = Xq1*Yq2 - Xq2*Yq1 ; + double Num = vcl_abs(xp*(Yq1-Yq2) + yp*(Xq2-Xq1) + CrossProduct); + + /** distance from Point P to Segment Q1Q2*/ + return (Num/SegmentLength); + } + + /** + * Method to compute the Bounding region of a rectangle + */ + + template<class TValue> + typename Rectangle<TValue>:: + RegionType + Rectangle<TValue> + ::GetBoundingRegion() const + { + VertexListConstIteratorType it = m_VertexList->Begin(); + + VertexType p1 = it.Value(); + ++it; + VertexType p2 = it.Value(); + + /** Compute the four corners of the recatangle*/ + double dx = vcl_cos(m_Orientation); + double dy = vcl_sin(m_Orientation); + double halfWidth = m_Width/2; + + + VertexListPointerType cornersVertex = VertexListType::New(); + VertexType tempCorner; + + tempCorner[0] = p1[0] + dy* halfWidth ; + tempCorner[1] = p1[1] - dx* halfWidth ; + cornersVertex->InsertElement(cornersVertex->Size(),tempCorner); + + tempCorner[0] = p1[0] - dy* halfWidth ; + tempCorner[1] = p1[1] + dx* halfWidth ; + cornersVertex->InsertElement(cornersVertex->Size(),tempCorner); + + tempCorner[0] = p2[0] + dy* halfWidth ; + tempCorner[1] = p2[1] - dx* halfWidth ; + cornersVertex->InsertElement(cornersVertex->Size(),tempCorner); + + tempCorner[0] = p2[0] - dy* halfWidth ; + tempCorner[1] = p2[1] + dx* halfWidth ; + cornersVertex->InsertElement(cornersVertex->Size(),tempCorner); + + /** Compute the bounding region*/ + RegionType region; + typename RegionType::SizeType size; + typename RegionType::IndexType index; + typename RegionType::IndexType maxId; + + size.Fill(0); + index.Fill(0); + maxId.Fill(0); + + VertexListConstIteratorType itCorners = cornersVertex->Begin(); + + long int x = 0,y = 0; + + if ( cornersVertex->Size()>0) + { + x = static_cast<long int>(itCorners.Value()[0]); + y = static_cast<long int>(itCorners.Value()[1]); + index[0] = x; + index[1] = y; + + ++itCorners; + while (itCorners != cornersVertex->End()) + { + x = static_cast<long int>(itCorners.Value()[0]); + y = static_cast<long int>(itCorners.Value()[1]); + + // Index search + if ( x < index[0] ) + { + index[0] = x; + } + if ( y < index[1] ) + { + index[1] = y; + } + // Max Id search for size computation + if ( x > maxId[0] ) + { + maxId[0] = x; + } + if ( y > maxId[1] ) + { + maxId[1] = y; + } + + ++itCorners; + } + + size[0] = maxId[0] - index[0]; + size[1] = maxId[1] - index[1]; + } + region.SetSize(size); + region.SetIndex(index); + return region; + + } + + /* + * Method add vertex + * + */ +template<class TValue> +void +Rectangle<TValue> +::AddVertex(const ContinuousIndexType &vertex) +{ + if(m_VertexList->Size() > 1 ) + itkGenericExceptionMacro(<<"Rectangle needs only TWO vertex, a width and an orientation "); + + m_VertexList->InsertElement(m_VertexList->Size(), vertex); +} + +/** + * PrintSelf Method + */ +template<class TValue> +void +Rectangle<TValue> +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + os << indent << "Verticies: " << m_VertexList << std::endl; +} + + +} // End namespace otb + +#endif diff --git a/Code/Common/otbVectorDataExtractROI.txx b/Code/Common/otbVectorDataExtractROI.txx index 703ba6f978b22a555ab833fa5b7e20ce031d2bb8..33ab70b64830cd8b91b9a1ac90f1a0926c75671c 100644 --- a/Code/Common/otbVectorDataExtractROI.txx +++ b/Code/Common/otbVectorDataExtractROI.txx @@ -23,6 +23,7 @@ #include "otbGenericMapProjection.h" #include "itkIdentityTransform.h" +#include "otbGenericRSTransform.h" #include "otbObjectList.h" #include "otbMacro.h" @@ -297,13 +298,20 @@ void VectorDataExtractROI<TVectorData> ::ProjectRegionToInputVectorProjection() { - typedef otb::GenericMapProjection<otb::INVERSE> ForwardMapProjectionType; - ForwardMapProjectionType::Pointer mapTransform = ForwardMapProjectionType::New(); - mapTransform->SetWkt(m_ROI.GetRegionProjection()); - - /*FORWARD ; From RegionProjection To long/lat Projection */ - typename VertexListType::Pointer regionCorners = VertexListType::New(); - ProjPointType point1, point2 , point3, point4; + /* Use the RS Generic projection */ + typedef otb::GenericRSTransform<> GenericRSTransformType; + typename GenericRSTransformType::Pointer genericTransform = GenericRSTransformType::New(); + + /** We need get the transformation*/ + const itk::MetaDataDictionary &inputDict = this->GetInput()->GetMetaDataDictionary(); + genericTransform->SetInputDictionary(inputDict); + genericTransform->SetInputProjectionRef( m_ROI.GetRegionProjection()); + genericTransform->SetOutputProjectionRef(this->GetInput()->GetProjectionRef() ); + genericTransform->InstanciateTransform(); + + + typename VertexListType::Pointer regionCorners = VertexListType::New(); + ProjPointType point1, point2 , point3, point4; /** Compute the extremities of the region*/ point1[0] = m_ROI.GetOrigin()[0]; @@ -317,44 +325,12 @@ VectorDataExtractROI<TVectorData> point4[0] = m_ROI.GetOrigin()[0]; point4[1] = m_ROI.GetOrigin()[1]+ m_ROI.GetSize()[1]; - - /** Aplly the transformation*/ - ProjPointType pGeo1 = mapTransform->TransformPoint(point1); - ProjPointType pGeo2 = mapTransform->TransformPoint(point2); - ProjPointType pGeo3 = mapTransform->TransformPoint(point3); - ProjPointType pGeo4 = mapTransform->TransformPoint(point4); - - /** INVERSE : From long/lat to InputVectorData projection*/ - typedef otb::GenericMapProjection<otb::FORWARD> InverseMapProjectionType; - InverseMapProjectionType::Pointer mapInverseTransform = InverseMapProjectionType::New(); - - typedef itk::Transform<double, 2, 2> GenericTransformType; - GenericTransformType::Pointer outputTransform ; - - if(!this->GetInput()->GetProjectionRef().empty()) - { - mapInverseTransform->SetWkt(this->GetInput()->GetProjectionRef()); - outputTransform = mapInverseTransform.GetPointer(); - - //Case of kml - if(mapInverseTransform->GetMapProjection() == NULL) - outputTransform = itk::IdentityTransform< double, 2 >::New(); - } - else - outputTransform = itk::IdentityTransform< double, 2 >::New(); - - - /** Finally Project the four corners in the InputDataVector Projection*/ - ProjPointType pCarto1 = /*mapInverseTransform*/ outputTransform->TransformPoint(pGeo1); - ProjPointType pCarto2 = /*mapInverseTransform*/ outputTransform->TransformPoint(pGeo2); - ProjPointType pCarto3 = /*mapInverseTransform*/ outputTransform ->TransformPoint(pGeo3); - ProjPointType pCarto4 = /*mapInverseTransform*/ outputTransform->TransformPoint(pGeo4); /** Fill the vertex List : First Convert Point To*/ - regionCorners->InsertElement(regionCorners->Size(),this->PointToContinuousIndex(pCarto1)); - regionCorners->InsertElement(regionCorners->Size(),this->PointToContinuousIndex(pCarto2)); - regionCorners->InsertElement(regionCorners->Size(),this->PointToContinuousIndex(pCarto3)); - regionCorners->InsertElement(regionCorners->Size(),this->PointToContinuousIndex(pCarto4)); + regionCorners->InsertElement(regionCorners->Size(),this->PointToContinuousIndex(genericTransform->TransformPoint(point1))); + regionCorners->InsertElement(regionCorners->Size(),this->PointToContinuousIndex(genericTransform->TransformPoint(point2))); + regionCorners->InsertElement(regionCorners->Size(),this->PointToContinuousIndex(genericTransform->TransformPoint(point3))); + regionCorners->InsertElement(regionCorners->Size(),this->PointToContinuousIndex(genericTransform->TransformPoint(point4))); /** Due to The projection : the Projected ROI can be rotated */ m_GeoROI = this->ComputeVertexListBoudingRegion(regionCorners.GetPointer()); diff --git a/Testing/Code/Common/CMakeLists.txt b/Testing/Code/Common/CMakeLists.txt index 7170f7c439d2956f659d76a4d5828355f0292537..e033f6c331210ff22eb18f880a2f4d68116838b6 100644 --- a/Testing/Code/Common/CMakeLists.txt +++ b/Testing/Code/Common/CMakeLists.txt @@ -503,7 +503,20 @@ ADD_TEST(coTvUnaryFunctorNeighborhoodImageFilter ${COMMON_TESTS6} 3 ) +#--------------------otb::Rectangle +ADD_TEST(coTuRectangleNew ${COMMON_TESTS6} + otbRectangleNew) +ADD_TEST(coTvRectangle ${COMMON_TESTS6} +--compare-ascii ${EPS} + ${BASELINE_FILES}/otbRectangleTest.txt + ${TEMP}/otbRectangleTest.txt + otbRectangle + ${TEMP}/otbRectangleTest.txt + 2.12 2.35 12.54 2.45 3. 0.15 1. 1. + + +) # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ otbCommonTests7 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -663,6 +676,8 @@ otbPolygonNew.cxx otbPolygon.cxx otbUnaryFunctorNeighborhoodImageFilterNew.cxx otbUnaryFunctorNeighborhoodImageFilter.cxx +otbRectangleNew.cxx +otbRectangle.cxx ) SET(BasicCommon_SRCS7 diff --git a/Testing/Code/Common/otbCommonTests6.cxx b/Testing/Code/Common/otbCommonTests6.cxx index 41ebe6c6c90f806ca7acbaf2ab669d040667c98d..cf321e7878af74e3ed8e9e7e5d8731f03087adab 100644 --- a/Testing/Code/Common/otbCommonTests6.cxx +++ b/Testing/Code/Common/otbCommonTests6.cxx @@ -37,4 +37,6 @@ void RegisterTests() REGISTER_TEST(otbPolygon); REGISTER_TEST(otbUnaryFunctorNeighborhoodImageFilterNew); REGISTER_TEST(otbUnaryFunctorNeighborhoodImageFilter); + REGISTER_TEST(otbRectangleNew); + REGISTER_TEST(otbRectangle); } diff --git a/Testing/Code/Common/otbRectangle.cxx b/Testing/Code/Common/otbRectangle.cxx new file mode 100644 index 0000000000000000000000000000000000000000..2d633e31631f93ef592ca33a241b450425f725e6 --- /dev/null +++ b/Testing/Code/Common/otbRectangle.cxx @@ -0,0 +1,72 @@ +/*========================================================================= + + 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 "otbRectangle.h" + +#include <iostream> +#include <fstream> + +int otbRectangle(int argc, char * argv[]) +{ + const char * outfname = argv[1]; + + typedef otb::Rectangle<> RectangleType; + typedef RectangleType::ContinuousIndexType ContinuousIndexType; + typedef RectangleType::VertexListType VertexListType; + typedef VertexListType::ConstIterator IteratorType; + + // Instantiating object + RectangleType::Pointer rectangle1 = RectangleType::New(); + + ContinuousIndexType newVertex; + + newVertex[0]=atof(argv[2]); + newVertex[1]=atof(argv[3]); + rectangle1->AddVertex(newVertex); + + newVertex[0]=atof(argv[4]); + newVertex[1]=atof(argv[5]); + rectangle1->AddVertex(newVertex); + + rectangle1->SetWidth(atof(argv[6])); + rectangle1->SetOrientation(atof(argv[7])); + + + /** Inside the rectangle test*/ + ContinuousIndexType InsideVertex; + InsideVertex[0] = atof(argv[8]); + InsideVertex[1] = atof(argv[9]); + + + std::ofstream outfile(outfname); + + if(rectangle1->IsInside(InsideVertex)) + outfile <<"The point " << InsideVertex << " Is Inside the rectangle" << std::endl; + else + outfile <<"The point " << InsideVertex << " Is Outside the rectangle" << std::endl; + + + + outfile<< "region Size" << rectangle1->GetBoundingRegion().GetSize()<<std::endl; + outfile<< "region Origin" << rectangle1->GetBoundingRegion().GetIndex()<<std::endl; + + + outfile.close(); + + + return EXIT_SUCCESS; +} diff --git a/Testing/Code/Common/otbRectangleNew.cxx b/Testing/Code/Common/otbRectangleNew.cxx new file mode 100644 index 0000000000000000000000000000000000000000..a94cb60e92c4a487f38f80712893181a59f4d679 --- /dev/null +++ b/Testing/Code/Common/otbRectangleNew.cxx @@ -0,0 +1,29 @@ +/*========================================================================= + + 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 "otbRectangle.h" + +int otbRectangleNew(int argc, char * argv[]) +{ + typedef otb::Rectangle<> RectangleType; + + //Instantiating object + RectangleType::Pointer polygon = RectangleType::New(); + + return EXIT_SUCCESS; +} diff --git a/Testing/Code/MultiScale/CMakeLists.txt b/Testing/Code/MultiScale/CMakeLists.txt index 42df66bb87f24d4bc0f6058aa3b47ffb5aaf2f07..3840d35d15468ac919214ad04ee1fb1e615d6cb5 100644 --- a/Testing/Code/MultiScale/CMakeLists.txt +++ b/Testing/Code/MultiScale/CMakeLists.txt @@ -356,6 +356,10 @@ ADD_TEST(msTvGeodesicMorphologyIterativeDecompositionImageFilter ${MULTISCALE_TE 2 ) +# ------- otb::StationaryFilterBankNew ---------- +ADD_TEST(msTuStationaryFilterBankNew ${MULTISCALE_TESTS4} + otbStationaryFilterBankNew) + # ------- otb::WaveletTransform ---------- ADD_TEST(msTuWaveletTransformNew ${MULTISCALE_TESTS4} otbWaveletTransformNew) @@ -380,6 +384,8 @@ ADD_TEST(msTvWaveletTransform ${MULTISCALE_TESTS4} ADD_TEST(msTuWaveletPacketTransformNew ${MULTISCALE_TESTS4} otbWaveletPacketTransformNew) + + # ------- Fichiers sources CXX ----------------------------------- SET(BasicMultiScale_SRCS1 otbMorphologicalPyramidResamplerNew.cxx @@ -418,6 +424,7 @@ otbGeodesicMorphologyIterativeDecompositionImageFilterNew.cxx otbGeodesicMorphologyIterativeDecompositionImageFilter.cxx ) SET(BasicMultiScale_SRCS4 +otbStationaryFilterBankNew.cxx otbWaveletTransformNew.cxx otbWaveletTransform.cxx otbWaveletPacketTransformNew.cxx diff --git a/Testing/Code/MultiScale/otbMultiScaleTests4.cxx b/Testing/Code/MultiScale/otbMultiScaleTests4.cxx index b83bc40f84426eeb7a0f91388239837da46e955b..d13c69dad4af998792caf43de8578a29ba1c4e08 100644 --- a/Testing/Code/MultiScale/otbMultiScaleTests4.cxx +++ b/Testing/Code/MultiScale/otbMultiScaleTests4.cxx @@ -26,6 +26,7 @@ void RegisterTests() { + REGISTER_TEST(otbStationaryFilterBankNew); REGISTER_TEST(otbWaveletTransformNew); REGISTER_TEST(otbWaveletTransform); REGISTER_TEST(otbWaveletPacketTransformNew); diff --git a/Testing/Code/MultiScale/otbStationaryFilterBankNew.cxx b/Testing/Code/MultiScale/otbStationaryFilterBankNew.cxx new file mode 100644 index 0000000000000000000000000000000000000000..e8e6ee39e1e0e3850dbf32e1dbaa0725baf5fb9c --- /dev/null +++ b/Testing/Code/MultiScale/otbStationaryFilterBankNew.cxx @@ -0,0 +1,35 @@ +/*========================================================================= + +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 "otbImage.h" +#include "otb_9_7_Operator.h" +#include "otbStationaryFilterBank.h" + +int otbStationaryFilterBankNew(int argc, char * argv[]) +{ + const int Dimension = 2; + typedef double PixelType; + typedef otb::Image< PixelType, Dimension > ImageType; + typedef otb::LowPass_9_7_Operator< PixelType, Dimension > LowPassOperator; + typedef otb::HighPass_9_7_Operator< PixelType, Dimension > HighPassOperator; + typedef otb::StationaryFilterBank< ImageType, ImageType, LowPassOperator, HighPassOperator > FilterType; + + FilterType::Pointer filter = FilterType::New(); + + return EXIT_SUCCESS; +}