diff --git a/Code/FeatureExtraction/otbTupinEdgeDetector.h b/Code/FeatureExtraction/otbTupinEdgeDetector.h
index fbe41344fb9e4a5d5e0eff8fcb73a4b1fe5d9f36..3a5fd9d15203c0ef3dfd38cf7ba6602f3c52a173 100755
--- a/Code/FeatureExtraction/otbTupinEdgeDetector.h
+++ b/Code/FeatureExtraction/otbTupinEdgeDetector.h
@@ -22,6 +22,7 @@
(_xout) = (_x)*cos(_theta) - (_y)*sin(_theta); \
(_yout) = (_x)*sin(_theta) + (_y)*cos(_theta)
+
namespace otb
{
@@ -141,10 +142,10 @@ private:
TupinEdgeDetector(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
- /** Length of the linear feature*/
+ /** Length of the linear feature = 2*m_LengthLine+1 */
unsigned int m_LengthLine;
- /** Width of the linear feature*/
+ /** Width of the linear feature = 2*m_WidthLine+1 */
unsigned int m_WidthLine;
/** Radius of the region*/
diff --git a/Code/FeatureExtraction/otbTupinEdgeDetector.txx b/Code/FeatureExtraction/otbTupinEdgeDetector.txx
index e5506c16ccefc07796f11b0e01670285796b6e1d..baaf1daeb46a4605a05ed1374586fc08f2881786 100755
--- a/Code/FeatureExtraction/otbTupinEdgeDetector.txx
+++ b/Code/FeatureExtraction/otbTupinEdgeDetector.txx
@@ -35,6 +35,9 @@ template <class TInputImage, class TOutputImage, class InterpolatorType >
TupinEdgeDetector<TInputImage, TOutputImage, InterpolatorType>::TupinEdgeDetector()
{
m_Radius.Fill(1);
+ m_LengthLine = 1;
+ m_WidthLine = 0;
+ m_FaceList.Fill(0);
m_Interpolator = InterpolatorType::New();
}
@@ -111,14 +114,14 @@ void TupinEdgeDetector< TInputImage, TOutputImage, InterpolatorType>
typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator fit;
// Define the size of the region by the radius
- m_Radius[0] = static_cast<unsigned int>(0.5*( (3*m_WidthLine + 2) - 1 ));
- m_Radius[1] = static_cast<unsigned int>(0.5*(m_LengthLine-1));
+ m_Radius[0] = static_cast<unsigned int>(3*m_WidthLine + 2);
+ m_Radius[1] = m_LengthLine ;
//std::cout <<"Radius " << m_Radius[0] << " " << m_Radius[1] << std::endl;
// Define the size of the facelist by taking into account the rotation of the region
- m_FaceList[0] = sqrt((m_Radius[0]*m_Radius[0])+ (m_Radius[1]*m_Radius[1]));
+ m_FaceList[0] = static_cast<unsigned int>( sqrt((m_Radius[0]*m_Radius[0])+ (m_Radius[1]*m_Radius[1]))+1 );
m_FaceList[1] = m_FaceList[0];
itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> bC;
@@ -147,7 +150,10 @@ void TupinEdgeDetector< TInputImage, TOutputImage, InterpolatorType>
Theta[3] = 3*M_PI / 4. ;
// Number of the zone
- int zone;
+ unsigned int zone;
+
+ // Pixel numbers in each zone
+ const int NbPixelZone = (2*m_WidthLine+1)*(2*m_LengthLine+1);
// Contains for the 4 directions the sum of the pixels belonging to each zone
double Sum[NB_DIR][NB_ZONE];
@@ -165,7 +171,7 @@ void TupinEdgeDetector< TInputImage, TOutputImage, InterpolatorType>
// Pixel location in the input image
int X, Y;
- // Pixel location in the system axis of the region after rotation
+ // Pixel location after rotation in the system axis of the region
double xout, yout;
// Pixel location in the input image after rotation
@@ -190,13 +196,11 @@ void TupinEdgeDetector< TInputImage, TOutputImage, InterpolatorType>
bit.OverrideBoundaryCondition(&nbc);
bit.GoToBegin();
-
-//std::cout <<"***"<< neighborhoodSize << std::endl;
while ( ! bit.IsAtEnd() )
{
-std::cout << "Xc,Yc " << bit.GetIndex() << std::endl;
+//std::cout << "Xc,Yc " << bit.GetIndex() << std::endl;
// Initialisations
for (int dir=0; dir<NB_DIR; dir++)
@@ -213,22 +217,21 @@ std::cout << "Xc,Yc " << bit.GetIndex() << std::endl;
Yc = bitIndex[1];
// Location of the central pixel between zone 1 and zone 2
- Xc12 = ( Xc - (m_WidthLine-1)/2 ) - 1 ;
+ Xc12 = Xc - m_WidthLine - 1;
// Location of the central pixel between zone 1 and zone 3
- Xc13 = ( Xc + (m_WidthLine-1)/2 ) + 1 ;
+ Xc13 = Xc + m_WidthLine + 1;
// Loop on the region
for (i = 0; i < neighborhoodSize; ++i)
{
-std::cout << "---"<< i <<" "<< bit.GetIndex(i)<< std::endl;
-std::cout << "val(X,Y) "<< static_cast<double>(bit.GetPixel(i)) << " " << std::endl;
+//std::cout << "---"<< i <<" "<< bit.GetIndex(i)<< std::endl;
+//std::cout << "val(X,Y) "<< static_cast<double>(bit.GetPixel(i)) << " " << std::endl;
bitIndex = bit.GetIndex(i);
X = bitIndex[0];
Y = bitIndex[1];
-
-//std::cout << "---"<< i <<" "<< static_cast<double>(bit.GetPixel(i))<< std::endl;
+
// We determine in the vertical direction with which zone the pixel belongs.
@@ -246,33 +249,29 @@ std::cout << "val(X,Y) "<< static_cast<double>(bit.GetPixel(i)) << " " << std::e
{
ROTATION( (X-Xc), (Y-Yc), Theta[dir], xout, yout);
- // point[0] = double(xout + Xc);
- // point[1] = double(yout + Yc);
-
Index[0] = static_cast<float>(xout + Xc);
Index[1] = static_cast<float>(yout + Yc);
-std::cout << "X' Y' "<< (xout + Xc) << " " << (yout + Yc) << std::endl;
-std::cout << "val(X',Y') "<< static_cast<double>(m_Interpolator->Evaluate( point )) << std::endl;
+//std::cout << "X' Y' "<< (xout + Xc) << " " << (yout + Yc) << std::endl;
+//std::cout << "val(X',Y') "<< static_cast<double>(m_Interpolator->EvaluateAtContinuousIndex( Index )) << std::endl;
- // Sum[dir][zone] += static_cast<double>(m_Interpolator->Evaluate( point ));
Sum[dir][zone] += static_cast<double>(m_Interpolator->EvaluateAtContinuousIndex( Index ));
-
-
+
}
- } // end of the loop on the pixels of the region
-
-//std::cout << static_cast<double>(Sum[0][0])/ double(WidthZone*HeightZone) << std::endl;
+ } // end of the loop on the pixels of the region
+
+ R12 = -1.;
+ R13 = -1.;
// Loop on the 4 directions
for ( int dir=0; dir<NB_DIR; dir++ )
{
-
+
// Calculation of the averages of the 3 zones
- M1 = Sum[dir][0] / double(m_WidthLine*m_LengthLine);
- M2 = Sum[dir][1] / double(m_WidthLine*m_LengthLine);
- M3 = Sum[dir][2] / double(m_WidthLine*m_LengthLine);
+ M1 = Sum[dir][0] / static_cast<double>(NbPixelZone);
+ M2 = Sum[dir][1] / static_cast<double>(NbPixelZone);
+ M3 = Sum[dir][2] / static_cast<double>(NbPixelZone);
// Calculation of the intensity of the linear feature
if (( M1 != 0 ) && (M2 != 0))
@@ -284,23 +283,15 @@ std::cout << "val(X',Y') "<< static_cast<double>(m_Interpolator->Evaluate( point
R13_theta[dir] = static_cast<double>( 1 - MIN( (M1/M3), (M3/M1) ) );
else
R13_theta[dir] = 0.;
+
+ // Determination of the maximum intensity of the linear feature
+ R12 = static_cast<double>( MAX( R12, R12_theta[dir] ) );
+ R13 = static_cast<double>( MAX( R13, R13_theta[dir] ) );
} // end of the loop on the directions
- // Determination of the maximum intensity of the linear feature
- R12 = R12_theta[0];
- R13 = R13_theta[0];
-
- for (int dir=1; dir<NB_DIR; dir++)
- {
- R12 = static_cast<double>( MAX( R12, R12_theta[dir] ) );
- R13 = static_cast<double>( MAX( R13, R13_theta[dir] ) );
- }
-
// Intensity of the linear feature
R = MIN ( R12, R13 );
-
-//std::cout << "R = " << R << std::endl;
// Assignment of this value to the output pixel
it.Set( static_cast<OutputPixelType>(R) );