diff --git a/Code/FeatureExtraction/otbFlst.h b/Code/FeatureExtraction/otbFlst.h index 3b9963e06eaa3fb1b78ec7ba931cb867134c64f2..3e4ddaa571575ee14ba7e6fe538cc906688b3c84 100644 --- a/Code/FeatureExtraction/otbFlst.h +++ b/Code/FeatureExtraction/otbFlst.h @@ -21,23 +21,6 @@ namespace otb * \brief Algorithme de Fast Level Sets Transform of an image * */ - -/** Définition des paramètres d'optimisation. - * Ils sont dépendant de la taille des images mais ces valeurs semblent - * convenir à la majoritée des images. */ -#define INIT_MAX_AREA 10 -#define STEP_MAX_AREA 8 - -/* A 'local configuration' of the pixel is the logical 'or' of these values, -stored in one byte. Corresponding bit is set if the neighbor is in region */ -#define EAST 1 -#define NORTH_EAST 2 -#define NORTH 4 -#define NORTH_WEST 8 -#define WEST 16 -#define SOUTH_WEST 32 -#define SOUTH 64 -#define SOUTH_EAST 128 template <class TInputImage, class TOutputTree> @@ -62,6 +45,7 @@ public: /** "typedef" pour simplifier la définition et la déclaration de variables. */ + typedef typename Superclass::InputImageType InputImageType; typedef typename Superclass::InputImagePointer InputImagePointer; typedef typename Superclass::InputImageRegionType InputImageRegionType; typedef typename Superclass::InputImagePixelType InputImagePixelType; @@ -78,17 +62,43 @@ protected: virtual ~FLST() {}; void PrintSelf(std::ostream& os, itk::Indent indent) const; + /** Définition des paramètres d'optimisation. + * Ils sont dépendant de la taille des images mais ces valeurs semblent + * convenir à la majoritée des images. */ + static const int INIT_MAX_AREA; + static const int STEP_MAX_AREA; + + /** configuration locale des pixels */ + static const int EAST; + static const int NORTH_EAST; + static const int NORTH; + static const int NORTH_WEST; + static const int WEST; + static const int SOUTH_WEST; + static const int SOUTH; + static const int SOUTH_EAST; + + static const char tabPatterns[2][256]; + + itkSetMacro(MinArea, int); + itkGetConstReferenceMacro(MinArea, int); + + void CalculateArea(int Width, int Height); void GenerateData(void); - + private: FLST(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented + int m_MinArea; + int m_AreaImage; + int m_HalfAreaImage; + int m_PerimeterImage; }; } // end namespace otb #ifndef OTB_MANUAL_INSTANTIATION -#include "otbFLST.txx" +#include "otbFlst.txx" #endif diff --git a/Code/FeatureExtraction/otbFlst.txx b/Code/FeatureExtraction/otbFlst.txx index 8d1c8b69c3fce7bea45c73efd06983e3c419a92f..5eed9751ad615ccc3903185ea49a2eb1996efdf3 100644 --- a/Code/FeatureExtraction/otbFlst.txx +++ b/Code/FeatureExtraction/otbFlst.txx @@ -1 +1,148 @@ +/*========================================================================= + + Programme : OTB (ORFEO ToolBox) + Auteurs : CS - P.Imbo + Language : C++ + Date : 21 fevrier 2006 + Role : + $Id$ + +=========================================================================*/ +#ifndef __otbFlst_txx +#define __otbFlst_txx + +#include "otbFlst.h" + +namespace otb +{ + +template< class TInputImage, class TOutputTree> +const int FLST<TInputImage,TOutputTree>::INIT_MAX_AREA = 10; + +template< class TInputImage, class TOutputTree> +const int FLST<TInputImage,TOutputTree>::STEP_MAX_AREA = 8; + +template< class TInputImage, class TOutputTree> +const int FLST<TInputImage,TOutputTree>::EAST = 1; + +template< class TInputImage, class TOutputTree> +const int FLST<TInputImage,TOutputTree>::NORTH_EAST = 2; + +template< class TInputImage, class TOutputTree> +const int FLST<TInputImage,TOutputTree>::NORTH = 4; + +template< class TInputImage, class TOutputTree> +const int FLST<TInputImage,TOutputTree>::NORTH_WEST = 8; + +template< class TInputImage, class TOutputTree> +const int FLST<TInputImage,TOutputTree>::WEST = 16; + +template< class TInputImage, class TOutputTree> +const int FLST<TInputImage,TOutputTree>::SOUTH_WEST = 32; + +template< class TInputImage, class TOutputTree> +const int FLST<TInputImage,TOutputTree>::SOUTH = 64; + +template< class TInputImage, class TOutputTree> +const int FLST<TInputImage,TOutputTree>::SOUTH_EAST =128; + +template< class TInputImage, class TOutputTree> +const char FLST<TInputImage,TOutputTree>::tabPatterns[2][256] = +{{ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, + 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, + 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, + 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, + 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 1, 2, 1, 1, + 1, 2, 1, 2, 2, 3, 2, 2, 1, 2, 1, 2, 1, 2, 1, 1, + 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 1, 2, 1, 1, + 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, + 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, + 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, + 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, + 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, + 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, + 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, + 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, + 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,-1}, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, + 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, + 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, + 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, + 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, + 0, 0, 1, 0, 0,-1, 0,-1, 0, 0, 1, 0, 0,-1, 0,-1, + 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, + 0, 0, 1, 0, 0,-1, 0,-1, 0, 0, 1, 0, 0,-1, 0,-1, + 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, + 1, 1, 2, 1, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, + 1, 1, 2, 1, 2, 1, 2, 1, 2, 2, 3, 2, 2, 1, 2, 1, + 1, 1, 2, 1, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, + 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, + 0, 0, 1, 0, 0,-1, 0,-1, 0, 0, 1, 0, 0,-1, 0,-1, + 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 2, 1, 1, 0, 1, 0, + 0, 0, 1, 0, 0,-1, 0,-1, 0, 0, 1, 0, 0,-1, 0,-1}}; + + +template< class TInputImage, class TOutputTree> +FLST<TInputImage,TOutputTree>::FLST() +{ + m_MinArea = 1; +} + +template< class TInputImage, class TOutputTree> +void FLST<TInputImage,TOutputTree>::CalculateArea(int Width, int Height) +{ + int iAreaImage; + + m_AreaImage = Width * Height; + m_HalfAreaImage = (m_AreaImage+1) / 2; + if(Width == 1) m_PerimeterImage = Height; + else if(Height == 1) m_PerimeterImage = Width; + else m_PerimeterImage = (Width+Height-2)*2; + + if(m_MinArea > m_AreaImage) std::cerr << "CalculateArea: min area > image"<< std::endl; + + +} + +template< class TInputImage, class TOutputTree> +void FLST<TInputImage,TOutputTree>::GenerateData(void) +{ + int iWidth, iHeight; + int iMinArea, iMaxArea, iAreaImage, iHalfAreaImage, iPerimeterImage; + int iExploration; /* Used to avoid reinitializing images */ +// static Point_plane tabPointsInShape; + static int** tabtabVisitedNeighbors; /* Exterior boundary */ +// static Shapes pGlobalTree; + static int iAtBorder; /* #points meeting at border of image */ + typename InputImageType::SizeType Taille; + + itkDebugMacro(<< "FLST::GenerateData() called"); + // Get the input and output pointers + const InputImageType *InputImage = this->GetImageInput(); + OutputTreeListType *OutputTree = this->GetOutput(); + + Taille = InputImage->GetLargestPossibleRegion().GetSize(); + iWidth = Taille[0]; + iHeight = Taille[1]; + + this->CalculateArea(iWidth, iHeight); + + itkDebugMacro(<< "FLST::GenerateData() finished"); +} + +/** + * Standard "PrintSelf" method + */ +template< class TInputImage, class TOutputTree> +void FLST<TInputImage,TOutputTree>::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + Superclass::PrintSelf( os, indent ); + os << indent << "MinArea: " << m_MinArea << std::endl; +} + + +} // end namespace otb + + +#endif diff --git a/Code/FeatureExtraction/otbImageToPathAlignFilter.h b/Code/FeatureExtraction/otbImageToPathAlignFilter.h deleted file mode 100644 index c68a132194b4260782ab712cffc87b83cf85a703..0000000000000000000000000000000000000000 --- a/Code/FeatureExtraction/otbImageToPathAlignFilter.h +++ /dev/null @@ -1,128 +0,0 @@ -/*========================================================================= - - Programme : OTB (ORFEO ToolBox) - Auteurs : CS - P.Imbo - Language : C++ - Date : 08 fevrier 2006 - Version : - Role : - $Id: $ - -=========================================================================*/ -#ifndef __otbImageToPathAlignFilter_h -#define __otbImageToPathAlignFilter_h - -#include "itkImageSource.h" -#include "itkConceptChecking.h" -#include "itkImage.h" - -namespace otb -{ - -/** \class ImageToPathAlignFilter - * \brief Base class for filters that take a Path as input and produce an image as output. - * Base class for filters that take a Path as input and produce an image as - * output. By default, if the user does not specify the size of the output - * image, the maximum size of the path's bounding box is used. The default - * spacing of the image is given by the spacing of the input path (currently - * assumed internally to be 1.0). - */ -template <class TInputImage, class TOutputPath> -class ImageToPathAlignFilter : public itk::ImageSource<TInputImage> -{ -public: - /** Standard class typedefs. */ - typedef ImageToPathAlignFilter Self; - typedef itk::ImageSource<TInputImage> Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; - - /** Method for creation through the object factory. */ - itkNewMacro(Self); - - /** Run-time type information (and related methods). */ - itkTypeMacro(ImageToPathAlignFilter,ImageSource); - - /** ImageDimension constants */ - itkStaticConstMacro(InputImageDimension, unsigned int, - TInputImage::ImageDimension); - - /** Some convenient typedefs. */ - typedef TOutputPath OutputPathType; - typedef typename OutputPathType::Pointer OutputPathPointer; - typedef typename OutputPathType::ConstPointer OutputPathConstPointer; - typedef TInputImage InputImageType; - typedef typename Superclass::InputImageRegionType InputImageRegionType; - typedef typename InputImageType::Pointer InputImagePointer; - typedef typename InputImageType::SizeType SizeType; - typedef typename InputImageType::ValueType ValueType; - typedef typename InputImageType::PixelType PixelType; - typedef double RealType; - typedef typename itk::Image<RealType,InputImageDimension> RealImageType; - typedef typename InputImageType::Pointer RealImagePointer; - - typedef typename itk::NumericTraits<PixelType>::RealType RealType; - - - - /** Spacing (size of a pixel) of the output image. The - * spacing is the geometric distance between image samples. - * It is stored internally as double, but may be set from - * float. \sa GetSpacing() */ - virtual void SetSpacing( const double* spacing); - virtual void SetSpacing( const float* spacing); - virtual const double* GetSpacing() const; - - /** Set/Get the value for pixels on and off the path. - * By default, this filter will return a "0" image with path pixels set to 1 */ - itkSetMacro(PathValue, ValueType); - itkGetMacro(PathValue, ValueType); - itkSetMacro(BackgroundValue, ValueType); - itkGetMacro(BackgroundValue, ValueType); - - /** The origin of the output image. The origin is the geometric - * coordinates of the index (0,0,...,0). It is stored internally - * as double but may be set from float. - * \sa GetOrigin() */ - virtual void SetOrigin( const double* origin); - virtual void SetOrigin( const float* origin); - virtual const double * GetOrigin() const; - - /** Set/Get Size */ - itkSetMacro(Size,SizeType); - itkGetMacro(Size,SizeType); - -protected: - ImageToPathAlignFilter(); - ~ImageToPathAlignFilter(); - - virtual void GenerateOutputInformation(){}; // do nothing - virtual void GenerateData(); - - - SizeType m_Size; - double m_Spacing[InputImageDimension]; - double m_Origin[InputImageDimension]; - ValueType m_PathValue; - ValueType m_BackgroundValue; - - virtual void PrintSelf(std::ostream& os, itk::Indent indent) const; - -private: - ImageToPathAlignFilter(const Self&); //purposely not implemented - void operator=(const Self&); //purposely not implemented - bool m_isMeaningfulSegment; /// to get all meaningful segments (maximal or not - int m_NbGradDirection; /// Number of allowed gradient direction, default 16 - int m_NbLineDirection; /// Number of line directions to scan, default 96) - double m_MinGradNorm; /// Minimum gradient norm to define a direction, default 2. - double m_Eps; /// -log10(max. number of false alarms), default 0 - -}; - -} // end namespace otb - -#ifndef OTB_MANUAL_INSTANTIATION -#include "otbImageToPathAlignFilter.txx" -#endif - -#endif diff --git a/Code/FeatureExtraction/otbImageToPathAlignFilter.txx b/Code/FeatureExtraction/otbImageToPathAlignFilter.txx deleted file mode 100644 index 04d87364b8d656b9e7bc9d8b444a4e25d9f91360..0000000000000000000000000000000000000000 --- a/Code/FeatureExtraction/otbImageToPathAlignFilter.txx +++ /dev/null @@ -1,544 +0,0 @@ -/*========================================================================= - - Programme : OTB (ORFEO ToolBox) - Auteurs : CS - P.Imbo - Language : C++ - Date : 08 fevrier 2006 - Version : - Role : - $Id: $ - -=========================================================================*/ -#ifndef __otbImageToPathAlignFilter_txx -#define __otbImageToPathAlignFilter_txx - -#include "otbImageToPathAlignFilter.h" -#include "itkImageRegionIteratorWithIndex.h" -#include "itkPathIterator.h" -#include "itkNumericTraits.h" -#include "itkImageRegionIteratorWithIndex.h" -#include "itkImageLinearConstIteratorWithIndex.h" -#include "itkImageLinearIteratorWithIndex.h" -#include <vector> - -struct one_segment -{ - short start; /* starting position (distance from border) */ - short end; /* ending position (hence, length is end-start+1) */ - double nfa; /* number of false alarms */ - char ok; -}; - - -namespace otb -{ - -/** Constructor */ -template <class TInputImage, class TOutputPath> -ImageToPathAlignFilter<TInputImage,TOutputPath> -::ImageToPathAlignFilter() -{ - this->SetNumberOfRequiredInputs(1); - m_Size.Fill(0); - - for (unsigned int i = 0; i < InputImageDimension; i++) - { - // Set an image spacing for the user - m_Spacing[i] = 1.0; - m_Origin[i] = 0; - } - - m_PathValue = itk::NumericTraits<ValueType>::One; - m_BackgroundValue = itk::NumericTraits<ValueType>::Zero; -} - -/** Destructor */ -template <class TInputImage, class TOutputPath> -ImageToPathAlignFilter<TInputImage,TOutputPath> -::~ImageToPathAlignFilter() -{ -} - - -//---------------------------------------------------------------------------- -template <class TInputImage, class TOutputPath> -void -ImageToPathAlignFilter<TInputImage,TOutputPath> -::SetSpacing(const double* spacing) -{ - unsigned int i; - for (i=0; i<InputImageDimension; i++) - { - if ( spacing[i] != m_Spacing[i] ) - { - break; - } - } - if ( i < InputImageDimension ) - { - for (i=0; i<InputImageDimension; i++) - { - m_Spacing[i] = spacing[i]; - } - } -} - -template <class TInputImage, class TOutputPath> -void -ImageToPathAlignFilter<TInputImage,TOutputPath> -::SetSpacing(const float* spacing) -{ - unsigned int i; - for (i=0; i<InputImageDimension; i++) - { - if ( (double)spacing[i] != m_Spacing[i] ) - { - break; - } - } - if ( i < InputImageDimension ) - { - for (i=0; i<InputImageDimension; i++) - { - m_Spacing[i] = spacing[i]; - } - } -} - -template <class TInputImage, class TOutputPath> -const double * -ImageToPathAlignFilter<TInputImage,TOutputPath> -::GetSpacing() const -{ - return m_Spacing; -} - -//---------------------------------------------------------------------------- -template <class TInputImage, class TOutputPath> -void -ImageToPathAlignFilter<TInputImage,TOutputPath> -::SetOrigin(const double* origin) -{ - unsigned int i; - for (i=0; i<InputImageDimension; i++) - { - if ( origin[i] != m_Origin[i] ) - { - break; - } - } - if ( i < InputImageDimension ) - { - for (i=0; i<InputImageDimension; i++) - { - m_Origin[i] = origin[i]; - } - } -} - -template <class TInputImage, class TOutputPath> -void -ImageToPathAlignFilter<TInputImage,TOutputPath> -::SetOrigin(const float* origin) -{ - unsigned int i; - for (i=0; i<InputImageDimension; i++) - { - if ( (double)origin[i] != m_Origin[i] ) - { - break; - } - } - if ( i < InputImageDimension ) - { - for (i=0; i<InputImageDimension; i++) - { - m_Origin[i] = origin[i]; - } - } -} - -template <class TInputImage, class TOutputPath> -const double * -ImageToPathAlignFilter<TInputImage,TOutputPath> -::GetOrigin() const -{ - return m_Origin; -} - -//---------------------------------------------------------------------------- - -template <class TInputImage, class TOutputPath> -void -ImageToPathAlignFilter<TInputImage,TOutputPath> -::GenerateData(void) -{ - - int i; - - itkDebugMacro(<< "ImageToPathAlignFilter::GenerateData() called"); - - // Get the input and output pointers - const InputImageType * InputImage = this->GetInput(); - OutputPathPointer OutputPath = this->GetOutput(); - - // Generate the image - double origin[InputImageDimension]; - SizeType size; - - for(i=0;i<InputImageDimension;i++) - { - // Set Image size to the size of the path's bounding box - //size[i] = (long unsigned int)(InputObject->GetBoundingBox()->GetMaximum()[i] - // - InputObject->GetBoundingBox()->GetMinimum()[i]); - size[i]=0; - origin[i]=0; - } - - typename InputImageType::IndexType index; - index.Fill(0); - typename InputImageType::RegionType region; - - // If the size of the output has been explicitly specified, the filter - // will set the output size to the explicit size, otherwise the size from the spatial - // paths's bounding box will be used as default. - - bool specified = false; - for (i = 0; i < InputImageDimension; i++) - { - if (m_Size[i] != 0) - { - specified = true; - break; - } - } - - if (specified) - { - region.SetSize( m_Size ); - } - else - { - itkExceptionMacro( << "Currently, the user MUST specify an image size" ) - //region.SetSize( size ); - } - region.SetIndex( index ); - - InputImage->SetLargestPossibleRegion( region); // - InputImage->SetBufferedRegion( region ); // set the region - InputImage->SetRequestedRegion( region ); // - - // If the spacing has been explicitly specified, the filter - // will set the output spacing to that explicit spacing, otherwise the spacing from - // the spatial object is used as default. - - specified = false; - for (i = 0; i < InputImageDimension; i++) - { - if (m_Spacing[i] != 0) - { - specified = true; - break; - } - } - - if (specified) - { - InputImage->SetSpacing(this->m_Spacing); // set spacing - } - else - { - itkExceptionMacro( << "Currently, the user MUST specify an image spacing" ) - //OutputImage->SetSpacing(InputObject->GetIndexToObjectTransform()->GetScaleComponent()); // set spacing - } - InputImage->SetOrigin(origin); // and origin - InputImage->Allocate(); // allocate the image - - RealImagePointer AngleImage = RealImageType::New(); - AngleImage->SetRegions( InputImage->GetRequestedRegion() ); - AngleImage->CopyInformation( InputImage ); - AngleImage->Allocate(); - - typedef itk::ImageLinearIteratorWithIndex< InputImageType > IteratorType; - typedef itk::ImageLinearConstIteratorWithIndex< RealType > ConstIteratorType; - - ConstIteratorType InputIt( InputImage, InputImage->GetRequestedRegion() ); - IteratorType AngleIt( AngleImage, InputImage->GetRequestedRegion() ); -// InputIt.SetDirection(0); -// AngleIt.SetDirection(0); - - /** 1ere ETAPE: Compute the direction of the level line at each point*/ - - typename RealImageType::IndexType requestedIndex = AngleImage->GetRequestedRegion().GetIndex(); - typename RealImageType::SizeType requestedSize = AngleImage->GetRequestedRegion().GetSize(); - typename RealImageType::IndexType idx = AngleIt.GetIndex(); - RealType com1,com2,gx,gy,norm; - - for ( AngleIt.GoToBegin(); !AngleIt.IsAtEnd(); ++AngleIt) - { - typename RealImageType::IndexType idx = AngleIt.GetIndex(); - if( (idx[0]!=requestedSize[0]) || (idx[1]!=requestedSize[1]) ) - { - RealType PixelA = static_cast<RealType>(InputImage->GetPixel(idx) ); - RealType PixelB = static_cast<RealType>(InputImage->GetPixel(idx+1) ); - idx[0] = requestedIndex[0] + idx[0]; - idx[1] = requestedIndex[1] + idx[1] +1; - RealType PixelC = static_cast<RealType>(InputImage->GetPixel(idx) ); - RealType PixelD = static_cast<RealType>(InputImage->GetPixel(idx+1) ); - RealType com1 = PixelD-PixelA; - RealType com2 = PixelB-PixelC; - gx = 0.5 * (com1+com2); - gy = 0.5 * (com1-com2); - norm = gx*gx + gy*gy; - if (norm <=m_MinGradNorm) - AngleIt.Set(static_cast<RealType>(-1000.0)); - else AngleIt.Set( static_cast<RealType>( atan2(gx,-gy)) ); - } - else - { - AngleIt.Set(static_cast<RealType>(-1000.0)); - } - } - - /** 2ème ETAPE: Compute the direction of the level line at each point*/ - RealType max_nfa = pow(10.0,-m_Eps); - int nx,ny; - nx = requestedSize[0]; - ny = requestedSize[1]; - - /* maximal length for a line */ - int n = (int)( ceil(hypot((double)nx,(double)ny))+1 ); - - /*** compute P(k,l) ***/ - std::vector<float> test; - test.resize(n+1); - int adr1,adr2,x,y; - float lambda,q; - float p = 1.0/ static_cast<float>(m_NbGradDirection); - float m = (float)(nx*ny)*(float)(nx*ny); - - q = 1.0-p; - - test[0] = 1.0; - for (y=1,adr2=0;y<=n;y++) { - adr1 = adr2; - adr2 += n+1; - test[adr2] = q*test[adr1]; - for (x=1;x<=y;x++) - test[adr2+x] = p*test[adr1+x-1] + q*test[adr1+x]; - } - - /*** sum to obtain proba (>=k among y) ***/ - for (y=1,adr1=n+1;y<=n;y++,adr1+=n+1) - for (x=y-1;x>=0;x--) - test[adr1+x] += test[adr1+x+1]; - - /*** multiply by m (number of segments) to obtain expectation***/ - for (adr1=(n+1)*(n+1);--adr1>=0;) - test[adr1] *= m; - -/* Main Procedure */ - - float prec = M_PI/static_cast<float>(m_NbGradDirection); - float ntheta = m_NbLineDirection/2; /* i.e. # directions of NON-ORIENTED lines */ - float dtheta = M_PI/(float)ntheta; - - int max_nblocs = n/2+1; /* maximal number of blocs */ - std::vector<int> count,startbloc,endbloc; - - count.resize(max_nblocs); - startbloc.resize(max_nblocs); - endbloc.resize(max_nblocs); - - int size_seg = 10000; /* initial allocation (may reallocate later) */ -// seg = (struct one_segment *)malloc(size_seg*sizeof(struct one_segment)); - - /* counter for recorded segments (seglist) */ - int iseglist = 0; - -/******************** first loop : the four sides ********************/ - int mx,my,ox,oy; - int xx,yy,pos,posmax,nblocs,inbloc; - int cur,j,side,tmp,l,lphase; - int itheta; - int iseg,size_seglist; - float theta,theta0,dx,dy,error; - double nfa; - - std::vector<one_segment> seg; - std::vector<float> seglist; - - seg.resize(size_seg); - seglist.resize(5*size_seglist); - - for (side=0;side<4;side++) { - - theta0 = 0.5*M_PI*(double)side; - mx = ((side==0 || side==2)?1:0); - my = ((side==1 || side==3)?1:0); - ox = ((side==1)?nx-1:0); - oy = ((side==2)?ny-1:0); - - posmax = nx*mx+ny*my; - - - /*** second loop : angles ***/ - for (itheta = 0; itheta<ntheta; itheta++) { - theta = theta0 + (float)(itheta)*dtheta; - dx = (float)cos((double)theta); - dy = (float)sin((double)theta); - - /*** third loop : start positions ***/ - for (pos=0;pos<posmax;pos++) { - - /* clear segment array */ - iseg = 0; - - /*** fourth loop : phase for two-spaced pixels ***/ - for (lphase=0;lphase<2;lphase++) { - - /*** detect aligned points by blocs ***/ - inbloc = nblocs = cur = l = count[0] = 0; - xx = ox+pos*mx + (int)(dx*(float)(l*2+lphase)); - yy = oy+pos*my + (int)(dy*(float)(l*2+lphase)); - - for (;xx>=0 && xx<nx && yy>=0 && yy<ny;) { - - error = static_cast<float>( AngleImage->GetPixel(yy*nx+xx) ); - if (error>-100.0) { - error -= theta; - while (error<=-M_PI) error += 2.0*M_PI; - while (error>M_PI) error -= 2.0*M_PI; - if (error<0.0) error = -error; - if (error<prec) { - cur++; - if (!inbloc) { - startbloc[nblocs]=l; - inbloc=1; - } - } else { - if (inbloc) { - endbloc[nblocs] = l-1; - nblocs++; - count[nblocs] = cur; - } - inbloc=0; - } - } - /* compute next point */ - l++; - xx = ox+pos*mx + (int)(dx*(float)(l*2+lphase)); - yy = oy+pos*my + (int)(dy*(float)(l*2+lphase)); - } - - /*** detect meaningful segments ***/ - for (i=0;i<nblocs;i++) - for (j=i;j<nblocs;j++) - if ((nfa = test[count[j+1]-count[i] - +(n+1)*(1+endbloc[j]-startbloc[i])]) < max_nfa) { - seg[iseg].start = startbloc[i]*2+lphase; - seg[iseg].end = endbloc[j]*2+lphase; - seg[iseg].nfa = nfa; - seg[iseg].ok = 1; - iseg++; - /* reallocate if necessary */ - if (iseg==size_seg) { - size_seg = (size_seg*3)/2; - seg.resize(size_seg); -// if (!seg) -// mwerror(FATAL,1,"Not enough memory."); - } - } - } - /*** end of phase loop ***/ - - /*** remove non-maximal segments ***/ - if (!m_isMeaningfulSegment) - for (i=0;i<iseg;i++) - for (j=0;j<iseg;j++) - if (i!=j) - - /* seg[i] is included in seg[j] ? */ - if (seg[i].start>=seg[j].start && seg[i].end<=seg[j].end) { - - /* remove the less meaningful of seg[i] and seg[j] */ - if (seg[i].nfa<seg[j].nfa) seg[j].ok=0; - else seg[i].ok=0; - - } - - /*** store detected segments ***/ - for (i=0;i<iseg;i++) - if (seg[i].ok) { - seglist[iseglist*5 ]=(float)(ox+pos*mx)+dx*(float)(seg[i].start); - seglist[iseglist*5+1]=(float)(oy+pos*my)+dy*(float)(seg[i].start); - seglist[iseglist*5+2]=(float)(ox+pos*mx)+dx*(float)(seg[i].end); - seglist[iseglist*5+3]=(float)(oy+pos*my)+dy*(float)(seg[i].end); - seglist[iseglist*5+4]=-(float)log10(seg[i].nfa); - iseglist++; - /* reallocate seglist if necessary */ - if (iseglist==size_seglist) { - size_seglist = (size_seglist*3)/2; - seglist.resize(size_seglist); -// if (!seglist) -// mwerror(FATAL,1,"Not enough memory."); - } - } - } - } - /*** end of second loop ***/ - - /******************** end of first loop ********************/ - - - /* build segments list */ - seglist.resize(5*iseglist); -// result = mw_new_flist(); -// result->size = result->max_size = iseglist; -// result->dim = 5; -// result->values = seglist; - - /* build curves if requested */ -// if (crv) { -// crv = mw_change_flists(crv,iseglist,iseglist); -// for (i=0;i<iseglist;i++) { -// crv->list[i] = mw_change_flist(NULL,2,2,2); -// for (j=0;j<4;j++) -// crv->list[i]->values[j] = seglist[i*5+j]+.5; -// } -// } - - - - } - - - - -// itk::PathIterator<InputImageType,OutputPathType> pathIt(InputImage,OutputPath); -// for( pathIt.GoToBegin(); !pathIt.IsAtEnd(); ++pathIt ) -// { -// pathIt.Set( m_PathValue ); -// } - - itkDebugMacro(<< "ImageToPathAlignFilter::GenerateData() finished"); - -} // end update function - - -template <class TInputImage, class TOutputPath> -void -ImageToPathAlignFilter<TInputImage,TOutputPath> -::PrintSelf(std::ostream& os, itk::Indent indent) const -{ - Superclass::PrintSelf(os, indent); - os << indent << "Size : " << m_Size << std::endl; - os << indent << "Path Value : " << m_PathValue << std::endl; - os << indent << "Background Value : " << m_BackgroundValue << std::endl; -} - - - -} // end namespace otb - -#endif diff --git a/Code/FeatureExtraction/otbTreeNeighborhood.h b/Code/FeatureExtraction/otbTreeNeighborhood.h index 1ffc6a761394adda0cb1be6bbd5d686ac2dad9e9..fa5603b9cff6c2af03d90cdee9f98fc03a1b1f34 100644 --- a/Code/FeatureExtraction/otbTreeNeighborhood.h +++ b/Code/FeatureExtraction/otbTreeNeighborhood.h @@ -4,7 +4,7 @@ Auteurs : CS - P.Imbo Language : C++ Date : 20 fevrier 2006 - Role : Definition de la classe TreeSource + Role : Definition de la classe TreeNeighborhood $Id$ =========================================================================*/ @@ -31,7 +31,7 @@ namespace otb * \ingroup TreeNeighborhoodObjects */ -class otbTreeNeighborhood : public itk::DataObject +class TreeNeighborhood : public itk::DataObject { public: /** Standard class typedefs. */ @@ -46,17 +46,32 @@ public: /** Run-time type information (and related methods). */ itkTypeMacro(TreeSource,itk::FunctionBase); - enum TypeOfTree {AMBIGUOUS, MIN, MAX, INVALID}; + typedef enum {AMBIGUOUS, MIN, MAX, INVALID} TypeOfTree; - itk::Point<float,2> PointType; - std::vector<PointType> PointListType; + typedef struct { + itk::Point<float,2> Point; + int x; + int y; + float Value; + } PointType; + + typedef std::vector<PointType> PointListType; protected: TreeNeighborhood(); - virtual ~TreeNeighborhood() {m_tabPoints.clear();} + virtual ~TreeNeighborhood() {} + void SetSize(int Taille); + const int ORDER_MAX(int k,int l); + const int ORDER_MIN(int k,int l); + const void SWAP(int k,int l); + const int ORDER_MAX2(int k,int l); + const int ORDER_MIN2(int k,int l); + void FixUp(); + void FixDown(); + void Add(int x, int y,float value); + void Remove(); void PrintSelf(std::ostream& os, itk::Indent indent) const; - enum TypeOfTree {AMBIGUOUS, MIN, MAX, INVALID}; private: TreeNeighborhood(const Self&); //purposely not implemented diff --git a/Code/FeatureExtraction/otbTreeNeighborhood.txx b/Code/FeatureExtraction/otbTreeNeighborhood.txx index 2f62a3967d0c7317ed458b2188c6183ddcc63b4a..570b07c25030359de45efe80b3dcff16003d3feb 100644 --- a/Code/FeatureExtraction/otbTreeNeighborhood.txx +++ b/Code/FeatureExtraction/otbTreeNeighborhood.txx @@ -22,8 +22,7 @@ namespace otb * */ -otbTreeNeighborhood -::TreeNeighborhood() +TreeNeighborhood::TreeNeighborhood() { m_tabPoints.clear(); m_iNbPoints = 0; @@ -31,21 +30,165 @@ otbTreeNeighborhood m_otherBound = 0.0; } -void -otbTreeNeighborhood -::SetSize(int Taille) +void TreeNeighborhood::SetSize(int Taille) { assert( Taille > 0 ); m_tabPoints.resize(Taille); } + +const int TreeNeighborhood::ORDER_MAX(int k,int l) +{ + int Taille = m_tabPoints.size(); + assert( k > 0); + assert( l > 0); + assert( k<= Taille); + assert( l>= Taille); + return (m_tabPoints[k].Value > m_tabPoints[l].Value); +} +const int TreeNeighborhood::ORDER_MIN(int k,int l) +{ + int Taille = m_tabPoints.size(); + assert( k > 0); + assert( l > 0); + assert( k<= Taille); + assert( l>= Taille); + return (m_tabPoints[k].Value < m_tabPoints[l].Value); +} + +const void TreeNeighborhood::SWAP(int k,int l) +{ + int Taille = m_tabPoints.size(); + assert( k > 0); + assert( l > 0); + assert( k<= Taille); + assert( l>= Taille); + m_tabPoints[0] = m_tabPoints[k]; + m_tabPoints[k] = m_tabPoints[l]; + m_tabPoints[l] = m_tabPoints[0]; +} + +const int TreeNeighborhood::ORDER_MAX2(int k,int l) +{ + int Taille = m_tabPoints.size(); + assert( k > 0); + assert( l > 0); + assert( k<= Taille); + assert( l>= Taille); + return (m_tabPoints[k].Value >= m_tabPoints[l].Value); + +} + +const int TreeNeighborhood::ORDER_MIN2(int k,int l) +{ + int Taille = m_tabPoints.size(); + assert( k > 0); + assert( l > 0); + assert( k<= Taille); + assert( l>= Taille); + return (m_tabPoints[k].Value <= m_tabPoints[l].Value); +} + +void TreeNeighborhood::FixUp() +{ + int k,l; + k = m_iNbPoints; + if (m_type == MAX) + while(k>1 && this->ORDER_MAX(k,l=k>>1)) + { + this->SWAP(k,l); + k = l; + } + else + while(k>1 && this->ORDER_MIN(k,l=k>>1)) + { + this->SWAP(k,l); + k=l; + } +} + +void TreeNeighborhood::FixDown() +{ + int k,l,N; + N = m_iNbPoints; + k = 1; + if (m_type == MAX) + while( (l = k << 1) <= N) + { + if(l < N && this->ORDER_MAX(l+1,l)) ++l; + if(this->ORDER_MAX2(k,l)) + break; + this->SWAP(k, l); + k = l; + } + else + while((l = k << 1) <= N) + { + if(l < N && this->ORDER_MIN(l+1,l)) ++l; + if(this->ORDER_MIN2(k,l)) + break; + this->SWAP(k, l); + k = l; + } +} + +void TreeNeighborhood::Add(int x, int y,float value) +{ + if(m_iNbPoints == 0) + m_otherBound = value; + else + switch(m_type) + { + case MIN: + if(value > m_otherBound) + m_otherBound = value; + else if(value < m_tabPoints[1].Value) + m_type = INVALID; + break; + case MAX: + if(value < m_otherBound) + m_otherBound = value; + else if(value > m_tabPoints[1].Value) + m_type = INVALID; + break; + case AMBIGUOUS: + if(value != m_tabPoints[1].Value) { + m_type = (value < m_tabPoints[1].Value)? MAX: MIN; + m_otherBound = value; + } + break; + } + if(m_type == INVALID) + return; + /* 2) Add the point in the heap and update it */ + ++m_iNbPoints; + m_tabPoints[m_iNbPoints].x = x; + m_tabPoints[m_iNbPoints].x = y; + m_tabPoints[m_iNbPoints].Value = value; + + this->FixUp(); /* Update the heap of neighbors */ +} + +void TreeNeighborhood::Remove() +{ + int Taille = m_tabPoints.size(); + assert(Taille >0); + assert(Taille = m_iNbPoints); + float value,valueTop; + value = m_tabPoints[1].Value; + + if(m_type == INVALID) + return; + valueTop = m_tabPoints[m_iNbPoints--].Value; + if(m_iNbPoints == 0) + return; + this->FixDown(); + if(value != valueTop && valueTop == m_otherBound) + m_type = AMBIGUOUS; +} -/** - * - */ -void -otbTreeNeighborhood +void TreeNeighborhood ::PrintSelf(std::ostream& os, itk::Indent indent) const { Superclass::PrintSelf(os,indent); diff --git a/Code/FeatureExtraction/otbTreeSource.h b/Code/FeatureExtraction/otbTreeSource.h index 3956d95b80d228240f98b504d09ce6797e078bf8..fd1b47836239afbeb8ddec7bed10d008759b94cc 100644 --- a/Code/FeatureExtraction/otbTreeSource.h +++ b/Code/FeatureExtraction/otbTreeSource.h @@ -17,6 +17,7 @@ #endif #include "itkProcessObject.h" +#include "otbTreeNeighborhood.h" namespace otb { @@ -45,10 +46,10 @@ public: /** Run-time type information (and related methods). */ itkTypeMacro(TreeSource,itk::ProcessObject); - typedef TOutputTree OutputTreeType; + typedef TreeNeighborhood OutputTreeType; typedef typename OutputTreeType::Pointer OutputTreePointerType; typedef typename OutputTreeType::ConstPointer OutputTreeConstPointerType; - typedef std::vector< OutputTreePointerType > OutputTreeListType; + typedef OutputTreeType OutputTreeListType; /** Get the tree output of this process object. */ OutputTreeListType * GetOutput(void) @@ -62,25 +63,6 @@ protected: TreeSource(); virtual ~TreeSource() {} void PrintSelf(std::ostream& os, itk::Indent indent) const; - - enum TypeOfTree {AMBIGUOUS, MIN, MAX, INVALID}; - - class point_plane { - short x,y; /*** Coordinates of the point */ - }; - - class Neighbor { - point_plane point; /*** Neighbor pixel */ - float value; /*** Its gray level */ - }; - - class Neighborhood { - Neighbor* tabPoints; /*** The array of neighbors, organized as a binary tree */ - int iNbPoints; /*** The size of the previous arrays */ - enum TypeOfTree type; /*** max- or min- oriented heap? */ - float otherBound; /*** Min gray level if max-oriented, max if min-oriented */ - }; - private: TreeSource(const Self&); //purposely not implemented diff --git a/Code/FeatureExtraction/otbTreeSource.txx b/Code/FeatureExtraction/otbTreeSource.txx index fb97e35a151740b84ea830174a19585bd42f837a..67f443589c52bd2aee42dc6149c8c56d28ab56fe 100644 --- a/Code/FeatureExtraction/otbTreeSource.txx +++ b/Code/FeatureExtraction/otbTreeSource.txx @@ -24,7 +24,6 @@ template<class TOutputTree> TreeSource<TOutputTree> ::TreeSource() { - m_TreeList.clear(); this->ProcessObject::SetNumberOfRequiredOutputs(1); diff --git a/Testing/Code/FeatureExtraction/otbFlstTest.cxx b/Testing/Code/FeatureExtraction/otbFlstTest.cxx index 1bc4412c970d27a00d764b0239b42d8158b18439..06a53af54b43c7cd7673a23064b398897b09ef31 100644 --- a/Testing/Code/FeatureExtraction/otbFlstTest.cxx +++ b/Testing/Code/FeatureExtraction/otbFlstTest.cxx @@ -19,6 +19,8 @@ #include "otbImageFileReader.h" #include "otbTreeSource.h" #include "otbImageToTreeFilter.h" +#include "otbTreeNeighborhood.h" +#include "otbFlst.h" int otbFlstTest( int argc, char ** argv ) { @@ -37,9 +39,13 @@ int otbFlstTest( int argc, char ** argv ) typedef itk::Image< InputPixelType, Dimension > InputImageType; typedef itk::Image< RealPixelType, Dimension > RealImageType; typedef otb::ImageFileReader< InputImageType > ReaderType; - + typedef otb::TreeNeighborhood TreeType; + typedef otb::FLST<InputImageType,TreeType> FlstType; + ReaderType::Pointer reader = ReaderType::New(); reader->SetFileName( inputFilename ); + + FlstType::Pointer tree; } catch( itk::ExceptionObject & err )