[Insight-users] Fwd: DeformableMesh3DFilter Problems

[Federico J. Lopez Bertoni]

Hi,

I've been trying to implement the DeformableMesh3DFilter and run into some
problems. I'm basically using a variation of the
DeformableModel1.cxxalgorithm in which I read in a vesselness image
(which I get using itk's
vesselness filter) and input his image to the gradient filter (I bypass the
gradient magnitude calculation since this image is giving a maximum value at
maximum vesselness measure). Instead of creating a mesh from a binary image,
I read in points in 3D and create a polyline mesh from them (this is my
initial mesh). I am  using an open mesh as an initial guess that should be
deformed to the vesselness image (which in my application corresponds to the
coronary arteries). The program compiles fine but when I run it, it crashes
when it attempts to do the deformable mesh filter.  Is there another way to
tackle this problem or solve the issue in this application? Anyway, below is
a copy of my code.

I really appreciate your help and time in advance.


Best Regards,

-- 
Federico J. Lopez Bertoni
Graduate Research Assistant
Duke University
Department of Biomedical Engineering
Office:  1373 CIEMAS
Phone:  919-660-5125
fjlb at duke.edu

#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include "itkMesh.h "
#include " itkVertexCell.h"
#include "itkLineCell.h"
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <stdio.h>
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkGradientRecursiveGaussianImageFilter.h"
//#include " itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkDeformableMesh3DFilter.h"
#include "itkImage.h"
#include "itkCovariantVector.h"
#include "itkPointSetToImageFilter.h "
using namespace std;



