[Insight-users] Help with extracting curves from an angiogram

[Luis Ibanez]

Hi Ifeoma,

Thanks for posting your code.
It makes a lot easier to look at your problem.

Here are some comments:

1) Since you are creating the image, you can initialize
    the iterator by doing:

    IteratorType iterator( m_InImage, region );

    instead of

    IteratorType iterator(m_InImage, m_InImage->GetRequestedRegion());


2) You don't seem to be setting the Origin or the Spacing
    of the image. This values are very important for curve
    extraction.


3) Please double check that in the for loop where you
    copy the pixel data, X corresponds to width, and Y
    corresponds to height, as well as, verify that the
    organization of the data is such that the Y axis
    is the one running faster:

    You have now:

 >    for(x=0; x<width; x++)
 >    {
 >       for(y=0; y<height; y++)
 >       {
 >          iterator.Set(pixel_data[x][y]);
 >          ++iterator;
 >       }
 >    }


     It is more common to have something like:


       for(y=0; y<height; y++)
         {
         for(x=0; x<width; x++)
           {
           iterator.Set( pixel_data[y][x] );
           ++iterator;
           }
         }


     Your code may be right, but is still worth to
     check because if you are swapping the dimensions,
     you may be going out of the bound in your pixel
     data array, and the consequences may show up in
     any random place of the execution later.



4)  Please try replacing

      m_ScalarProduct->SetInput2( m_Eigen->GetMaxEigenVector() );

     (the line that crashes)

      with

      typedef EigenVectorFilterType::EigenVectorImageType
                                             EigenVectorImageType;

      EigenVectorImageType * eigenV1 = m_Eigen->GetMaxEigenVector();

      m_ScalarProduct->SetInput2( eigenV1 );


      Run it and check in what line it crashes now.
      (and let us know what you find).

      The only reason for a crash in GetMaxEigenVector() is if
      the internal dynamic_cast<> fails. Which is unlikely to
      occurr.



Still, the main suspect is your process for copying the data into
        m_InImage.

A simple way to verify this is to *comment out* the copying of
pixel data, and simply call


       m_InImage->FillBuffer( 0 );

Then run the program agin. If the program doesn't crash, then you
will know that the pixel data copying is the culprit. It the program
still crashes... then we have to look for other suspects.






     Regards,


        Luis




