diff --git a/Code/Common/otbLabelObjectToPolygonFunctor.txx b/Code/Common/otbLabelObjectToPolygonFunctor.txx
index 3b1cdc85c0f3c6954ba589519bb402798e033705..e51f63c1f801c6dabd94fdbb67da2c343dec7a53 100644
--- a/Code/Common/otbLabelObjectToPolygonFunctor.txx
+++ b/Code/Common/otbLabelObjectToPolygonFunctor.txx
@@ -472,7 +472,7 @@ void
LabelObjectToPolygonFunctor<TLabelObject,TPolygon>
::WalkLeft(unsigned int line,const IndexType & startPoint, const IndexType & endPoint, PolygonType * polygon, const StateType state)
{
- //assert( vcl_abs(static_cast<long int>(line+m_LineOffset-endPoint[1]))<=1 && "End point not with +/-1 line from line");
+ assert( vcl_abs(static_cast<long int>(line+m_LineOffset-endPoint[1]))<=1 && "End point not with +/-1 line from line");
typename PolygonType::VertexType::VectorType offset;
@@ -538,7 +538,7 @@ void
LabelObjectToPolygonFunctor<TLabelObject,TPolygon>
::WalkRight(unsigned int line,const IndexType & startPoint, const IndexType & endPoint, PolygonType * polygon, const StateType state)
{
- //assert( vcl_abs(static_cast<long int>(line+m_LineOffset-endPoint[1]))<=1 && "End point not with +/-1 line from line");
+ assert( vcl_abs(static_cast<long int>(line+m_LineOffset-endPoint[1]))<=1 && "End point not with +/-1 line from line");
typename PolygonType::VertexType::VectorType offset;
diff --git a/Code/FeatureExtraction/otbImageToHessianDeterminantImageFilter.h b/Code/FeatureExtraction/otbImageToHessianDeterminantImageFilter.h
index 07e3169616b6ab03963a95fe39df9d8a05558880..1e97f295a210ddf2dc9e111ced1a60c69f3585bc 100644
--- a/Code/FeatureExtraction/otbImageToHessianDeterminantImageFilter.h
+++ b/Code/FeatureExtraction/otbImageToHessianDeterminantImageFilter.h
@@ -56,10 +56,7 @@ public:
*/
inline TOutput operator()(const TInput& input)
{
- if(input[0]*input[1] - 0.9*input[2]*input[2] < 1e-7)
- return 0.;
-
- return static_cast<TOutput>(input[0]*input[1] - 0.9*input[2]*input[2]);
+ return static_cast<TOutput>(input[0]*input[1] - input[2]*input[2]);
}
diff --git a/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilter.h b/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilter.h
index 6dabbb21cebafb7028beb4358dcc8db539517fe9..018dcfeefa13e641b13ae8971be8ebc6b8b43d7a 100644
--- a/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilter.h
+++ b/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilter.h
@@ -171,7 +171,7 @@ protected:
* \return key point orientation
*/
virtual double AssignOrientation(const NeighborhoodType& neigh,
- double S, IndexType cindex);
+ double S);
/** ComputeDescriptor
*
@@ -183,29 +183,12 @@ protected:
*/
virtual VectorType ComputeDescriptor(const NeighborhoodType& neigh,
double O,
- double S, IndexType cindex);
+ double S);
/**
* Compute min a b c
*/
virtual int GetMin( int a , int b , int c);
- /**
- * Compute Integral Images
- */
- virtual void IntegralImage();
-
- /**
- * HaarX & HaarY : Response To HaarWavelet
- */
- virtual double HaarX(unsigned int row , unsigned int col,double S);
- virtual double HaarY(unsigned int row , unsigned int col,double S);
-
- virtual double BoxIntegral( int row, int col, int nbRow , int nbCol);
-
- /**
- * Getangle
- */
- virtual double GetAngle(unsigned int X, unsigned int Y);
private:
ImageToSURFKeyPointSetFilter(const Self&); //purposely not implemented
@@ -225,7 +208,6 @@ private:
InputImagePointerType m_ImageCurrent;
InputImagePointerType m_ImageMovedPrev;
InputImagePointerType m_ImageMovedNext;
- InputImagePointerType m_IntegralImage;
/** ImageToDeterminantHessianFilter filter */
DetHessianPointerFilter m_DetHessianFilter;
@@ -244,4 +226,3 @@ private:
#endif
-
diff --git a/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilter.txx b/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilter.txx
index 96f40fc6cbbaee0f0961fdfaefc1ef75b8467e8d..6d9b45aa205ab858228db126efac565e330925df 100644
--- a/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilter.txx
+++ b/Code/FeatureExtraction/otbImageToSURFKeyPointSetFilter.txx
@@ -20,8 +20,6 @@ PURPOSE. See the above copyright notices for more information.
#include "otbImageToSURFKeyPointSetFilter.h"
#include "itkCenteredRigid2DTransform.h"
-#include "itkImageRegionIterator.h"
-#include "otbImageFileWriter.h"
namespace otb
{
@@ -36,7 +34,6 @@ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet >
m_ScalesNumber = 3;
m_NumberOfPoints = 0;
m_DetHessianFilter = ImageToDetHessianImageType::New();
- m_IntegralImage = InputImageType::New();
}
@@ -58,19 +55,11 @@ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
::GenerateData(void)
{
- /**Integral Image */
- std::cout << "Begin Integral Image"<< std::endl;
- this->IntegralImage();
- std::cout << "Done"<< std::endl;
double k;
double sigma_in = 2.;
SizeType radius;
- /*New Adds*/
- double fscale = 1.;
- double initSigma = 1.6;
-
/*Output*/
OutputPointSetPointerType outputPointSet = this->GetOutput();
@@ -84,13 +73,12 @@ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
/* Computation loop over octaves*/
- for (int i = 0; i < m_OctavesNumber; ++i )
+ for (int i = 0; i < m_OctavesNumber; i++ )
{
sigma_in = 2.;
m_ImageList = ImageListType::New();
-
/*--------------------------------------------------------
Octave per octave
--------------------------------------------------------*/
@@ -120,7 +108,7 @@ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
}
- for (int j = 0; j < m_ScalesNumber; ++j )
+ for (int j = 0; j < m_ScalesNumber; j++ )
{
/** Incrementation of the gaussian width
* the width of the gaussian have to be doubled for
@@ -159,7 +147,7 @@ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
/* extremum Search over octave's scales */
/*----------------------------------------------------*/
- for (int jj = 1; jj < (int)(m_ImageList->Size() - 1 ); ++jj)
+ for (int jj = 1; jj < (int)(m_ImageList->Size() - 1 ); jj++)
{
m_ImageCurrent = m_ImageList->GetNthElement(jj);
m_ImageMovedPrev = m_ImageList->GetNthElement(jj-1);
@@ -194,24 +182,25 @@ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
keyPoint[1] = it.GetIndex()[1];
//keyPoint[2] = sigma_in/pow(k,(double)jj)*pow(2.,(double)i);
- //double sigmaDetected = sigma_in/pow(k,(double)jj)*pow(2.,(double)i);
+ double sigmaDetected = sigma_in/pow(k,(double)jj)*pow(2.,(double)i);
+
+ radius.Fill(GetMin((int)(this->GetInput()->GetLargestPossibleRegion().GetSize()[0] - keyPoint[0]),
+ (int)(this->GetInput()->GetLargestPossibleRegion().GetSize()[1] - keyPoint[1]),
+ (int)(6*sigmaDetected) ) ); // changer le sigma detected par keypoint[2]
- double fSize = initSigma*pow(2.,((double)i+1.)/m_OctavesNumber);
- double sigmaDetected = fSize * fscale;
- radius.Fill((int)sigmaDetected*6);
/*
Computing the orientation of the key point detected
*/
- NeighborhoodIteratorType itNeighOrientation(radius,m_IntegralImage/*this->GetInput()*/ ,
- /*this->GetInput()*/m_IntegralImage->GetLargestPossibleRegion());
+ NeighborhoodIteratorType itNeighOrientation(radius,this->GetInput() ,
+ this->GetInput()->GetLargestPossibleRegion());
itNeighOrientation.SetLocation(it.GetIndex());
/** TO DO*/
//keyPoint[3] = AssignOrientation( itNeighOrientation.GetNeighborhood() ,keyPoint[2] );
- double orientationDetected = AssignOrientation( itNeighOrientation.GetNeighborhood() , sigmaDetected , it.GetIndex());
+ double orientationDetected = AssignOrientation( itNeighOrientation.GetNeighborhood() , sigmaDetected );
/*Filling the Point pointSet Part*/
outputPointSet->SetPoint(m_NumberOfPoints, keyPoint);
@@ -220,7 +209,10 @@ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
/*----------------------------------------*/
/* Descriptor Computation */
/*----------------------------------------*/
- radius.Fill(vcl_floor(sigmaDetected*10));
+
+ radius.Fill(GetMin((int)(this->GetInput()->GetLargestPossibleRegion().GetSize()[0] - keyPoint[0]),
+ (int)(this->GetInput()->GetLargestPossibleRegion().GetSize()[1] - keyPoint[1]),
+ (int)(10*sigmaDetected))); // TODO a changer sigmaDetected par Keypoint[2]
NeighborhoodIteratorType itNeighDescriptor(radius,this->GetInput(),
this->GetInput()->GetLargestPossibleRegion());
@@ -228,7 +220,9 @@ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
VectorType descriptor;
descriptor.resize(64);
//descriptor = ComputeDescriptor(itNeighDescriptor.GetNeighborhood(),keyPoint[3],keyPoint[2]);
- descriptor = ComputeDescriptor(itNeighDescriptor.GetNeighborhood(),orientationDetected,sigmaDetected, it.GetIndex());
+ descriptor = ComputeDescriptor(itNeighDescriptor.GetNeighborhood(),orientationDetected,sigmaDetected);
+
+
/*Updating the pointset data with values of the descriptor*/
OutputPixelType data;
@@ -242,12 +236,12 @@ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
while (itDescriptor != descriptor.end())
{
data.SetElement(IndDescriptor, *itDescriptor);
- ++IndDescriptor;
- ++itDescriptor;
+ IndDescriptor++;
+ itDescriptor++;
}
outputPointSet->SetPointData(m_NumberOfPoints, data);
- ++m_NumberOfPoints;
+ m_NumberOfPoints++;
}
++it;
++itNeighPrev;
@@ -255,8 +249,6 @@ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
}/*End while for extremum search*/
- fscale += fscale;
-
} /*End Iteration over scales */
m_ImageList->Clear();
@@ -356,126 +348,81 @@ ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
template <class TInputImage, class TOutputPointSet>
double
ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
-::AssignOrientation(const NeighborhoodType& neigh , double S, IndexType cindex)
+::AssignOrientation(const NeighborhoodType& neigh , double S)
{
- SizeType size = this->GetInput()->GetLargestPossibleRegion().GetSize();
-
- unsigned int scale = vcl_floor(S + 0.5 );
- std::vector<double> resX, resY, Angs;
- for(int i = -6; i <= 6; ++i)
+ int i= 0;
+ int pas =( (i+S)-(int)(i+S) > 0.5 )?((int)S+1):(int)S;
+ int Largeur = 2*neigh.GetRadius()[0]+1; // Width & length of a neighborhood
+ int rayon = neigh.GetRadius()[0]; // radius of the neigh
+ int col, raw;
+ double dist;
+ double w; // weight of the circular gaussian
+
+ OutputPointType pt;
+
+ // Gradient orientation histogram
+ double angle;
+ int bin = 0;
+ int Pi = 180;
+ int LengthBin = 60;
+ int NbBins = (2*Pi/LengthBin);
+ std::vector<double> tab(NbBins*2 , 0.);
+
+ while (i < (int)neigh.Size())
{
- for(int j = -6; j <= 6; ++j)
+ col = i%Largeur - rayon;
+ raw = i/Largeur - rayon;
+ dist = vcl_sqrt(static_cast<double>(col *col + raw * raw) );
+ col +=rayon;
+ raw +=rayon; // Backup to the image coordinate axes
+
+ if (dist < 6*S)
{
- if(i*i + j*j < 36)
+ // Haar Wavelets responses accumulated in an histogram with Pi/3 precison
+ if (( col > pas && col < Largeur - pas ) && ( raw > pas && raw < Largeur - pas) )
{
- double w = 1./(2.*CONST_PI*4*S*S)*vcl_exp(-( i*i+j*j ) / (2.*4.*S*S) );
- resX.push_back( w*this->HaarX(cindex[1] + j*scale ,cindex[0] + i*scale , 4*scale));
- resY.push_back( w*this->HaarY(cindex[1] + j*scale ,cindex[0] + i*scale , 4*scale));
- Angs.push_back(GetAngle(resX.back(), resY.back()));
- }
- }
- }
- //Calculate the dominant direction
- double sumX, sumY;
- double max=0., orientation = 0.;
- double ang1, ang2, ang;
+ w = vcl_exp(-((col-rayon)*(col-rayon) + (raw-rayon)*(raw-rayon))/(2*2.5*2.5*S*S) );
+ pt[0] = (neigh[(col+pas) + raw * Largeur] - neigh[(col-pas) + raw *Largeur ]) * w;
+ pt[1] = (neigh[col + (raw+pas)* Largeur ] - neigh[col + (raw-pas)*Largeur]) * w;
- // loop slides pi/3 window around feature point
- for(ang1 = 0; ang1 < 2*CONST_PI; ang1 += 0.15)
- {
- ang2 = ( (ang1 + CONST_PI/3. > 2*CONST_PI ) ? ang1 - 5.*CONST_PI/3. : ang1+CONST_PI/3.0f);
- sumX = sumY = 0;
- for(unsigned int k = 0; k < Angs.size(); k++)
- {
- // get angle from the x-axis of the sample point
- ang = Angs[k];
-
- // determine whether the point is within the window
- if (ang1 < ang2 && ang1 < ang && ang < ang2)
- {
- sumX+=resX[k];
- sumY+=resY[k];
- }
- else if (ang2 < ang1 &&
- ((ang > 0 && ang < ang2) || (ang > ang1 && ang < 2*CONST_PI) ))
+ if (pt[0] + pt[1] != 0)
{
- sumX+=resX[k];
- sumY+=resY[k];
+ angle = atan( pt[0]/pt[1] )*( Pi/M_PI);
+ if (angle < 0 )
+ angle += 2*Pi;
+
+ bin = (int)(angle/LengthBin);
+
+ if ( bin <= NbBins-1 || bin >= 0 )
+ {
+ tab[2*bin] += pt[0];
+ tab[2*bin+1] += pt[1];
+ }
}
+ }
}
+ i+= pas;
+ }
- // if the vector produced from this window is longer than all
- // previous vectors then this forms the new dominant direction
- if (sumX*sumX + sumY*sumY > max)
+ //Find Orientation
+ double indice = 0;
+ double max = 0;
+ double length = 0;
+
+ //Detection de l'orientation du point courant
+ for (int i = 0; i < NbBins*2; i = i+2)
+ {
+ length = vcl_sqrt( tab[i]*tab[i] + tab[i+1]*tab[i+1] );
+ if ( length > max)
{
- // store largest orientation
- max = sumX*sumX + sumY*sumY;
- orientation = GetAngle(sumX, sumY);
+ max = length;
+ indice = i/2;
}
- }
- return orientation;
-}
-
-
-/*-----------------------------------------------------------
- * HaarX wavelet reponse in the X axis
- -----------------------------------------------------------*/
-
-template <class TInputImage, class TOutputPointSet>
-double
-ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
-::HaarX( unsigned int row , unsigned int col, double S )
-{
- return this->BoxIntegral(row-S/2, col, S, S/2) -1 * this->BoxIntegral(row-S/2, col-S/2, S, S/2);
-}
-
-// /*-----------------------------------------------------------
-// * HaarY wavelet reponse in the Y axis
-// -----------------------------------------------------------*/
-
-template <class TInputImage, class TOutputPointSet>
-double
-ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
-::HaarY(unsigned int row , unsigned int col, double S )
-{
- return BoxIntegral(row, col-S/2, S/2, S) -1 * BoxIntegral(row-S/2, col-S/2, S/2, S);
+ }
-}
-
-
-/*BoxIntegral */
-
-template <class TInputImage, class TOutputPointSet>
-double
-ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
-::BoxIntegral( int row, int col, int nbRow , int nbCol)
-{
- //Size of the image
- InputImagePointerType intimage = m_IntegralImage;
- SizeType size = intimage->GetLargestPossibleRegion().GetSize();
-
- // The subtraction by one for row/col is because row/col is inclusive.
- int r1 = std::min((double)row, (double)size[1]) - 1;
- int c1 = std::min((double)col, (double)size[0]) - 1;
- int r2 = std::min((double)row + nbRow, (double)size[1]) - 1;
- int c2 = std::min((double)col + (double)nbCol, (double)size[0]) - 1;
-
- double A(0.), B(0.), C(0.), D(0.);
-
- IndexType i1,i2,i3,i4;
- i1[1] = r1; i1[0] = c1;
- i2[1] = r1; i2[0] = c2;
- i3[1] = r2; i3[0] = c1;
- i4[1] = r2; i4[0] = c2;
-
- if (r1 >= 0 && c1 >= 0) A = intimage->GetPixel(i1);//data[r1 * size[0] + c1];
- if (r1 >= 0 && c2 >= 0) B = intimage->GetPixel(i2);//data[r1 * size[0] + c2];
- if (r2 >= 0 && c1 >= 0) C = intimage->GetPixel(i3);//data[r2 * size[0] + c1];
- if (r2 >= 0 && c2 >= 0) D = intimage->GetPixel(i4);//data[r2 * size[0] + c2];
-
- return std::max(0., A - B - C + D);
+ return (indice+0.5)*LengthBin;
}
/*-----------------------------------------------------------
@@ -486,348 +433,98 @@ template <class TInputImage, class TOutputPointSet>
typename ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
::VectorType
ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
-::ComputeDescriptor(const NeighborhoodType& neigh , double O , double S , IndexType cindex)
+::ComputeDescriptor(const NeighborhoodType& neigh , double O , double S )
{
-// int y, x, sample_x, sample_y, count=0;
-// int i = 0, ix = 0, j = 0, jx = 0, xs = 0, ys = 0;
-// double dx, dy, mdx, mdy, co, si;
-
-// VectorType descriptorVector;
-// descriptorVector.resize(64);
-
-// double gauss_s1 = 0.f, gauss_s2 = 0.f;
-// double rx = 0.f, ry = 0.f, rrx = 0.f, rry = 0.f, len = 0.f;
-// double cx = -0.5f, cy = 0.f; //Subregion centers for the 4x4 gaussian weighting
-
-
-
-// x = cindex[0];
-// y = cindex[1];
-// co = vcl_cos(O);
-// si = vcl_sin(O);
-
-// i = -8;
-
-// //Calculate descriptorVectorriptor for this interest point
-// //Area of size 24 s x 24 s
-// //***********************************************
-// while(i < 12)
-// {
-// j = -8;
-// i = i-4;
-
-// cx += 1.f;
-// cy = -0.5f;
-
-// while(j < 12)
-// {
-// dx=dy=mdx=mdy=0.f;
-// cy += 1.f;
-
-// j = j - 4;
-
-// ix = i + 5;
-// jx = j + 5;
-
-// xs = vcl_floor(x + ( -jx*sigma*si + ix*sigma*co) +0.5);
-// ys = vcl_floor(y + ( jx*sigma*co + ix*sigma*si) +0.5 );
-
-// for (int k = i; k < i + 9; ++k)
-// {
-// for (int l = j; l < j + 9; ++l)
-// {
-// //Get coords of sample point on the rotated axis
-// sample_x =vcl_floor(x + (-l*sigma*si + k*sigma*co) + 0.5 );
-// sample_y = vcl_floor(y + ( l*sigma*co + k*sigma*si) + 0.5 );
-
-// //Get the gaussian weighted x and y responses
-// gauss_s1= 1./(2.0f*CONST_PI*sigma*sigma)*vcl_exp(-((xs-sample_x)*(xs-sample_x) + (ys-sample_y)*(ys-sample_y) ) / (2*2.5*2.5*sigma*sigma) );
-
-// rx = this->HaarX(sample_y, sample_x, 2*sigma);
-// ry = this->HaarY(sample_y, sample_x, 2*sigma);
-
-// //Get the gaussian weighted x and y responses on rotated axis
-// // rrx = -rx*si + ry*co;
-// // rry = rx*co + ry*si;
-
-// // rrx = gauss_s1*rrx;
-// // rry = gauss_s1*rry;
-
-// rrx = gauss_s1*rx;
-// rry = gauss_s1*ry;
-
-// dx += rrx;
-// dy += rry;
-// mdx += fabs(rrx);
-// mdy += fabs(rry);
-
-// }
-// }
-
-// //Add the values to the descriptorVectorriptor vector
-// //gauss_s2 = 1./(2.0f*CONST_PI*1.5*1.5)*vcl_exp(-((cx-2.)*(cx-2.) + (cy-2.)*(cy-2.) ) / (2*1.5*1.5) );
-// gauss_s2 = 1.;
-// descriptorVector[count++] = dx*gauss_s2;
-// descriptorVector[count++] = dy*gauss_s2;
-// descriptorVector[count++] = mdx*gauss_s2;
-// descriptorVector[count++] = mdy*gauss_s2;
-
-// len += dx*dx + dy*dy + mdx*mdx + mdy*mdy;
-
-// j += 9;
-// }
-// i += 9;
-// }
-
-// //Convert to Unit Vector
-// len = sqrt(len);
-// for(int i = 0; i < 64; i++)
-// {
-// descriptorVector[i] /= len;
-// }
-
-// return descriptorVector;
-
-
-
typedef itk::CenteredRigid2DTransform<double> TransformType;
TransformType::Pointer eulerTransform = TransformType::New();
TransformType::ParametersType ParamVec(5);
PointImageType pSrc , pDst;
+ double angle = O * M_PI / 180;
+
+
+
+ int i = 0, col, raw , Nbin, pas = 1;
+ double xx = 0, yy = 0;
+ double dx, dy , w;
+ int Largeur = 2*neigh.GetRadius()[0]+1;
+ double rayon = static_cast<double>(Largeur)/4.;
+ double r = neigh.GetRadius()[0];
+ double dist = 0;
+ double x0 = neigh.GetCenterNeighborhoodIndex()% Largeur;
+ double y0 = neigh.GetCenterNeighborhoodIndex()/ Largeur;
+
+ //std::cout << " x0 " << x0 << " y0 " << y0 << angle << std::endl;
VectorType descriptorVector;
descriptorVector.resize(64);
/** Parameters of the transformation*/
- ParamVec[0] = O;
- ParamVec[1] = cindex[0];
- ParamVec[2] = cindex[1];
+ ParamVec[0] = angle;
+ ParamVec[1] = x0;
+ ParamVec[2] = y0;
ParamVec[3] = 0;
ParamVec[4] = 0;
eulerTransform->SetParameters(ParamVec);
- int sigma = vcl_floor(S+0.5);
- unsigned int count = 0;
- double cx = -0.5;
- double cy = 0.;
- double si = vcl_sin(O);
- double co = vcl_cos(O);
-
- for(int i = -10 * sigma ; i < 10 * sigma ; i+= 5*sigma )
+ while (i < (int)neigh.Size())
{
- cx += 1.;
- cy = -0.5;
- for(int j = -10 * sigma ; j < 10 * sigma ; j+= 5*sigma )
- {
- cy += 1.;
- double dx = 0.,dy =0., mdx = 0.,mdy = 0.;
-
- for(int k = i ; k< i+5*sigma ; k+= sigma )
- for(int l = j ; l< j+5*sigma ; l+= sigma )
- {
- double distance = vcl_sqrt(static_cast<double>(k*k + i*i));
- if(distance < 10*S)
- {
- typename InputImageType::OffsetType offset ;
- offset[0] = k;
- offset[1] = l;
-
- IndexType currentI = cindex + offset;
- //m_IntegralImage->TransformIndexToPhysicalPoint(currentI,pSrc);
- //pDst = eulerTransform->TransformPoint(pSrc);
- //Get coords of sample point on the rotated axis
- pDst[0] = vcl_floor( (cindex[0] + (-l*si + k*co)) + 0.5 );
- pDst[1] = vcl_floor( (cindex[1] + ( l*co + k*si)) + 0.5 );
-
- double w = 1./(2.*3.3*CONST_PI*S*S)*vcl_exp(-( offset[0]*offset[0]+ offset[1]*offset[1]) / (2.*3.3*3.3*S*S) );
- double rx = this->HaarX(pDst[1],pDst[0],2*sigma);
- double ry = this->HaarY(pDst[1],pDst[0],2*sigma);
-
- //Get the gaussian weighted x and y responses on rotated axis
- double rrx = -rx*vcl_sin(O) + ry*vcl_cos(O);
- double rry = rx*vcl_cos(O) + ry*vcl_sin(O);
- rrx *= w;
- rry *= w;
-
- dx += rrx;
- dy += rry;
- mdx += vcl_abs(rrx);
- mdy += vcl_abs(rry);
- }
- }
-
- //Edit the descriptor
- double w2 = 1./(2.*CONST_PI*1.5*1.5)*vcl_exp(-( (cx -2.)*(cx -2.)+ (cy-2.)*(cy-2.)) / (2.*1.5*1.5 ));
- descriptorVector[count++] = w2*dx;
- descriptorVector[count++] = w2*dy;
- descriptorVector[count++] = w2*mdx;
- descriptorVector[count++] = w2*mdy;
- }
- }
-
-// int i = 0, col, raw , Nbin, pas = 1;
-// double xx = 0, yy = 0;
-// double dx, dy , w;
-// int Largeur = 2*neigh.GetRadius()[0]+1;
-// double rayon = static_cast<double>(Largeur)/4.;
-// double r = neigh.GetRadius()[0];
-// double dist = 0;
-// double x0 = neigh.GetCenterNeighborhoodIndex()% Largeur;
-// double y0 = neigh.GetCenterNeighborhoodIndex()/ Largeur;
-
-
-// VectorType descriptorVector;
-// descriptorVector.resize(64);
-
-// /** Parameters of the transformation*/
-// ParamVec[0] = angle;
-// ParamVec[1] = x0;
-// ParamVec[2] = y0;
-// ParamVec[3] = 0;
-// ParamVec[4] = 0;
-// eulerTransform->SetParameters(ParamVec);
-
-// while (i < (int)neigh.Size())
-// {
-// col = i % Largeur;
-// raw = i / Largeur;
-
-// if (( col > pas && col < Largeur - pas ) && ( raw > pas && raw < Largeur - pas) )
-// {
-// double distanceX = (raw-r);
-// double distanceY = (col-r);
-// dist = vcl_sqrt(distanceX*distanceX + distanceY*distanceY);
-
-// if (dist <= r )
-// {
-// /* Transform point to compensate the rotation the orientation */
-// pDst[0] = col;
-// pDst[1] = raw;
-// pSrc = eulerTransform->TransformPoint(pDst);
-
-// /** New Coordinates (rotated) */
-// col = static_cast<int>(vcl_floor(pSrc[0]));
-// raw = static_cast<int>(vcl_floor(pSrc[1]));
-
-// if (raw==0) raw =+1;
-// if (col ==0) col +=1;
-
-// xx = static_cast<int> (pSrc[1]/rayon);
-// yy = static_cast<int> (pSrc[0]/rayon);
-// Nbin = static_cast<int> (xx + 4*yy);
-
-// if ( Nbin < 16) //because 64 descriptor length
-// {
-// double distanceXcompensee_2 = (pSrc[0] - r)*(pSrc[0] - r);
-// double distanceYcompensee_2 = (pSrc[1] - r)*(pSrc[1] - r);
-
-// w = 1./(2.*CONST_PI*S*S)*vcl_exp(-( distanceXcompensee_2 + distanceYcompensee_2 ) / (2*3.3*3.3*S*S) );
-
-// dx = 0.5 * (neigh[(col+pas) + raw * Largeur] - neigh[(col-pas) + raw *Largeur]) * w;
-// dy = 0.5 * (neigh[col + (raw+ pas)* Largeur] - neigh[col + (raw-pas)*Largeur]) * w;
-
-// //dx = this->HaarX(pSrc[1] ,pSrc[0] , 2*vcl_floor(S+0.5));
-// //dy = this->HaarY(pSrc[1] ,pSrc[0] , 2*vcl_floor(S+0.5));
-
-// descriptorVector[4*Nbin ] += dx;
-// descriptorVector[4*Nbin+1] += dy;
-// descriptorVector[4*Nbin+2] += vcl_abs(dx);
-// descriptorVector[4*Nbin+3] += vcl_abs(dy);
-// }
-// }
-// }
-// ++i;
-// }
-
- double accu = 0;
- for (int i = 0; i < 64; ++i)
- accu += descriptorVector[i]*descriptorVector[i];
-
- for (int j = 0; j < 64; ++j)
- descriptorVector[j] /= vcl_sqrt(accu);
+ col = i % Largeur;
+ raw = i / Largeur;
- return descriptorVector;
-}
+ if (( col > pas && col < Largeur - pas ) && ( raw > pas && raw < Largeur - pas) )
+ {
+ double distanceX = (raw-r);
+ double distanceY = (col-r);
+ dist = vcl_sqrt(distanceX*distanceX + distanceY*distanceY);
-/*----------------------------------------------------------------
- Compute Integral Images
- -----------------------------------------------------------------*/
-template <class TInputImage, class TOutputPointSet >
-void
-ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
-::IntegralImage()
-{
- //InputImagePointerType m_IntegralImage = InputImageType::New();
- InputImagePointerType sourceImage = const_cast<InputImageType*>(this->GetInput());
+ if (dist <= r )
+ {
+ /* Transform point to compensate the rotation the orientation */
+ pDst[0] = col;
+ pDst[1] = raw;
+ pSrc = eulerTransform->TransformPoint(pDst);
- m_IntegralImage->SetRegions(sourceImage->GetLargestPossibleRegion());
- m_IntegralImage->Allocate();
- m_IntegralImage->FillBuffer(0.);
+ /** New Coordinates (rotated) */
+ col = static_cast<int>(vcl_floor(pSrc[0]));
+ raw = static_cast<int>(vcl_floor(pSrc[1]));
- itk::ImageRegionIterator<InputImageType> itSource(sourceImage,
- sourceImage->GetLargestPossibleRegion());
+ if (raw==0) raw =+1;
+ if (col ==0) col +=1;
- itk::ImageRegionIterator<InputImageType> itInt(m_IntegralImage,
- m_IntegralImage->GetLargestPossibleRegion());
+ xx = static_cast<int> (pSrc[1]/rayon);
+ yy = static_cast<int> (pSrc[0]/rayon);
+ Nbin = static_cast<int> (xx + 4*yy);
- itSource.GoToBegin();
- itInt.GoToBegin();
- double rowSum = 0.;
+ if ( Nbin < 16) //because 64 descriptor length
+ {
+ double distanceXcompensee_2 = (pSrc[0] - r)*(pSrc[0] - r);
+ double distanceYcompensee_2 = (pSrc[1] - r)*(pSrc[1] - r);
- while(!itSource.IsAtEnd())
- {
- IndexType index = itSource.GetIndex();
+ w = vcl_exp(-( distanceXcompensee_2 + distanceYcompensee_2 ) / (2*3.3*3.3*S*S) );
- if(index[0] == 0)
- rowSum = 0.;
+ dx = 0.5 * (neigh[(col+pas) + raw * Largeur] - neigh[(col-pas) + raw *Largeur]) * w;
+ dy = 0.5 * (neigh[col + (raw+ pas)* Largeur] - neigh[col + (raw-pas)*Largeur]) * w;
- if(index[1] == 0 )
- {
- rowSum += itSource.Get();
- itInt.Set(rowSum);
- }
- else
- {
- //Index of the upper pixel (raw -1 & col)
- IndexType prevIndex = index;
- prevIndex[1] -= 1;
- itInt.SetIndex(prevIndex);
- double prevValue = itInt.Get();
-
- //Edit the integral image
- itInt.SetIndex(index);
- rowSum += itSource.Get();
- itInt.Set(prevValue + rowSum);
- }
- ++itSource;
- ++itInt;
+ descriptorVector[4*Nbin ] += dx;
+ descriptorVector[4*Nbin+1] += dy;
+ descriptorVector[4*Nbin+2] += vcl_abs(dx);
+ descriptorVector[4*Nbin+3] += vcl_abs(dy);
+ }
+ }
}
+ i++;
+ }
-}
-
-/*----------------------------------------------------------------
- * GetAngle
- * -----------------------------------------------------------------*/
-template <class TInputImage, class TOutputPointSet >
-double
-ImageToSURFKeyPointSetFilter< TInputImage, TOutputPointSet>
-::GetAngle(unsigned int X, unsigned int Y)
-{
- if(X > 0 && Y >= 0)
- return vcl_atan2(Y,X);
-
- if(X < 0 && Y >= 0)
- return CONST_PI - vcl_atan2(-Y,X);
-
- if(X < 0 && Y < 0)
- return CONST_PI + vcl_atan2(Y,X);
+ double accu = 0;
+ for (int i = 0; i < 64; i++)
+ accu += descriptorVector[i]*descriptorVector[i];
- if(X > 0 && Y < 0)
- return 2*CONST_PI - vcl_atan2(-Y,X);
+ for (int j = 0; j < 64; j++)
+ descriptorVector[j] /= vcl_sqrt(accu);
+ return descriptorVector;
- return 0;
}
/*----------------------------------------------------------------
diff --git a/Code/FeatureExtraction/otbLineSegmentDetector.txx b/Code/FeatureExtraction/otbLineSegmentDetector.txx
index e35d4e44953fa4353cf7f0002491a8dcab38593a..1d222944210b3e1833a968509bb74c223caa2f21 100644
--- a/Code/FeatureExtraction/otbLineSegmentDetector.txx
+++ b/Code/FeatureExtraction/otbLineSegmentDetector.txx
@@ -851,7 +851,7 @@ LineSegmentDetector<TInputImage, TPrecision>
/** Get The Bounding Region*/
OutputImageDirRegionType region = rectangle->GetBoundingRegion();
- region.Crop(m_OrientationFilter->GetOutput()->GetRequestedRegion());
+ region.Crop( m_OrientationFilter->GetOutput()->GetLargestPossibleRegion() );
itk::ImageRegionIterator<OutputImageDirType> it(m_OrientationFilter->GetOutput(), region/*m_OrientationFilter->GetOutput()->GetRequestedRegion()*/);
it.GoToBegin();
@@ -862,7 +862,7 @@ LineSegmentDetector<TInputImage, TPrecision>
if( rectangle->IsInside( it.GetIndex()) && m_OrientationFilter->GetOutput()->GetRequestedRegion().IsInside( it.GetIndex()) )
{
++pts;
-
+
if(this->IsAligned(it.Get(), rec[5] /*theta*/ ,rec[6] /*Prec*/))
NbAligned++;
}
diff --git a/Testing/Code/FeatureExtraction/otbLineSegmentDetector.cxx b/Testing/Code/FeatureExtraction/otbLineSegmentDetector.cxx
index 32f2aace426bb3f2f65edb1aa11b726423754a31..74ebc650609558bfb16d409a1d74b2592a0693fc 100644
--- a/Testing/Code/FeatureExtraction/otbLineSegmentDetector.cxx
+++ b/Testing/Code/FeatureExtraction/otbLineSegmentDetector.cxx
@@ -62,12 +62,10 @@ int otbLineSegmentDetector( int argc, char * argv[] )
reader->SetFileName(infname);
reader->GenerateOutputInformation();
- reader->Update();
-
-// lsdFilter->SetInput(reader->GetOutput());
+ lsdFilter->SetInput(reader->GetOutput());
-// drawLineFilter->SetInput(reader->GetOutput());
-// drawLineFilter->SetInputLineSpatialObjectList(lsdFilter->GetOutput());
+ drawLineFilter->SetInput(reader->GetOutput());
+ drawLineFilter->SetInputLineSpatialObjectList(lsdFilter->GetOutput());
// LinesListType::const_iterator it = lsdFilter->GetOutput()->begin();
@@ -101,13 +99,13 @@ int otbLineSegmentDetector( int argc, char * argv[] )
-// /** Write The Output Image*/
-// WriterType::Pointer writer = WriterType::New();
-// writer->SetFileName(outfname);
-// writer->SetInput(drawLineFilter->GetOutput());
-// writer->Update();
+ /** Write The Output Image*/
+ WriterType::Pointer writer = WriterType::New();
+ writer->SetFileName(outfname);
+ writer->SetInput(drawLineFilter->GetOutput());
+ writer->Update();
-// std::cout << " lsdFilter Output Size" << lsdFilter->GetOutput()->size() <<std::endl;
+ std::cout << " lsdFilter Output Size" << lsdFilter->GetOutput()->size() <<std::endl;
return EXIT_SUCCESS;
diff --git a/Utilities/otbsiftfast/libsiftfast.cpp b/Utilities/otbsiftfast/libsiftfast.cpp
index b0a7ce68aee65a8a523068e6fb3955a2989d4598..baf2d4b01f3292703bbd2478da9b26c4a9cf6ba8 100755
--- a/Utilities/otbsiftfast/libsiftfast.cpp
+++ b/Utilities/otbsiftfast/libsiftfast.cpp
@@ -1616,41 +1616,40 @@ void PlaceInIndex(float* fdesc, float mag, float ori, float rx, float cx)
neworient = (int)oribin;
ofrac = oribin-(float)neworient;
- if( newrow >= -1 && newrow < 4 && neworient >= 0 && neworient <= 8 && rfrac >= 0 && rfrac < 1)
- {
- for(int i = 0; i < 2; ++i) {
- if( (unsigned int)(i+newrow) >= 4 )
+ assert( newrow >= -1 && newrow < 4 && neworient >= 0 && neworient <= 8 && rfrac >= 0 && rfrac < 1);
+
+ for(int i = 0; i < 2; ++i) {
+ if( (unsigned int)(i+newrow) >= 4 )
continue;
- float frowgrad;
- if( i == 0 )
+ float frowgrad;
+ if( i == 0 )
frowgrad = mag*(1-rfrac);
- else
+ else
frowgrad = mag*rfrac;
- for(int j = 0; j < 2; ++j) {
+ for(int j = 0; j < 2; ++j) {
if( (unsigned int)(j+newcol) >= 4 )
- continue;
+ continue;
float fcolgrad;
if( j == 0 )
- fcolgrad = frowgrad*(1-cfrac);
+ fcolgrad = frowgrad*(1-cfrac);
else
- fcolgrad = frowgrad*cfrac;
+ fcolgrad = frowgrad*cfrac;
float* pfdescorient = fdesc + 8*(4*(i+newrow)+j+newcol);
for(int k = 0; k < 2; ++k) {
- float forigrad;
- if( k == 0 )
- forigrad = fcolgrad*(1-ofrac);
- else
- forigrad = fcolgrad*ofrac;
+ float forigrad;
+ if( k == 0 )
+ forigrad = fcolgrad*(1-ofrac);
+ else
+ forigrad = fcolgrad*ofrac;
- pfdescorient[(neworient+k)&7] += forigrad;
+ pfdescorient[(neworient+k)&7] += forigrad;
}
- }
- }
- }
+ }
+ }
}
void FreeKeypoints(Keypoint keypt)