[Insight-users] vessel enhancing diffusion filter (release 200)
Luis Ibanez
luis.ibanez at kitware.com
Wed Jun 18 12:12:27 EDT 2008
Hi Oleksandr,
1) The Hessian matrix is held as an image of Tensors.
It may be possible to rework this code in order to
compute the Hessian values, per pixel as needed.
It will require some coffee and persistence.
You are welcome to give it a try to identifying
the locations of the code that should be modified.
2) It will be interesting to identify the maximum
size of image that currently can be processed in
your machine with this pipeline.
Then you could explore the use of the "ReleaseDataFlag"
as a mechanism for releasing the memory of filter
outputs that are no longer needed.
3) To be pragmatic,
Here is what you may want to do:
a) Build your ITK with ITK_USE_REVIEW enabled.
b) In your test program insert itkMemoryProbes.
Their use is as follows:
#include "itkMemoryProbesCollectorBase.h"
itk::MemoryProbesCollectorBase memorymeter
memorymeter.Start( "Filter1");
filter1->Update();
memorymeter.Stop( "Filter1");
memorymeter.Start( "Filter2");
filter2->Update();
memorymeter.Stop( "Filter2");
memorymeter.Start( "Filter3");
filter3->Update();
memorymeter.Stop( "Filter3");
memorymeter.Report(std::cout);
Note that the order of the filters is important.
the code above assumes that the filters are connected
filter2->SetInput( filter1->GetOutput() );
filter3->SetInput( filter2->GetOutput() );
...
c) run the test again and post the results to the list.
4) The memory leaks reported by Valgrind may not be
real, given that the program crashed and therefore
the normal destructor of ITK classes were not invoked.
Of course, this doesn't mean that there are no memory
leaks at all.... so, let's keep this in mind as a test
to be run again once you get the pipeline to run without
crashing.
Regards,
Luis
------------------------
Oleksandr Dzyubak wrote:
> Hi Luis,
>
> Of course, I meant to thank all the authors/developers. Thanks.
>
> 1) Well, number 6 comes from the Hessian matrix since it is symmetric.
> I meant are there any other arrays should be held/allocated in memory or
> lets put it that way, is it necessary to keep the whole arrays
> permanently?
>
> OK. I agree with you that since I run out of resources, some workarounds
> could be
> a) process ROI
> and/or
> b) sacrifice an accuracy (is that so?) and use float or even short int.
>
> 3) Following your advice a wrote a simple code snippet
> to test memory allocation and the results are below.
>
> ********Begin Test
> dzyubak at debian: /Memory$ ./a.out
> How much memory (in MB) would you like to allocate?
> 2900
> You are able to allocate: 2900 MB memory chunk
> dzyubak at debian: /Memory$ ./a.out 2990
> How much memory (in MB) would you like to allocate?
> 2990
> Error: memory could not be allocated
> dzyubak at debian: /Memory$
> *******End Test*******
>
> Summing up.
> Even though I have 2GB + 2GB(swap), it looks like my limit is 2900MB.
> From my previous calculation (~ 1.2397 GB), I should be good, right?
> Am I missing something?
>
> 5) Memory request.
> As to me, it looks like the memory request was "1,840,443,761 bytes".
> Should be OK, it is not?
> The filter itself does not give the details of the "memory allocation
> crash"
> so I ran it with valgrind (see below). Looks like there is some memory
> leak, right?
>
>
> ***********Begin Valgring********
>
> dzyubak at debian: /Images$ valgrind
> ./itkAnisotropicDiffusionVesselEnhancementImageFilterTest
> H61_8um_100Slices.hdr H61_8um_100Slices_diff_enh_1_1_2.hdr 1 1 2
> ==15766== Memcheck, a memory error detector.
> ==15766== Copyright (C) 2002-2006, and GNU GPL'd, by Julian Seward et al.
> ==15766== Using LibVEX rev 1658, a library for dynamic binary translation.
> ==15766== Copyright (C) 2004-2006, and GNU GPL'd, by OpenWorks LLP.
> ==15766== Using valgrind-3.2.1-Debian, a dynamic binary instrumentation
> framework.
> ==15766== Copyright (C) 2000-2006, and GNU GPL'd, by Julian Seward et al.
> ==15766== For more details, rerun with: -v
> ==15766==
> Reading input image : H61_8um_100Slices.hdr
> ==15766== Warning: set address range perms: large range 198880000
> (undefined)
> Enhancing vessels.........: H61_8um_100Slices.hdr
> ==15766== Warning: set address range perms: large range 198880000
> (undefined)
> ==15766== Warning: set address range perms: large range 198880000
> (undefined)
> ==15766== Warning: set address range perms: large range 1193280004
> (undefined)
> Iteration: 0
> **15766** new/new[] failed and should throw an exception, but Valgrind
> cannot throw exceptions and so is aborting instead. Sorry.
> ==15766== at 0x401D5AD: VALGRIND_PRINTF_BACKTRACE (valgrind.h:2366)
> ==15766== by 0x401D7F5: operator new[](unsigned)
> (vg_replace_malloc.c:195)
> ==15766== by 0x80E2D88: itk::ImportImageContainer<unsigned long,
> itk::SymmetricSecondRankTensor<double, 3> >::AllocateElements(unsigned
> long) const (in
> /mnt/Public/ITK_VTK_Test/Vessel_Enhancement_Diffusion/Images/itkAnisotropicDiffusionVesselEnhancementImageFilterTest)
>
> ==15766== by 0x80C7A7B: itk::ImportImageContainer<unsigned long,
> itk::SymmetricSecondRankTensor<double, 3> >::Reserve(unsigned long) (in
> /mnt/Public/ITK_VTK_Test/Vessel_Enhancement_Diffusion/Images/itkAnisotropicDiffusionVesselEnhancementImageFilterTest)
>
> ==15766== by 0x80C7AF6:
> itk::Image<itk::SymmetricSecondRankTensor<double, 3>, 3>::Allocate() (in
> /mnt/Public/ITK_VTK_Test/Vessel_Enhancement_Diffusion/Images/itkAnisotropicDiffusionVesselEnhancementImageFilterTest)
>
> ==15766== by 0x80C7BF9:
> itk::ImageAdaptor<itk::Image<itk::SymmetricSecondRankTensor<double, 3>,
> 3>, itk::NthElementPixelAccessor<float,
> itk::SymmetricSecondRankTensor<double, 3> > >::Allocate() (in
> /mnt/Public/ITK_VTK_Test/Vessel_Enhancement_Diffusion/Images/itkAnisotropicDiffusionVesselEnhancementImageFilterTest)
>
> ==15766== by 0x80F268C:
> itk::HessianRecursiveGaussianImageFilter<itk::Image<double, 3>,
> itk::Image<itk::SymmetricSecondRankTensor<double, 3>, 3>
> >::GenerateData() (in
> /mnt/Public/ITK_VTK_Test/Vessel_Enhancement_Diffusion/Images/itkAnisotropicDiffusionVesselEnhancementImageFilterTest)
>
> ==15766== by 0x43D5B0E:
> itk::ProcessObject::UpdateOutputData(itk::DataObject*) (in
> /usr/local/lib/InsightToolkit/libITKCommon.so.3.7.0)
> ==15766== by 0x438F750: itk::DataObject::UpdateOutputData() (in
> /usr/local/lib/InsightToolkit/libITKCommon.so.3.7.0)
> ==15766== by 0x438F673: itk::DataObject::Update() (in
> /usr/local/lib/InsightToolkit/libITKCommon.so.3.7.0)
> ==15766== by 0x43D4FDE: itk::ProcessObject::Update() (in
> /usr/local/lib/InsightToolkit/libITKCommon.so.3.7.0)
> ==15766== by 0x80EF733:
> itk::AnisotropicDiffusionVesselEnhancementImageFilter<itk::Image<double,
> 3>, itk::Image<double, 3> >::UpdateDiffusionTensorImage() (in
> /mnt/Public/ITK_VTK_Test/Vessel_Enhancement_Diffusion/Images/itkAnisotropicDiffusionVesselEnhancementImageFilterTest)
>
> ==15766==
> ==15766== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 69 from 1)
> ==15766== malloc/free: in use at exit: 716,721,809 bytes in 14,113 blocks.
> ==15766== malloc/free: 17,083 allocs, 2,969 frees, 1,840,443,761 bytes
> allocated.
> ==15766== For counts of detected errors, rerun with: -v
> ==15766== searching for pointers to 14,113 not-freed blocks.
> ==15766== checked 518,558,692 bytes.
> ==15766==
> ==15766== LEAK SUMMARY:
> ==15766== definitely lost: 0 bytes in 0 blocks.
> ==15766== possibly lost: 119,782,223 bytes in 10,850 blocks.
> ==15766== still reachable: 596,939,586 bytes in 3,263 blocks.
> ==15766== suppressed: 0 bytes in 0 blocks.
> ==15766== Reachable blocks (those to which a pointer was found) are not
> shown.
> ==15766== To see them, rerun with: --show-reachable=yes
> dzyubak at debian: /Images$
>
> *********End Valgrind**********
>
>
> Luis Ibanez wrote:
>
>>
>> Hi Oleksandr,
>>
>>
>> We are glad to hear that you found this Insight Journal paper useful.
>>
>>
>> Just to make the credits straight:
>>
>> The orginal filter was proposed in:
>>
>> "Vessel enhancing diffusion:
>> A scale space representation of vessel structures."
>> by R. Mannieshing, M.A. Viergever, and W. J . Niessen.
>> Medical Image Analysis, 2006.
>>
>> The implementation of this filter in ITK code was done by
>> Andinet Enquobahrie (at Kitware) and supervised by Stephen Aylward
>> (at Kitware).
>>
>>
>> ---------
>>
>> About your question regarding memory requirements:
>>
>> 1) The number "6" in the formula comes from the representation
>> of the Hessian matrix in ITK. Hessian matrices are symmetric
>> and therfore, instead of 9 independent components, they only
>> have 6 independent components. Therefore ITK use only 6 doubles
>> for storing them.
>>
>> 2) Yes, there are more intermediate structures used internally,
>> than the ones used by the Hessian filter.
>>
>>
>> 3) The swap memory should count. That is, you should be able
>> to allocate 4Gb in your machine.
>>
>> That's actually an interesting test to perform before we
>> move further.
>>
>> Can you check in a simple C++ program that you can allocate
>> an array of 4Gb (well... maybe 3.5Gb) of memory ?
>>
>>
>> 4) When the filter was tested, images of 100x100x100 were used.
>> We will have to profile memory use starting from there...
>>
>>
>> 5) It will be useful to know what the actuall memory request was.
>> Did you got that number as part of the error message ?
>>
>> If so, could you please post it to the list ?
>>
>>
>>
>> Thanks
>>
>>
>> Luis
>>
>>
>>
>>
>> --------------------------
>> Oleksandr Dzyubak wrote:
>>
>>> Hi Luis,
>>>
>>> First of all many thanks all of you for such a good filter, "Vessel
>>> enhancing diffusion filter".
>>> With great pleasure I read the article from the the IJ distribution
>>> and the references therein.
>>>
>>> Since I am working on vessels, of course, I was tempted to give that
>>> filter a try. So I did.
>>> As well as Laura, I got the same error message and following your advice
>>> and formula, I calculated the memory request.
>>>
>>> I tested the filter using an image 565x440x100 pixels.
>>> BTW, why 6 in your formula "sizeof(double) x 6 bytes per pixel"?
>>> Do you store some intermediate results all the time?
>>> Lets calculate memory. In my case sizeof(double)=8.
>>>
>>> octave:2> 8*6*(565*440*100)
>>> ans = 1193280000
>>>
>>> OK. Filter + ImageItself = 1.19 + 0.0497 ~ 1.2397 GB
>>>
>>> I have 2 GB + 2 GB (swap). As you say, I almost hit the limit but
>>> some piece is still left.
>>> Does a swap part count? I cropped the original image (which is 100
>>> times larger then the one I used)
>>> down to 47MB just to test the filter and even with such a small image
>>> size the filter fails to allocate memory?
>>>
>>> Does it mean that this filter has no use for boxes with limited
>>> resources?
>>>
>>> BTW, the only
>>> "itkAnisotropicDiffusionVesselEnhancementImageFilterTest" fails.
>>> The other one,
>>> "itkMultiScaleHessianSmoothed3DToVesselnessMeasureImageFilterTest"
>>> works fine.
>>>
>>> Just by luck or you implemented another memory model?
>>>
>>> Thanks,
>>>
>>> Alex
>>>
>>> Luis Ibanez wrote:
>>>
>>>>
>>>> Hi Laura,
>>>>
>>>> This is annoying, but normal.
>>>>
>>>> This code computes Hessians of images, which in 3D requires
>>>> the allocation of sizeof(double) x 6 bytes per pixel.
>>>>
>>>> That is, you will need 48 bytes per pixel of your image
>>>> in order to store the resulting Hessian alone. There will be
>>>> of course additional intermediate allocations.
>>>>
>>>> If you need to process large images you may need a 64bits
>>>> machine with a larger memory...
>>>>
>>>> What is the actual size (in pixels) of the image that
>>>> you are processing ?
>>>>
>>>> Can you process selected regions of the image ?
>>>>
>>>> Usually there is a lot of empty (or at least, non interesting)
>>>> space in medical images. You could use the RegionOfInterest
>>>> filter to reduce the vessel enhancing processing to smaller
>>>> section of the image.
>>>>
>>>>
>>>> Please let us know,
>>>>
>>>>
>>>> Thanks
>>>>
>>>>
>>>> Luis
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> ---------------------------------
>>>> Laura Fernandez de Manuel wrote:
>>>>
>>>>> Hi all,
>>>>>
>>>>> We have been checking the implementation of the "Vessel Enhancing
>>>>> Diffusion Filter" depicted here:
>>>>>
>>>>> http://insight-journal.org/midas/handle.php?handle=1926/558
>>>>>
>>>>> although it works properly with the example 3D images provided
>>>>> (ranging from around 30 to 250 KB) we didn't succeed to make it
>>>>> work in images any larger (5MB images failed already for instance).
>>>>> We work in a system with 4GB RAM so we don't know which can be the
>>>>> source of the "Failed to allocate memory for image" errors that we
>>>>> get. Here, I attach the error message we get:
>>>>>
>>>>> ------------------------------------------------------------------------------------------------------
>>>>>
>>>>> ./itkAnisotropicDiffusionVesselEnhancementImageFilterTest.exe
>>>>> image00_2.mhd image_2Enhanced.mhd
>>>>> Reading input image : image00_2.mhd
>>>>> Enhancing vessels.........: image00_2.mhd
>>>>> Iteration: 0
>>>>> Computing vesselness for scale with sigma= 0.2
>>>>> Exception caught:
>>>>> itk::ExceptionObject (0138FB20)
>>>>> Location: "class itk::SymmetricSecondRankTensor<double,3>
>>>>> *__thiscall itk::ImportImageContainer<unsigned long,class
>>>>> itk::SymmetricSecondRankTensor<double,3> >:: AllocateElements
>>>>> (unsigned long) const"
>>>>> File: itk3.6.0\code\common\itkImportImageContainer.txx
>>>>> Line: 193
>>>>> Description: Failed to allocate memory for image.
>>>>> ------------------------------------------------------------------------------------------------------
>>>>>
>>>>>
>>>>> Thanks a lot!
>>>>>
>>>>> carlos & laura
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> ------------------------------------------------------------------------
>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> Insight-users mailing list
>>>>> Insight-users at itk.org
>>>>> http://www.itk.org/mailman/listinfo/insight-users
>>>>
>>>>
>>>> _______________________________________________
>>>> Insight-users mailing list
>>>> Insight-users at itk.org
>>>> http://www.itk.org/mailman/listinfo/insight-users
>>>
>>>
>>>
>>>
>
>
More information about the Insight-users
mailing list