ITK/Examples/ImageSegmentation/ExtractContourWithSnakes: Difference between revisions

From KitwarePublic
< ITK‎ | Examples
Jump to navigationJump to search
(Deprecated content that is moved to sphinx)
 
(7 intermediate revisions by 3 users not shown)
Line 1: Line 1:
==ExtractContourWithSnakes.cxx==
{{warning|1=The media wiki content on this page is no longer maintained. The examples presented on the https://itk.org/Wiki/* pages likely require ITK version 4.13 or earlier releases.   In many cases, the examples on this page no longer conform to the best practices for modern ITK versions.}}
<source lang="cpp">
#include "itkImage.h"
#include "itkRandomImageSource.h"
#include "itkGradientRecursiveGaussianImageFilter.h"
#include "itkGradientMagnitudeImageFilter.h"
#include "itkImageFileReader.h"
 
#include <vnl/vnl_matrix.h>
#include "vnl/algo/vnl_determinant.h"
#include "vnl/algo/vnl_matrix_inverse.h"
#include <vnl/vnl_vector.h>
 
#include <iostream>
 
typedef itk::Image< unsigned char, 2 >  ImageType;
typedef itk::Image< float,  2 >        FloatImageType;
typedef ImageType::IndexType IndexType;
typedef itk::CovariantVector< float, 2  >      OutputPixelType;
typedef itk::Image< OutputPixelType, 2 >    OutputImageType;
typedef itk::GradientRecursiveGaussianImageFilter<FloatImageType, OutputImageType> FilterType;
typedef itk::GradientMagnitudeImageFilter<ImageType, FloatImageType >  gradMagfilterType;
typedef itk::ImageFileReader< ImageType >  readerType;
 
vnl_vector<double> generateCircle(double cx, double cy, double rx, double ry, int n);
void createImage(ImageType::Pointer image, int w, int h, double cx, double cy, double rx, double ry);
vnl_matrix<double> computeP(double alpha, double beta, double gamma, double N) throw (int);
vnl_vector<double> sampleImage(vnl_vector<double> x, vnl_vector<double> y, OutputImageType::Pointer gradient, int position);
 
 
int main( int argc, char* argv[] )
{
//Image dimensions
int w = 300;
int h = 300;
ImageType::Pointer image;
if (argc < 7)
{
std::cout << "Usage " << argv[0] << " points alpha beta gamma sigma iterations [image]" << std::endl;
return 0;
} else if (argc < 8) {
//Synthesize the image
image = ImageType::New();
createImage(image, w, h, 150, 150, 50, 50);
} else if (argc == 8) {
//Open the image
readerType::Pointer reader = readerType::New();
reader->SetFileName( argv[7] );
try {
reader->Update();
image = reader->GetOutput();
w = image->GetLargestPossibleRegion().GetSize()[0];
h = image->GetLargestPossibleRegion().GetSize()[1];
} catch( itk::ExceptionObject & err ) {
return 3;
}
}
 
//Snake parameters
double alpha = 0.001;
double beta = 0.4;
double gamma = 100;
double iterations = 1;
int nPoints = 20;
double sigma;
 
nPoints = atoi(argv[1]);
alpha = atof(argv[2]);
beta = atof(argv[3]);
gamma = atof(argv[4]);
sigma = atof(argv[5]);
iterations = atoi(argv[6]);
 
//Temporal variables
vnl_matrix<double> P;
vnl_vector<double> v;
double N;
 
//Generate initial snake circle
v = generateCircle(130, 130, 50, 50, nPoints);
 
//Computes P matrix.
N = v.size()/2;
try {
P = computeP(alpha, beta, gamma, N);
} catch (int n) {
return n;
}
 
//Computes the magnitude gradient
gradMagfilterType::Pointer gradientMagnitudeFilter = gradMagfilterType::New();
gradientMagnitudeFilter->SetInput( image );
gradientMagnitudeFilter->Update();
 
//Computes the gradient of the gradient magnitude
FilterType::Pointer gradientFilter = FilterType::New();
gradientFilter->SetInput( gradientMagnitudeFilter->GetOutput() );
gradientFilter->SetSigma( sigma );
gradientFilter->Update();
 
//Loop
vnl_vector<double> x(N);
vnl_vector<double> y(N);
 
std::cout << "Initial snake" << std::endl;
for (int i=0; i<N; i++) {
x[i] = v[2*i];
y[i] = v[2*i+1];
std::cout << "(" << x[i] << ", " << y[i] << ")" << std::endl;
}
 
for (int i=0; i<iterations; i++) {
vnl_vector<double> fex;
vnl_vector<double> fey;
fex = sampleImage(x, y, gradientFilter->GetOutput(), 0);
fey = sampleImage(x, y, gradientFilter->GetOutput(), 1);
 
x = (x+gamma*fex).post_multiply(P);
y = (y+gamma*fey).post_multiply(P);
}
 
//Display the answer
std::cout << "Final snake after " << iterations << " iterations" << std::endl;
vnl_vector<double> v2(2*N);
for (int i=0; i<N; i++) {
v2[2*i] = x[i];
v2[2*i+1] = y[i];
std::cout << "(" << x[i] << ", " << y[i] << ")" << std::endl;
}
 
return 0;
}
 
vnl_vector<double> generateCircle(double cx, double cy, double rx, double ry, int n) {
vnl_vector<double> v(2*(n+1));
 
for (int i=0; i<n; i++) {
v[2*i] = cx + rx*cos(2*M_PI*i/n);
v[2*i+1] = cy + ry*sin(2*M_PI*i/n);
}
v[2*n]=v[0];
v[2*n+1]=v[1];
return v;
}
 
void createImage(ImageType::Pointer image, int w, int h, double cx, double cy, double rx, double ry) {
 
itk::Size<2> size;
size[0] = w;
size[1] = h;
 
itk::RandomImageSource<ImageType>::Pointer randomImageSource = itk::RandomImageSource<ImageType>::New();
randomImageSource->SetNumberOfThreads(1); // to produce non-random results
randomImageSource->SetSize(size);
randomImageSource->SetMin(200);
randomImageSource->SetMax(255);
randomImageSource->Update();
 
image->SetRegions(randomImageSource->GetOutput()->GetLargestPossibleRegion());
image->Allocate();
 
IndexType index;
 
//Draw oval.
for (int i=0; i<w; i++) {
for (int j=0; j<h; j++) {
index[0] = i; index[1] = j;
if ( ((i-cx)*(i-cx)/(rx*rx) + (j-cy)*(j-cy)/(ry*ry) ) < 1) {
image->SetPixel(index, randomImageSource->GetOutput()->GetPixel(index)-100);
} else {
image->SetPixel(index, randomImageSource->GetOutput()->GetPixel(index));
}
 
}
}
}
 
vnl_matrix<double> computeP(double alpha, double beta, double gamma, double N) throw (int) {
 
double a = gamma*(2*alpha+6*beta)+1;
double b = gamma*(-alpha-4*beta);
double c = gamma*beta;
 
vnl_matrix<double> P(N,N);
 
P.fill(0);
 
//fill diagonal
P.fill_diagonal(a);
 
//fill next two diagonals
for (int i=0; i<(N-1); i++) {
P(i+1,i) = b;
P(i,i+1) = b;
}
//Moreover
P(0, N-1)=b;
P(N-1, 0)=b;
 
//fill next two diagonals
for (int i=0; i<(N-2); i++) {
P(i+2,i) = c;
P(i,i+2) = c;
}
//Moreover
P(0, N-2)=c;
P(1, N-1)=c;
P(N-2, 0)=c;
P(N-1, 1)=c;
 
if ( vnl_determinant(P) == 0.0 ) {
std::cerr << "Singular matrix. Determinant is 0." << std::endl;
throw 2;
}
 
//Compute the inverse of the matrix P
vnl_matrix< double > Pinv;
Pinv = vnl_matrix_inverse< double >(P);
 
return Pinv;
}
 
vnl_vector<double> sampleImage(vnl_vector<double> x, vnl_vector<double> y, OutputImageType::Pointer gradient, int position) {
int size;
size = x.size();
vnl_vector<double> ans(size);
 
IndexType index;
for (int i=0; i<size; i++) {
index[0] = x[i];
index[1] = y[i];
ans[i] = gradient->GetPixel(index)[position];
}
return ans;
}
</source>
 
{{ITKCMakeLists|{{SUBPAGENAME}}}}

Latest revision as of 16:32, 7 June 2019

Warning: The media wiki content on this page is no longer maintained. The examples presented on the https://itk.org/Wiki/* pages likely require ITK version 4.13 or earlier releases. In many cases, the examples on this page no longer conform to the best practices for modern ITK versions.