int main( int argc, char *argv[] )
{
    if( argc < 4 )
    {
        std::cerr << "Missing Parameters " << std::endl;
        std::cerr << "Usage: VesselnessImage RCApoints ImageOut" <<
std::endl;
        return EXIT_FAILURE;
    }

    //Read RCA, LCA points in.
    ifstream infile;
    infile.open( argv[2] , ios::in );
    if (!infile.good())
    {
        cout << "Error Opening file " << endl;
        exit (1);
    }
    int X[200], Y[200], Z[200], rad[200], i=0;
    while(!infile.eof())
    {
        infile >> X[i] >> Y[i] >> Z[i] >> rad[i];
        if (X[i] == -1 )
            break;
        else
        //cout << i << "\t" << X[i] << "\t" << Y[i] << "\t" << Z[i] << "\t"
<< rad[i] << endl;
        i++;
    }

    /*ifstream infile1;
    infile1.open( argv[3] , ios::in );
    if (!infile1.good())
    {
        cout << "Error Opening file "<< endl;
        exit (1);
    }
    int X1[200], Y1[200], Z1[200], rad1[200], i1=0;
    while(!infile1.eof())
    {
        infile1 >> X1[i1] >> Y1[i1] >> Z1[i1] >> rad1[i1];
        if (X1[i1] == -1 )
            break;
        else
        cout << i1 << "\t" << X1[i1] << "\t" << Y1[i1] << "\t" << Z1[i1] <<
"\t" << rad1[i1] << endl;
        i1++;
    }*/

    //typedef unsigned short PixelType;
    typedef double PixelType;
    const unsigned int Dimension = 3;
    typedef itk::Mesh< PixelType, Dimension >         MeshType;
    typedef MeshType::CellType                CellType;
    typedef itk::VertexCell< CellType >        VertexType;
    typedef itk::LineCell< CellType >       LineType;
    typedef itk::Image< PixelType, Dimension >      ImageType;

    MeshType::Pointer mesh = MeshType::New();
    MeshType::PointType    points;
    //Create the mesh
    for(unsigned int j=0; j<92; j++)
    {
        points[0] = X[j];
        points[1] = Y[j];
        points[2] = Z[j];
        mesh->SetPoint( j, points );
    }
    /*for(unsigned int j=0; j<92; j++)
    {
        points[0][0] = X[j];
        points[0][1] = Y[j];
        points[0][2] = Z[j];
        mesh->SetPoint( j, points );
    }*/

    CellType::CellAutoPointer cellpointer;
    for(unsigned int j=0; j<91; j++)
    {
        cellpointer.TakeOwnership( new LineType );
        cellpointer->SetPointId( 0, j );
        cellpointer->SetPointId( 0, j+1 );
        mesh->SetCell( j, cellpointer );

        cellpointer.TakeOwnership( new VertexType );
        cellpointer->SetPointId( 0, j );
        mesh->SetCell( j + 90, cellpointer );
    }

    // Print out the number of points and the number of cells.
    std::cout << "# Points= " << mesh->GetNumberOfPoints() << std::endl;
    std::cout << "# Cell  = " << mesh->GetNumberOfCells() << std::endl;

    typedef MeshType::PointsContainer::ConstIterator  PointIterator;
    PointIterator pointIterator = mesh->GetPoints()->Begin();
    PointIterator pointEnd      = mesh->GetPoints()->End();

    while( pointIterator != pointEnd )
    {
        //std::cout << pointIterator.Value() << std::endl;
        ++pointIterator;
    }

    typedef MeshType::CellsContainer::ConstIterator  CellIterator;

    CellIterator cellIterator = mesh->GetCells()->Begin();
    CellIterator cellEnd      = mesh->GetCells()->End();

    while( cellIterator != cellEnd )
    {
        CellType * cell = cellIterator.Value();
        //std::cout << cell->GetNumberOfPoints() << std::endl;
        ++cellIterator;
    }

    cellIterator = mesh->GetCells()->Begin();
    cellEnd      = mesh->GetCells()->End();

    while( cellIterator != cellEnd )
    {
        CellType * cell = cellIterator.Value();
        //std::cout << "cell with " << cell->GetNumberOfPoints();
        //std::cout << " points   " << std::endl;

    // Software Guide : BeginCodeSnippet
        typedef CellType::PointIdIterator     PointIdIterator;

        PointIdIterator pointIditer = cell->PointIdsBegin();
        PointIdIterator pointIdend  = cell->PointIdsEnd();

    while( pointIditer != pointIdend )
    {
        //std::cout << *pointIditer << std::endl;
        ++pointIditer;
    }
    // Software Guide : EndCodeSnippet

        ++cellIterator;
    }

    //Mesh creation done at this point.
    std::cout << "PolyLine mesh created..." << std::endl;
    std::cout << "Read in Input Image..." << std::endl;
    //Now read in Vesselness image.
    //Read in Vesselness (Cost) Image
    typedef itk::ImageFileReader< ImageType > ReaderType;
    ReaderType::Pointer reader = ReaderType::New();
    reader->SetFileName( argv[1] );
    //Set image to variable
    ImageType::Pointer image = ImageType::New();
    image = reader->GetOutput();
    //Define Gradient Filter Type
    typedef itk::CovariantVector< double, Dimension >    GradientPixelType;
    typedef itk::Image< GradientPixelType, Dimension >    GradientImageType;
    typedef itk::GradientRecursiveGaussianImageFilter< ImageType,
GradientImageType >        GradientFilterType;
    //typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<
ImageType,ImageType >    GradientMagnitudeFilterType;
    GradientFilterType::Pointer gradientMapFilter =
GradientFilterType::New();
    //GradientMagnitudeFilterType::Pointer gradientMapFilter =
GradientMagnitudeFilterType::New();
    //Set input of gradient filter to be vesselnes image
    gradientMapFilter->SetInput( reader->GetOutput() );
    gradientMapFilter->SetSigma( 1.0 );
    //trigger execution
    try
    {
        gradientMapFilter->Update();
    }
    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 was created successfully" << std::endl;

    // Mesh Source filter
    //MeshSourceType::Pointer meshSource = MeshSourceType::New();

    //Deformable model filter
    typedef itk::DeformableMesh3DFilter< MeshType, MeshType >
DeformableFilterType;
    DeformableFilterType::Pointer deformableModelFilter =
DeformableFilterType::New();
    //Set input to deformable model filter
    deformableModelFilter->SetGradient( gradientMapFilter->GetOutput() );
    deformableModelFilter->SetInput( mesh );
    // Set parameters for deformable mesh
    typedef itk::CovariantVector< double, 2 >    double2DVector;
    typedef itk::CovariantVector< double, 3 >    double3DVector;

    double2DVector stiffness;
    stiffness[0] = 0.01;
    stiffness[1] = 0.1;

    double3DVector scale;
    scale[0] = 1.0;
    scale[1] = 1.0;
    scale[2] = 1.0;

    deformableModelFilter->SetStiffness( stiffness );
    deformableModelFilter->SetScale( scale );
    deformableModelFilter->SetGradientMagnitude( 0.8 );
    deformableModelFilter->SetTimeStep( 0.01 );
    deformableModelFilter->SetStepThreshold( 60 );

    std::cout << "Deformable mesh fitting..." << std::endl;

    //Trigger execution of deformable mesh filter
    try
    {
        deformableModelFilter->Update();
    }
    catch( itk::ExceptionObject & excep )
    {
        std::cerr << "Exception caught while using the deformable mesh
filter" << std::endl;
        std::cerr << excep << std::endl;
    }

    std::cout << "Mesh Source: " << std::endl;

    typedef itk::PointSetToImageFilter< MeshType, ImageType >
MeshFilterType;
    MeshFilterType::Pointer meshFilter = MeshFilterType::New();
    meshFilter->SetOrigin( image->GetOrigin() );
    meshFilter->SetSize( image->GetLargestPossibleRegion().GetSize() );
    meshFilter->SetSpacing( image->GetSpacing() );
    meshFilter->SetInput( mesh );

    try
    {
        meshFilter->Update();
    }
    catch( itk::ExceptionObject & excep )
    {
        std::cerr << "Exception Caught!" << std::endl;
        std::cerr << excep << std::endl;
    }

    typedef itk::ImageFileWriter< ImageType > WriterType;
    WriterType::Pointer writer = WriterType::New();
    writer->SetInput( meshFilter->GetOutput() );
    writer->SetFileName( argv[3] );
    writer->Update();

    return EXIT_SUCCESS;
}




-- 
Federico J. Lopez Bertoni
Graduate Research Assistant
Duke University
Department of Biomedical Engineering
Office:  1373 CIEMAS
Phone:  919-660-5125
fjlb at duke.edu



-- 
Federico J. Lopez Bertoni
Graduate Research Assistant
Duke University
Department of Biomedical Engineering
Office:  1373 CIEMAS
Phone:  919-660-5125
fjlb at duke.edu
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20071204/954b8574/attachment.html