[Insight-users] Re: Help with extracting curves from an angiogram
Luis Ibanez
luis.ibanez at kitware.com
Thu Dec 8 21:44:30 EST 2005
Hi Ifeoma,
Thanks for doing the suggested tests.
This looks quite strange.
Does it happens with any image ?
Please do the following:
Replace the Image that you are creating with an ImageFileReader
and feed your program with an image from the Examples directory
in ITK.
Note that you will have to comment out any references to width
and height.
Please let us know what happens when you use a reader
and a standard image as the input for your code.
Thanks
Luis
---------------------
Ifeoma Nwogu wrote:
> Hello,
>
> I had sent a message out earlier in the week, still no luck - my windows
> application (VC++ 7.1 compiler) constantly crashes when I call a filter->Update
> (). Is this a problem you have encoutered with folks in the past?
>
> I tried reinstalling an older version of ITK - from 2.2.0, I downgraded to
> 1.8.1 but no difference. The application still crashes at the same two points
> filter->Update() or filter->GetMaxEigenVector();
>
> Any more suggestions or ideas are very welcome. Thanks much.
>
> Ifeoma
>
> P.S. I did the steps suggested below successfully.
>
>
>>
>>Quoting Luis Ibanez <luis.ibanez at kitware.com>:
>>
>>
>>>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();
>>>>}
>>>
>>
>>
>>
>
>
>
>
>
More information about the Insight-users
mailing list