Skip to content
Snippets Groups Projects
Commit 1cd74a22 authored by Patrick Imbo's avatar Patrick Imbo
Browse files

MAJ -> OK

parent 58444d9f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Code/FeatureExtraction/otbImageToPathAlignFilter.h
+ 27
− 13
View file @ 1cd74a22
......@@ -14,11 +14,12 @@
#include "itkImageSource.h"
#include "itkConceptChecking.h"
#include "itkImage.h"
namespace otb
{
/** \class PathToImageAlignFilter
/** \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
......@@ -27,11 +28,11 @@ namespace otb
* assumed internally to be 1.0).
*/
template <class TInputImage, class TOutputPath>
class PathToImageAlignFilter : public itk::ImageSource<TInputImage>
class ImageToPathAlignFilter : public itk::ImageSource<TInputImage>
{
public:
/** Standard class typedefs. */
typedef PathToImageAlignFilter Self;
typedef ImageToPathAlignFilter Self;
typedef itk::ImageSource<TInputImage> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
......@@ -40,22 +41,30 @@ public:
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(PathToImageAlignFilter,ImageSource);
itkTypeMacro(ImageToPathAlignFilter,ImageSource);
/** ImageDimension constants */
itkStaticConstMacro(InputImageDimension, unsigned int,
TInputImage::ImageDimension);
/** Some convenient typedefs. */
typedef typename Superclass::InputImageRegionType InputImageRegionType;
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;
/** ImageDimension constants */
itkStaticConstMacro(InputImageDimension, unsigned int,
TInputImage::ImageDimension);
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
......@@ -84,8 +93,8 @@ public:
itkGetMacro(Size,SizeType);
protected:
PathToImageAlignFilter();
~PathToImageAlignFilter();
ImageToPathAlignFilter();
~ImageToPathAlignFilter();
virtual void GenerateOutputInformation(){}; // do nothing
virtual void GenerateData();
......@@ -100,8 +109,13 @@ protected:
virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
private:
PathToImageAlignFilter(const Self&); //purposely not implemented
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
};
......
Code/FeatureExtraction/otbImageToPathAlignFilter.txx
+ 293
− 23
View file @ 1cd74a22
......@@ -16,14 +16,27 @@
#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>
PathToImageAlignFilter<TInputImage,TOutputPath>
::PathToImageAlignFilter()
ImageToPathAlignFilter<TInputImage,TOutputPath>
::ImageToPathAlignFilter()
{
this->SetNumberOfRequiredInputs(1);
m_Size.Fill(0);
......@@ -41,8 +54,8 @@ PathToImageAlignFilter<TInputImage,TOutputPath>
/** Destructor */
template <class TInputImage, class TOutputPath>
PathToImageAlignFilter<TInputImage,TOutputPath>
::~PathToImageAlignFilter()
ImageToPathAlignFilter<TInputImage,TOutputPath>
::~ImageToPathAlignFilter()
{
}
......@@ -50,7 +63,7 @@ PathToImageAlignFilter<TInputImage,TOutputPath>
//----------------------------------------------------------------------------
template <class TInputImage, class TOutputPath>
void
PathToImageAlignFilter<TInputImage,TOutputPath>
ImageToPathAlignFilter<TInputImage,TOutputPath>
::SetSpacing(const double* spacing)
{
unsigned int i;
......@@ -72,7 +85,7 @@ PathToImageAlignFilter<TInputImage,TOutputPath>
template <class TInputImage, class TOutputPath>
void
PathToImageAlignFilter<TInputImage,TOutputPath>
ImageToPathAlignFilter<TInputImage,TOutputPath>
::SetSpacing(const float* spacing)
{
unsigned int i;
......@@ -94,7 +107,7 @@ PathToImageAlignFilter<TInputImage,TOutputPath>
template <class TInputImage, class TOutputPath>
const double *
PathToImageAlignFilter<TInputImage,TOutputPath>
ImageToPathAlignFilter<TInputImage,TOutputPath>
::GetSpacing() const
{
return m_Spacing;
......@@ -103,7 +116,7 @@ PathToImageAlignFilter<TInputImage,TOutputPath>
//----------------------------------------------------------------------------
template <class TInputImage, class TOutputPath>
void
PathToImageAlignFilter<TInputImage,TOutputPath>
ImageToPathAlignFilter<TInputImage,TOutputPath>
::SetOrigin(const double* origin)
{
unsigned int i;
......@@ -125,7 +138,7 @@ PathToImageAlignFilter<TInputImage,TOutputPath>
template <class TInputImage, class TOutputPath>
void
PathToImageAlignFilter<TInputImage,TOutputPath>
ImageToPathAlignFilter<TInputImage,TOutputPath>
::SetOrigin(const float* origin)
{
unsigned int i;
......@@ -147,7 +160,7 @@ PathToImageAlignFilter<TInputImage,TOutputPath>
template <class TInputImage, class TOutputPath>
const double *
PathToImageAlignFilter<TInputImage,TOutputPath>
ImageToPathAlignFilter<TInputImage,TOutputPath>
::GetOrigin() const
{
return m_Origin;
......@@ -157,10 +170,12 @@ PathToImageAlignFilter<TInputImage,TOutputPath>
template <class TInputImage, class TOutputPath>
void
PathToImageAlignFilter<TInputImage,TOutputPath>
ImageToPathAlignFilter<TInputImage,TOutputPath>
::GenerateData(void)
{
unsigned int i;
int i;
itkDebugMacro(<< "ImageToPathAlignFilter::GenerateData() called");
// Get the input and output pointers
......@@ -211,7 +226,7 @@ PathToImageAlignFilter<TInputImage,TOutputPath>
InputImage->SetLargestPossibleRegion( region); //
InputImage->SetBufferedRegion( region ); // set the region
InputImage->SetRequestedRegion( 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
......@@ -239,17 +254,272 @@ PathToImageAlignFilter<TInputImage,TOutputPath>
InputImage->SetOrigin(origin); // and origin
InputImage->Allocate(); // allocate the image
itk::ImageRegionIteratorWithIndex<InputImageType> imageIt(InputImage,region);
for( imageIt.GoToBegin(); !imageIt.IsAtEnd(); ++imageIt )
{
imageIt.Set(m_BackgroundValue);
}
RealImagePointer AngleImage = RealImageType::New();
AngleImage->SetRegions( InputImage->GetRequestedRegion() );
AngleImage->CopyInformation( InputImage );
AngleImage->Allocate();
typedef itk::ImageLinearIteratorWithIndex< InputImageType > IteratorType;
typedef itk::ImageLinearConstIteratorWithIndex< RealType > ConstIteratorType;
itk::PathIterator<InputImageType,OutputPathType> pathIt(InputImage,OutputPath);
for( pathIt.GoToBegin(); !pathIt.IsAtEnd(); ++pathIt )
{
pathIt.Set( m_PathValue );
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");
......@@ -258,7 +528,7 @@ PathToImageAlignFilter<TInputImage,TOutputPath>
template <class TInputImage, class TOutputPath>
void
PathToImageAlignFilter<TInputImage,TOutputPath>
ImageToPathAlignFilter<TInputImage,TOutputPath>
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment