[Insight-users] Question on DeformableModel1.cxx
John Drozd
john.drozd at gmail.com
Tue Jun 29 12:13:23 EDT 2010
Hello,
I have a binary mask dicom image with orientation 0 0 -1 / 0 1 0
I ran it through the DeformableModel1.cxx program,
to create and deform a mesh over the mask,
but my outputted deformed dicom image has the orientation 1 0 0 / 0 1 0
I tried copying the MetaDataDictionary from the mask to the mesh,
but this did not help.
I tried setting the direction of the outputted dicom image,
but when viewed in 3D Slicer, there was no image (all black).
I suppose the points on the mesh would have to be transformed to fit the 0 0
-1 / 0 1 0 orientation.
Can anyone help?
Below is my code: ( a slightly modified version of DeformableModel1.cxx)
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: DeformableModel1.cxx,v $
Language: C++
Date: $Date: 2009-03-17 21:44:42 $
Version: $Revision: 1.30 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#ifdef __BORLANDC__
#define ITK_LEAN_AND_MEAN
#endif
// Software Guide : BeginLatex
//
// This example illustrates the use of the \doxygen{DeformableMesh3DFilter}
// and \doxygen{BinaryMask3DMeshSource} in the hybrid segmentation
framework.
//
// \begin{figure} \center
//
\includegraphics[width=\textwidth]{DeformableModelCollaborationDiagram.eps}
// \itkcaption[Deformable model collaboration diagram]{Collaboration
// diagram for the DeformableMesh3DFilter applied to a segmentation task.}
// \label{fig:DeformableModelCollaborationDiagram}
// \end{figure}
//
// The purpose of the DeformableMesh3DFilter is to take an initial surface
// described by an \doxygen{Mesh} and deform it in order to adapt it to the
// shape of an anatomical structure in an image.
// Figure~\ref{fig:DeformableModelCollaborationDiagram} illustrates a
typical
// setup for a segmentation method based on deformable models. First, an
// initial mesh is generated using a binary mask and an isocontouring
// algorithm (such as marching cubes) to produce an initial mesh. The binary
// mask used here contains a simple shape which vaguely resembles the
// anatomical structure that we want to segment. The application of the
// isocontouring algorithm produces a $3D$ mesh that has the shape of this
// initial structure. This initial mesh is passed as input to the deformable
// model which will apply forces to the mesh points in order to reshape the
// surface until make it fit to the anatomical structures in the image.
//
// The forces to be applied on the surface are computed from an approximate
// physical model that simulates an elastic deformation. Among the forces to
// be applied we need one that will pull the surface to the position of the
// edges in the anatomical structure. This force component is represented
// here in the form of a vector field and is computed as illustrated in the
// lower left of Figure~\ref{fig:DeformableModelCollaborationDiagram}. The
// input image is passed to a
// \doxygen{GradientMagnitudeRecursiveGaussianImageFilter}, which computes
// the magnitude of the image gradient. This scalar image is then passed to
// another gradient filter
// (\doxygen{GradientRecursiveGaussianImageFilter}). The output of this
// second gradient filter is a vector field in which every vector points to
// the closest edge in the image and has a magnitude proportional to the
// second derivative of the image intensity along the direction of the
// gradient. Since this vector field is computed using Gaussian derivatives,
// it is possible to regulate the smoothness of the vector field by playing
// with the value of sigma assigned to the Gaussian. Large values of sigma
// will result in a large capture radius, but will have poor precision in
the
// location of the edges. A reasonable strategy may involve the use of large
// sigmas for the initial iterations of the model and small sigmas to refine
// the model when it is close to the edges. A similar effect could be
// achieved using multiresolution and taking advantage of the image pyramid
// structures already illustrated in the registration framework.
//
// \index{Deformable Models}
// \index{DeformableMesh3DFilter}
//
// Software Guide : EndLatex
#include <iostream>
// Software Guide : BeginLatex
//
// We start by including the headers of the main classes required for this
// example. The BinaryMask3DMeshSource is used to produce an initial
// approximation of the shape to be segmented. This filter takes a binary
// image as input and produces a Mesh as output using the marching cube
// isocontouring algorithm.
//
// \index{itk::BinaryMask3DMeshSource!Header}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkBinaryMask3DMeshSource.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Then we include the header of the DeformableMesh3DFilter that
// implements the deformable model algorithm.
//
// \index{itk::DeformableMesh3DFilter!Header}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkDeformableMesh3DFilter.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We also need the headers of the gradient filters that will be used for
// computing the vector field. In our case they are the
// GradientMagnitudeRecursiveGaussianImageFilter and
// GradientRecursiveGaussianImageFilter.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkGradientRecursiveGaussianImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The main data structures required in this example are the Image
// and the Mesh classes. The deformable model \emph{per se} is
// represented as a Mesh.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkImage.h"
#include "itkMesh.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The \code{PixelType} of the image derivatives is represented with a
// \doxygen{CovariantVector}. We include its header in the following line.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkCovariantVector.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The deformed mesh is converted into a binary image using the
// \doxygen{PointSetToImageFilter}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkPointSetToImageFilter.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// In order to read both the input image and the mask image, we need the
// \doxygen{ImageFileReader} class. We also need the
\doxygen{ImageFileWriter}
// to save the resulting deformed mask image.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
// Software Guide : EndCodeSnippet
int main( int argc, char *argv[] )
{
if( argc < 4 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " InputImage BinaryImage DeformedMaskImage" << std::endl;
return 1;
}
// Software Guide : BeginLatex
//
// Here we declare the type of the image to be processed. This implies a
// decision about the \code{PixelType} and the dimension. The
// DeformableMesh3DFilter is specialized for $3D$, so the choice
// is clear in our case.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
const unsigned int Dimension = 3;
//typedef double PixelType;
typedef signed short PixelType;
typedef itk::Image<PixelType, Dimension> ImageType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The input to BinaryMask3DMeshSource is a binary mask that we
// will read from a file. This mask could be the result of a rough
// segmentation algorithm applied previously to the same anatomical
// structure. We declare below the type of the binary mask image.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
//typedef itk::Image< unsigned char, Dimension > BinaryImageType;
typedef itk::Image< signed short, Dimension > BinaryImageType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Then we define the type of the deformable mesh. We represent the
// deformable model using the Mesh class. The \code{double} type used as
// template parameter here is to be used to assign values to every point
// of the Mesh.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
//typedef itk::Mesh<double> MeshType;
typedef itk::Mesh<signed short> MeshType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The following lines declare the type of the gradient image:
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
//typedef itk::CovariantVector< double, Dimension > GradientPixelType;
typedef itk::CovariantVector< signed short, Dimension >
GradientPixelType;
typedef itk::Image< GradientPixelType, Dimension > GradientImageType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// With it we can declare the type of the gradient filter and the
gradient
// magnitude filter:
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::GradientRecursiveGaussianImageFilter<ImageType,
GradientImageType>
GradientFilterType;
typedef
itk::GradientMagnitudeRecursiveGaussianImageFilter<ImageType,ImageType>
GradientMagnitudeFilterType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The filter implementing the isocontouring algorithm is the
// BinaryMask3DMeshSource filter.
//
// \index{itk::BinaryMask3DMeshSource!Instantiation}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::BinaryMask3DMeshSource< BinaryImageType, MeshType >
MeshSourceType;
// Software Guide : EndCodeSnippet
// typedef itk::BinaryMaskToNarrowBandPointSetFilter<
// BinaryImageType, MeshType > MeshSourceType;
// Software Guide : BeginLatex
//
// Now we instantiate the type of the DeformableMesh3DFilter that
// implements the deformable model algorithm. Note that both the input
and
// output types of this filter are \doxygen{Mesh} classes.
//
// \index{DeformableMesh3DFilter!Instantiation}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::DeformableMesh3DFilter<MeshType,MeshType>
DeformableFilterType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Let's declare two readers. The first will read the image to be
// segmented. The second will read the binary mask containing a first
// approximation of the segmentation that will be used to initialize a
// mesh for the deformable model.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::ImageFileReader< ImageType > ReaderType;
typedef itk::ImageFileReader< BinaryImageType > BinaryReaderType;
ReaderType::Pointer imageReader = ReaderType::New();
BinaryReaderType::Pointer maskReader = BinaryReaderType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// In this example we take the filenames of the input image and the
binary
// mask from the command line arguments.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
imageReader->SetFileName( argv[1] );
maskReader->SetFileName( argv[2] );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We create here the GradientMagnitudeRecursiveGaussianImageFilter that
// will be used to compute the magnitude of the input image gradient. As
// usual, we invoke its \code{New()} method and assign the result to a
// \doxygen{SmartPointer}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
GradientMagnitudeFilterType::Pointer gradientMagnitudeFilter
=
GradientMagnitudeFilterType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The output of the image reader is connected as input to the gradient
// magnitude filter. Then the value of sigma used to blur the image is
// selected using the method \code{SetSigma()}.
//
//
\index{itk::Gradient\-Magnitude\-Recursive\-Gaussian\-Image\-Filter!SetInput()}
//
\index{itk::Gradient\-Magnitude\-Recursive\-Gaussian\-Image\-Filter!SetSigma()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
gradientMagnitudeFilter->SetInput( imageReader->GetOutput() );
gradientMagnitudeFilter->SetSigma( 1.0 );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// In the following line, we construct the gradient filter that will take
// the gradient magnitude of the input image that will be passed to the
// deformable model algorithm.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
GradientFilterType::Pointer gradientMapFilter = GradientFilterType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The magnitude of the gradient is now passed to the next step of
// gradient computation. This allows us to obtain a second derivative of
// the initial image with the gradient vector pointing to the maxima of
// the input image gradient. This gradient map will have the properties
// desirable for attracting the deformable model to the edges of the
// anatomical structure on the image. Once again we must select the value
// of sigma to be used in the blurring process.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
gradientMapFilter->SetInput( gradientMagnitudeFilter->GetOutput());
gradientMapFilter->SetSigma( 1.0 );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// At this point, we are ready to compute the vector field. This is done
// simply by invoking the \code{Update()} method on the second derivative
// filter. This was illustrated in
// Figure~\ref{fig:DeformableModelCollaborationDiagram}.
//
// Software Guide : EndLatex
try
{
// Software Guide : BeginCodeSnippet
gradientMapFilter->Update();
// Software Guide : EndCodeSnippet
}
catch( itk::ExceptionObject & e )
{
std::cerr << "Exception caught when updating gradientMapFilter " <<
std::endl;
std::cerr << e << std::endl;
return -1;
}
std::cout << "The gradient map created!" << std::endl;
// Software Guide : BeginLatex
//
// Now we can construct the mesh source filter that implements the
// isocontouring algorithm.
//
// \index{BinaryMask3DMeshSource!New()}
// \index{BinaryMask3DMeshSource!Pointer}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
MeshSourceType::Pointer meshSource = MeshSourceType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Then we create the filter implementing the deformable model and set
its
// input to the output of the binary mask mesh source. We also set the
// vector field using the \code{SetGradient()} method.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
DeformableFilterType::Pointer deformableModelFilter =
DeformableFilterType::New();
deformableModelFilter->SetGradient( gradientMapFilter->GetOutput() );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Here we connect the output of the binary mask reader to the input of
// the BinaryMask3DMeshSource that will apply the isocontouring algorithm
// and generate the initial mesh to be deformed. We must also select the
// value to be used for representing the binary object in the image. In
// this case we select the value $200$ and pass it to the filter using
its
// method \code{SetObjectValue()}.
//
// \index{itk::BinaryMask3DMeshSource!SetInput()}
// \index{itk::BinaryMask3DMeshSource!SetObjectValue()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
BinaryImageType::Pointer mask = maskReader->GetOutput();
meshSource->SetInput( mask );
//meshSource->SetObjectValue( 200 );
meshSource->SetObjectValue( 1 );
std::cout << "Creating mesh..." << std::endl;
try
{
meshSource->Update();
}
catch( itk::ExceptionObject & excep )
{
std::cerr << "Exception Caught !" << std::endl;
std::cerr << excep << std::endl;
}
deformableModelFilter->SetInput( meshSource->GetOutput() );
// Software Guide : EndCodeSnippet
meshSource->GetOutput()->Print(std::cout);
std::cout << "Deformable mesh created using Marching Cube!" << std::endl;
// Software Guide : BeginLatex
//
// Next, we set the parameters of the deformable model computation.
// \code{Stiffness} defines the model stiffness in the vertical and
// horizontal directions on the deformable surface. \code{Scale} helps to
// accommodate the deformable mesh to gradient maps of different size.
//
// \index{itk::DeformableMesh3DFilter!SetStiffness()}
// \index{itk::DeformableMesh3DFilter!SetScale()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::CovariantVector<double, 2> double2DVector;
typedef itk::CovariantVector<double, 3> double3DVector;
double2DVector stiffness;
stiffness[0] = 0.0001;
stiffness[1] = 0.1;
double3DVector scale;
scale[0] = 1.0;
scale[1] = 1.0;
scale[2] = 1.0;
deformableModelFilter->SetStiffness( stiffness );
deformableModelFilter->SetScale( scale );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Other parameters to be set are the gradient magnitude, the time step
and
// the step threshold. The gradient magnitude controls the magnitude of
the
// external force. The time step controls the length of each step during
// deformation. Step threshold is the number of the steps the model will
// deform.
//
// \index{itk::DeformableMesh3DFilter!SetGradientMagnitude()}
// \index{itk::DeformableMesh3DFilter!SetTimeStep()}
// \index{itk::DeformableMesh3DFilter!SetStepThreshold()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
deformableModelFilter->SetGradientMagnitude( 0.8 );
deformableModelFilter->SetTimeStep( 0.01 );
deformableModelFilter->SetStepThreshold( 60 );
// Software Guide : EndCodeSnippet
std::cout << "Deformable mesh fitting...";
// Software Guide : BeginLatex
//
// Finally, we trigger the execution of the deformable model computation
// using the \code{Update()} method of the DeformableMesh3DFilter. As
// usual, the call to \code{Update()} should be placed in a
// \code{try/catch} block in case any exceptions are thrown.
//
// \index{DeformableMesh3DFilter!Update()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
try
{
deformableModelFilter->Update();
}
catch( itk::ExceptionObject & excep )
{
std::cerr << "Exception Caught !" << std::endl;
std::cerr << excep << std::endl;
}
// Software Guide : EndCodeSnippet
std::cout << "Mesh Source: " << meshSource;
// Software Guide : BeginLatex
//
// The \doxygen{PointSetToImageFilter} takes the deformed
// mesh and produce a binary image corresponding to the node
// of the mesh. Note that only the nodes are producing the image
// and not the cells. See the section on SpatialObjects to produce
// a complete binary image from cells using the
\doxygen{MeshSpatialObject}
// combined with the \doxygen{SpatialObjectToImageFilter}.
// However, using SpatialObjects is computationally more expensive.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::PointSetToImageFilter<MeshType,ImageType> MeshFilterType;
MeshFilterType::Pointer meshFilter = MeshFilterType::New();
meshFilter->SetOrigin(mask->GetOrigin());
meshFilter->SetSize(mask->GetLargestPossibleRegion().GetSize());
meshFilter->SetSpacing(mask->GetSpacing());
//added
//meshFilter->SetDirection(mask->GetDirection());
meshFilter->SetInput(meshSource->GetOutput());
try
{
meshFilter->Update();
}
catch( itk::ExceptionObject & excep )
{
std::cerr << "Exception Caught !" << std::endl;
std::cerr << excep << std::endl;
}
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The resulting deformed binary mask can be written on disk
// using the \doxygen{ImageFileWriter}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::ImageFileWriter<ImageType> WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetInput(meshFilter->GetOutput());
writer->SetFileName(argv[3]);
writer->Update();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Note that in order to successfully segment images, input
// parameters must be adjusted to reflect the characteristics of the
// data. The output of the filter is an Mesh. Users can use
// their own visualization packages to see the segmentation results.
//
// Software Guide : EndLatex
return EXIT_SUCCESS;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100629/572f076e/attachment-0001.htm>
More information about the Insight-users
mailing list