-------------------
Ifeoma Nwogu wrote:
> Hi Luis,
> 
> Sorry to burden you - I already tried what you suggested although not with the 
> ImageFileWriter. After loading the image, I simply passed through its contents 
> with an iterator and checked iterator.Value() and it was correct all the way.
> 
> But I'm attaching my file, hopefully it should be easy to read through since 
> its mostly the 2D Curve Extraction code. Thanks so much for your help, I was 
> at quite a loss.
> 
> I will also try the Try/Catch block but I'm not sure if that will tell me 
> exactly what the problem is. Thanks.
> 
> Ifeoma
> 
> 
> Quoting Luis Ibanez <luis.ibanez at kitware.com>:
> 
> 
>>Hi Ifeoma,
>>
>>Since the only thing that you are doing differently is to write
>>your image data into an image, that's probably what you are
>>doing wrong.
>>
>>
>>
>>Could you please post the code that you are using for loading
>>your image into an ITK image ?
>>
>>
>>
>>Also, Please add a ImageFileWriter and try to save your ITK image
>>as soon as you are done with the process of importing it from
>>memory. In this way you will be able to double check if the
>>image was imported correctly.
>>
>>
>>You should also add a Try/Catch block around the line where you
>>are detecting the crash. It is quite likely that an exception is
>>being thrown.  You will find many examples on how to write try/catch
>>block in the ITK Software Guide
>>
>>         http://www.itk.org/ItkSoftwareGuide.pdf
>>
>>
>>
>>
>>BTW:
>>
>>
>>           WHAT VERSION OF ITK ARE YOU USING ?
>>
>>
>>
>>
>>Please let us know,
>>
>>
>>
>>     Thanks,
>>
>>
>>        Luis
>>
>>
>>--------------------
>>Ifeoma Nwogu wrote:
>>
>>>Hello,
>>>
>>>I'm trying to use the code that was written code for 2D Curve Extraction 
>>>application (ceExtractorConsoleBase and ceExtractorConsoleBase classes) but
>>
>>>once I try to call the line 
>>>
>>>m_ScalarProduct->SetInput2(m_Eigen->GetMaxEigenVector()); 
>>>and m_Eigen is of type 
>>>itk::EigenAnalysis2DImageFilter< ImageType, ImageType, VectorImageType > 
>>>  
>>>my application crashes because of the m_Eigen filter! I am using VC++ 7.1
>>
>>>compiler, I went through the tutorials online and successfully implemented
>>
>>a 
>>
>>>simpe ITK program.
>>>
>>>The only thing I think I'm doing differently from the code online is
>>
>>writing 
>>
>>>my pixel values into an image pointer instead of using a reader (my images
>>
>>are 
>>
>>>already loaded in memory). 
>>>
>>>Please let me know if you have any ideas on what could be wrong here. Thank
>>
>>>you!
>>>
>>>Ifeoma
>>>
>>>
>>>
>>>
>>>
>>>_______________________________________________
>>>Insight-users mailing list
>>>Insight-users at itk.org
>>>http://www.itk.org/mailman/listinfo/insight-users
>>>
>>>
>>
> 
> 
> 
> 
> ------------------------------------------------------------------------
> 
> #if defined(_MSC_VER)
> #pragma warning ( disable : 4786 )
> #endif
> 
> #include "ITKImageResponseMap.h"
> 
> /************************************
>  *  Constructor
>  ***********************************/
> ITKImageResponseMap::ITKImageResponseMap(ImageData *img)
> {
>    int x, y;
>    float** pixel_data = img->getPixelData();
>    int width = img->getWidth();
>    int height = img->getHeight();
> 
>    m_RescaleIntensitySmoothed   = RescaleIntensityFilterType::New();
>    m_RescaleIntensityMedialness = RescaleIntensityFilterType::New();
>    m_RescaleIntensityMaxEigen   = RescaleIntensityFilterType::New();*/
> 
>    // Initialize the zeroth order directional Gaussian filters
>    m_Hx = InputGaussianFilterType::New();
>    m_Hy = InputGaussianFilterType::New();
>    m_Hx->SetDirection(0);
>    m_Hy->SetDirection(1);
>    
>    m_Hxy = GaussianFilterType::New();
>    m_Hxy->SetDirection(1);
> 
>    // Initialize the first order directional Gaussian filters
>    m_H1x = GaussianFilterType::New();
>    m_H1y = GaussianFilterType::New();
>    m_H1x->SetDirection(0);
>    m_H1y->SetDirection(1);
>    m_H1x->SetOrder(GaussianFilterType::FirstOrder);
>    m_H1y->SetOrder(GaussianFilterType::FirstOrder);
>    m_H1xy = GaussianFilterType::New();
>    m_H1xy->SetDirection(1);
>    m_H1xy->SetOrder(GaussianFilterType::FirstOrder);
>       
>    // Initialize the second order directional Gaussian filters
>    m_H2x = GaussianFilterType::New();
>    m_H2y = GaussianFilterType::New();
>    m_H2x->SetDirection( 0 );
>    m_H2y->SetDirection( 1 );
>    m_H2x->SetOrder( GaussianFilterType::SecondOrder );
>    m_H2y->SetOrder( GaussianFilterType::SecondOrder );
> 
>    //Loading image data
>    // Create an ITK input image and write contents from ImageData*
>    m_InImage = ImageType::New();
>    ImageType::IndexType start;
>    ImageType::SizeType size;
>    start[0] = 0;
>    start[1] = 0;
>    size[0] = width;
>    size[1] = height;
> 
>    ImageType::RegionType region;
>    region.SetSize(size);
>    region.SetIndex(start);
>    m_InImage->SetRegions(region);
>    m_InImage->Allocate();
> 
>    IteratorType iterator(m_InImage, m_InImage->GetRequestedRegion());
> 
>    //populate itk image with loaded image data
>    for(x=0; x<width; x++)
>    {
>       for(y=0; y<height; y++)
>       {
>          iterator.Set(pixel_data[x][y]);
>          ++iterator;
>       }
>    }
> 
>    // Set filter sigma value
>    this->SetSigma(6);
> 
>    // Applying the filters on the input image and its derivatives
>    m_Hx->SetInputImage(m_InImage);
>    m_Hy->SetInputImage(m_InImage);
>    m_Hxy->SetInputImage(m_Hx->GetOutput());
> 
>    m_H1x->SetInputImage(m_Hy->GetOutput());
>    m_H1y->SetInputImage(m_Hx->GetOutput());
> 
>    m_H2x->SetInputImage(m_Hy->GetOutput());
>    m_H2y->SetInputImage(m_Hx->GetOutput());
>    
>    m_H1xy->SetInputImage(m_H1x->GetOutput());
> 
>    //Combine the xx and yy filters to obtain part of the Hessian matrix
>    m_Add = AddFilterType::New();
>    m_Add->SetInput1( m_H2x->GetOutput() );
>    m_Add->SetInput2( m_H2y->GetOutput() );
> 
>    //Define the mixed derivatives 
>    m_Modulus = ModulusFilterType::New();
>    m_Modulus->SetInput1(m_H1x->GetOutput());
>    m_Modulus->SetInput2(m_H1y->GetOutput());
>   
>    //The Gradient filters are defined, to get the gradient projection on the eigen vector
>    m_Gradient = GradientFilterType::New();  
>    m_Gradient->SetInput(m_InImage);
> 
>    //The EigenAnalysis2D filters are used to compute pixel-wise eigen values and vectors
>    m_Eigen = EigenFilterType::New();
>    m_Eigen->SetInput1(m_H2x->GetOutput());
>    m_Eigen->SetInput2(m_H1xy->GetOutput());
>    m_Eigen->SetInput3(m_H2y->GetOutput());
>       
>    //The join filter combines two images, to get an image in which each pixel contains components 
>    // of the first image followed by the components of the second image
>    m_Join = JoinFilterType::New();
>    m_Join->SetInput1( m_H1x->GetOutput() );
>    m_Join->SetInput2( m_H1y->GetOutput() );
>   
>    // Scalar product filters show the gradient projected on the eigen vectors
>    m_ScalarProduct = ScalarProductFilterType::New();
>    m_ScalarProduct->SetInput1(m_Join->GetOutput());
>    /m_ScalarProduct->SetInput2(m_Eigen->GetMaxEigenVector()); //CRASHES HERE!
>  
>    // Normalize the parametric space
>    m_RescaleIntensitySmoothed->SetInput(m_Hxy->GetOutput());
>    m_RescaleIntensitySmoothed->SetOutputMinimum( -1.0 );
>    m_RescaleIntensitySmoothed->SetOutputMaximum(  1.0 );
> 
>    m_RescaleIntensityMedialness->SetInput( m_ScalarProduct->GetOutput() );
>    m_RescaleIntensityMedialness->SetOutputMinimum( -1.0 );
>    m_RescaleIntensityMedialness->SetOutputMaximum(  1.0 );
> 
>    m_RescaleIntensityMaxEigen->SetInput(m_Eigen->GetMaxEigenValue());
>    m_RescaleIntensityMaxEigen->SetOutputMinimum(0.0);
>    m_RescaleIntensityMaxEigen->SetOutputMaximum(1.0);
> 
>    //Create a structure to store the max eigen values from the eigen filter
>    ImageType::Pointer eigenImage = ImageType::New();
>    ImageType::IndexType start;
>    ImageType::SizeType size;
>    start[0] = 0;
>    start[1] = 0;
>    size[0] = width;
>    size[1] = height;
>    ImageType::RegionType region;
>    region.SetSize(size);
>    region.SetIndex(start);
>    eigenImage->SetRegions(region);
>    eigenImage->Allocate();
>    eigenImage = m_Eigen->GetMaxEigenValue();*/
>    ImageType::Pointer eigenImage = m_RescaleIntensityMaxEigen->GetOutput();
> 
>    IteratorType eigenIt(eigenImage, eigenImage->GetRequestedRegion());
>    eigenIt.GoToBegin();
> 
>    //populate image data pixels with max eigen values
>    for(x=0; x<width; x++)
>    {
>       for(y=0; y<height; y++)
>       {
>          ImageType::PixelType pixel2D = eigenIt.Value(); //CRASH OCCURS HERE!
>          printf("Eigen data for x=%d, y=%d is %f\n",x,y, (double)pixel2D);
>          ++eigenIt;
>       }
>    }
> }
> 
> void ITKImageResponseMap::SetSigma(RealType value)
> {  
>    m_Hx->SetSigma(value);
>    m_Hy->SetSigma(value);
> 
>    m_Hxy->SetSigma(value);
> 
>    m_H1x->SetSigma(value);
>    m_H1y->SetSigma(value);
> 
>    m_H1xy->SetSigma(value);
> 
>    m_H2x->SetSigma(value);
>    m_H2y->SetSigma(value);
> 
>    //m_Gradient->SetSigma(value);
> }
> 
> void ITKImageResponseMap::Execute()
> {  
>    m_Hxy->Update();
>    m_H1xy->Update();
>    m_Add->Update();
>    m_Modulus->Update(); 
> }