[Insight-users] itkStreamingImageFilter with itkHessianSmoothed3DToVesselnessMeasureImageFilter
Iván Macía
imacia at vicomtech.org
Thu Nov 6 04:41:18 EST 2008
Hi again Dan,
I was thinking that maybe calculations are correct but we have a problem
with the sampling (Nyquist theorem).
If the variance of a Gaussian tends to zero, the Gaussian tends to a pulse
function at the origin, and you are calculating a second order discrete
derivative of this.
Maybe the result is correct but useless. Just wondering, I could be wrong
and there is really an error in the calculations (or both things!).
Regards
Ivan
-----Mensaje original-----
De: Dan Mueller [mailto:dan.muel at gmail.com]
Enviado el: miércoles, 05 de noviembre de 2008 10:14
Para: Iván Macía
CC: insight-users at itk.org
Asunto: Re: [Insight-users] itkStreamingImageFilter with
itkHessianSmoothed3DToVesselnessMeasureImageFilter
Hi Iván,
I have been doing some tests with the discrete Gaussian operator and
image functions you contributed.
> I just wanted to know if someone has made any test to see if the problems
> with the derivatives and Hessian calculation have been fixed after latest
> changes and results are coherent.
I am encountering some problems with GaussianDerivativeOperator
related to image spacing/scale. I have put together a simple test
program (see below) to demonstrate the issue. The issue can be
observed for images of small spacing. Notice in the last line of the
output below all the operator coefficients are negative, which in my
limited understanding is a problem. If you -- or anyone else -- have
any ideas I'd appreciate to hear from you.
D:\Temp\HessianTest\Build\release>OperatorTest
order=2 spacing=5 sigma=5 variance=25 variance(adjusted)=1
0.402 1.83 4.95 4.99 -0.466 4.99 4.95 1.83 0.402
order=2 spacing=4 sigma=4 variance=16 variance(adjusted)=1
0.257 1.17 3.15 3.12 -0.466 3.12 3.15 1.17 0.257
order=2 spacing=3 sigma=3 variance=9 variance(adjusted)=1
0.144 0.653 1.75 1.66 -0.466 1.66 1.75 0.653 0.144
order=2 spacing=2 sigma=2 variance=4 variance(adjusted)=1
0.0635 0.285 0.749 0.624 -0.466 0.624 0.749 0.285 0.0635
order=2 spacing=1 sigma=1 variance=1 variance(adjusted)=1
0.0151 0.0653 0.15 0 -0.466 0 0.15 0.0653 0.0151
order=2 spacing=0.5 sigma=0.5 variance=0.25 variance(adjusted)=1
0.00302 0.0102 0 -0.156 -0.466 -0.156 0 0.0102 0.00302
order=2 spacing=0.1 sigma=0.1 variance=0.01 variance(adjusted)=1
-0.000846 -0.00742 -0.048 -0.206 -0.466 -0.206 -0.048 -0.00742
-0.000846 <<< all negative!
> Could someone also test the performance?
I put together a simple program (see below) to test the performance of
DiscreteHessianGaussianImageFunction (the test image is in the ITK
Testing\Data\Input folder). As you assumed the recursive approach is
much faster:
HessianTest HeadMRVolumeCompressed.mha
Filter took 0.0792656 seconds
Function took 17.0815 seconds
Image Size = [48, 62, 42]
Image Spacing = [4.0, 4.0, 4.0]
Sigma = 1.0
MaxKernelWidth = 32
MaximumError = 0.01
NormalizeAcrossScale = false
UseImageSpacing = false
Speed Factor = 17.08 / 0.079 = 216
(Windows Vista SP1, Visual Studio 2005, Release build, Intel Core 2
Duo 2GHz, 3 GB RAM)
So to make the new function pay off (in terms of speed) it seems a
user will have to restrict the computation to an image region ~1/250
of the whole input image. For my use cases this is a certainty, but I
can image for some (most?) this will not be the case.
Any feedback would be appreciated.
Regards, Dan
==== CMakeLists.txt ====
PROJECT( DiscreteGaussianTest )
CMAKE_MINIMUM_REQUIRED( VERSION 2.6 )
# Find required packages
FIND_PACKAGE(ITK REQUIRED)
# Use required packages
INCLUDE(${ITK_USE_FILE})
# Add this as include directory
INCLUDE_DIRECTORIES(
${CMAKE_SOURCE_DIR}
${SOURCE_PATH}
${VXL_INCLUDE_DIRS}
)
# Executables
ADD_EXECUTABLE(OperatorTest test_operator.cxx)
TARGET_LINK_LIBRARIES(OperatorTest ${ITK_LIBRARIES})
ADD_EXECUTABLE(HessianTest test_hessian.cxx)
TARGET_LINK_LIBRARIES(HessianTest ${ITK_LIBRARIES})
==== test_operator.cxx ====
#include <iostream>
#include <iomanip>
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkGaussianDerivativeOperator.h"
int main(int argc, char* argv[])
{
// Typedefs
const unsigned int Dimension = 2;
typedef float PixelType;
typedef itk::Image< PixelType, Dimension > ImageType;
typedef itk::GaussianDerivativeOperator< PixelType, Dimension >
OperatorType;
// Create operator
OperatorType op;
op.SetDirection( 0 );
op.SetMaximumKernelWidth( 64 );
op.SetMaximumError( 0.001 );
op.SetNormalizeAcrossScale( true );
op.SetUseDerivativeOperator( false );
op.SetOrder( 2 );
double spacings[] = { 5.0, 4.0, 3.0, 2.0, 1.0, 0.5, 0.1 };
for ( unsigned int j=0; j<7; j++)
{
op.SetVariance( spacings[j]*spacings[j] );
op.SetSpacing( spacings[j] );
op.CreateDirectional( );
// Output operator info
std::cout << std::setprecision( 3 );
std::cout << "order=" << op.GetOrder();
std::cout << " spacing=" << spacings[j];
std::cout << " sigma=" << spacings[j];
std::cout << " variance=" << spacings[j]*spacings[j];
std::cout << " variance(adjusted)=" << op.GetVariance();
std::cout << std::endl;
for (unsigned int i=0; i<op.Size(); i++)
std::cout << op[i] << " ";
std::cout << std::endl;
}
return EXIT_SUCCESS;
}
==== test_hessian.cxx ====
#include <iostream>
#include <iomanip>
#include "itkTimeProbe.h"
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkDiscreteHessianGaussianImageFunction.h"
#include "itkHessianRecursiveGaussianImageFilter.h"
int main(int argc, char* argv[])
{
// Typedefs
const unsigned int Dimension = 3;
typedef float PixelType;
typedef itk::Image< PixelType, Dimension > ImageType;
typedef itk::ImageFileReader< ImageType > ReaderType;
typedef itk::ImageRegionIteratorWithIndex< ImageType > IteratorType;
typedef itk::DiscreteHessianGaussianImageFunction< ImageType >
HessianFunctionType;
typedef itk::HessianRecursiveGaussianImageFilter< ImageType >
HessianFilterType;
// Parameters
if ( argc != 2 )
std::cout << "USAGE: " << argv[0] << " InputFilename" << std::endl;
char * InputFilename = argv[1];
double sigma = 1.0;
// Read input image
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName( InputFilename );
reader->Update();
ImageType::Pointer input = reader->GetOutput();
input->DisconnectPipeline();
// Compute hessian using recursive filter
itk::TimeProbe probe1;
HessianFilterType::Pointer hessianFilter = HessianFilterType::New();
hessianFilter->SetNormalizeAcrossScale( false );
hessianFilter->SetSigma( sigma );
hessianFilter->SetInput( input );
probe1.Start( );
hessianFilter->Update();
probe1.Stop( );
input->DisconnectPipeline( );
std::cout << "Filter took " << probe1.GetMeanTime() << " seconds"
<< std::endl;
// Compute hessian using function
itk::TimeProbe probe2;
HessianFunctionType::Pointer hessianFunction =
HessianFunctionType::New();
hessianFunction->SetInputImage( input );
hessianFunction->SetSigma( sigma );
hessianFunction->SetMaximumError( 0.01 );
hessianFunction->SetMaximumKernelWidth( 32 );
hessianFunction->SetNormalizeAcrossScale( false );
hessianFunction->SetUseImageSpacing( false );
hessianFunction->Initialize( );
HessianFunctionType::OutputType out;
IteratorType it( input, input->GetLargestPossibleRegion() );
probe2.Start( );
for (it.Begin(); !it.IsAtEnd(); ++it)
out = hessianFunction->EvaluateAtIndex( it.GetIndex() );
probe2.Stop( );
std::cout << "Function took " << probe2.GetMeanTime() << "
seconds" << std::endl;
return EXIT_SUCCESS;
}
2008/11/4 Iván Macía <imacia at vicomtech.org>:
> Dear all,
>
> I have been out of the office the whole month (honeymoon ;) ) so i could
not
> follow the conversation. This last week I have been in touch with Dan so
he
> could move the contribution to the Review directory.
>
> I just wanted to know if someone has made any test to see if the problems
> with the derivatives and Hessian calculation have been fixed after latest
> changes and results are coherent. As I commented to Moti thanks to his
> comments I found a bug in the Hessian calculation.
>
> Could someone also test the performance? It is obvious that the recursive
> Gaussian should be faster but takes a lot of memory and when calculating
> only a few points the discrete approach should be better. Calculation
times
> depend also considerably on the size of the kernel but this affects also
the
> precision. A trade-off would be the right choice.
>
> Thanks for your help
>
> Ivan
>
>
> -----Mensaje original-----
> De: Karthik Krishnan [mailto:karthik.krishnan at kitware.com]
> Enviado el: lunes, 13 de octubre de 2008 15:07
> Para: Dan Mueller
> CC: Moti Freiman; insight-users at itk.org; Iván Macía
> Asunto: Re: [Insight-users] itkStreamingImageFilter with
> itkHessianSmoothed3DToVesselnessMeasureImageFilter
>
> Dan Mueller wrote:
>> Hi Karthik, Moti, others,
>>
>> Out of interest, I am in the process of transferring Ivan's
>> contribution to the ITK Review folder: I already have the filters and
>> some tests in my local copy, I just have to add a few more tests. I
>> expect to commit the files by week's end.
>>
> Perfect. Thanks a lot Dan. This is very helpful.
>
> Moti: I guess when the filter is in, you could re-run the test code
> you've mailed about, testing Ivan's Hessian operator, to verify that the
> first two eigen values are indeed equal and negative and the third eigen
> value is nearly 0 (for bright cylindrical objects in a dark background).
>
> Thanks again.
> --
> karthik
>> Regards, Dan
>>
>> 2008/10/12 Karthik Krishnan <karthik.krishnan at kitware.com>:
>>
>>> Moti:
>>>
>>> Yes, you can make the hessian filters in ITK streamable. It has been on
>>> my TODO list, since I have a need to detect needles in CT in realtime,
>>> for which I'd like to do a localized hessian analysis, which I cannot
>>> do with the existing framework.
>>>
>>> The Hessian filter in ITK, uses the RecursiveGaussianFilter to compute
>>> scale-space derivatives (see Deriche). This is fast, however its an IIR
>>> approximation of the gaussian and cannot be streamed.
>>>
>>> ITK provides a discrete gaussian operator for smoothing, (approximated
>>> from weighted sums of bessel functions). This is used in the
>>> DiscreteGaussianImageFilter. These are streamable.
>>>
>>> Ivan Macia, (see contribution on Insight-journal below [1]) has very
>>> nicely provided implementation of operators for
>>>
>>> * Local gaussian derivatives using the discrete gaussian
>>> * Local hessian using the above
>>>
>>> They are provided both in the form of operators and image functions,
>>> which you can quite trivially wrap into filters. These filters should
>>> be streamable.
>>>
>>> Unfortunately, the RecursiveGaussian has been used as the
>>> building block for most filters in ITK, so at some point, I would
>>> envision having duplicate filters using these.
>>>
>>> ------
>>> [1] http://www.insight-journal.org/browse/publication/179
>>> "Generalized Computation of Gaussian Derivatives Using ITK", Ivan
>>> Macia, VICOMTech.
>>>
>>>
>>> Thanks
>>> --
>>> Karthik Krishnan
>>> R&D Engineer,
>>> Kitware Inc.
>>> Ph: 518 371 3971 x119
>>> Fax: 518 371 3971
>>>
>>>
>>> On Sun, Oct 12, 2008 at 1:55 AM, Moti Freiman <freiman at cs.huji.ac.il>
> wrote:
>>>
>>>> Hi Luis,
>>>> Thanks for your detailed answer.
>>>> According to your answer, the Hessian filter itself is not streamable,
>>>> However, the main cause for memory insufficient is this filter,
>>>> since it produces images that are 6-times (in 3D) larger than the
> original
>>>> file.
>>>> If I'm right, streaming the EigenAnalysisFilter and the vesselness
will
> not
>>>> solve the problem of large files. ?
>>>> Moti
>>>>
>>>> On Mon, Oct 6, 2008 at 3:34 PM, Luis Ibanez <luis.ibanez at kitware.com>
> wrote:
>>>>
>>>>> Hi Moti,
>>>>>
>>>>> Looking at the code of this filter, the piece that looks suspicious is
>>>>> the internal use of the Symmetric Eigen Analysis filter.
>>>>>
>>>>> Both of these filters should be able to stream data, but, in the way
>>>>> they are nested, the EigenAnalysys filter may be trying to update the
>>>>> entire image unnecessarily.
>>>>>
>>>>> Please try the following:
>>>>>
>>>>> 1) Remove the EigenAnalysisFilter from inside the
>>>>> itkHessianSmoothed3DToVesselnessMeasureImageFilter
>>>>>
>>>>> 2) create a sequential pipeline where you connect:
>>>>>
>>>>> Hessian --> Symm.EigenAnalysis --> VesselnessFilter
>>>>>
>>>>> The Hessian filter itself is not streamable, but the
>>>>> Symm.EigenAnalysis and the Vesselness ones, are both
>>>>> pixel-wise filter, and they should be streamable.
>>>>>
>>>>>
>>>>> Please let us know if you encounter any problem.
>>>>>
>>>>> BTW: Once we sort this out, it will be great if you
>>>>> could add the conclusions of your observations
>>>>> as a review to this paper, and if you found it
>>>>> useful, we could try moving this contribution
>>>>> to the Review directory, in order to include it
>>>>> with the upcoming release of ITK 3.10.
>>>>>
>>>>>
>>>>> Regards,
>>>>>
>>>>>
>>>>> Luis
>>>>>
>>>>>
>>>>> ---------------------
>>>>> Moti Freiman wrote:
>>>>>
>>>>>> Hi,
>>>>>> I want to compute a vesselness measure on a large CT volume
> (512*512*512,
>>>>>> 16bit per voxel) using the
>>>>>> itkHessianSmoothed3DToVesselnessMeasureImageFilter
>>>>>> published in the InsightJournal
>>>>>> (http://www.insight-journal.org/browse/publication/163).
>>>>>> Since my machine memory is limited I tried to use the
>>>>>> itkStreamingImageFilter to divide the computation for several small
> parts
>>>>>> which do not require so much memory, and then merge them.
>>>>>> However it seems that the streaming filter do not divide the
> computation
>>>>>> using the itkHessianSmoothed3DToVesselnessMeasureImageFilter.
>>>>>> Any reason, idea?
>>>>>> Thanks,
>>>>>> Moti
>>>>>> --
>>>>>> __
>>>>>> Moti Freiman, Ph.D Student.
>>>>>> Medical Image Processing and Computer-Assisted Surgery Laboratory.
>>>>>> School of Computer Science and Engineering.
>>>>>> The Hebrew University of Jerusalem Givat Ram, Jerusalem 91904, Israel
>>>>>> Phone: +(972)-2-658-5371 (laboratory)
>>>>>> WWW site: http://www.cs.huji.ac.il/~freiman
>>>>>>
>>>>>>
>>>>>>
> ------------------------------------------------------------------------
>>>>>>
>>>>>> _______________________________________________
>>>>>> Insight-users mailing list
>>>>>> Insight-users at itk.org
>>>>>> http://www.itk.org/mailman/listinfo/insight-users
>>>>>>
>>>>
>>>> --
>>>> __
>>>> Moti Freiman, Ph.D Student.
>>>> Medical Image Processing and Computer-Assisted Surgery Laboratory.
>>>> School of Computer Science and Engineering.
>>>> The Hebrew University of Jerusalem Givat Ram, Jerusalem 91904, Israel
>>>> Phone: +(972)-2-658-5371 (laboratory)
>>>> WWW site: http://www.cs.huji.ac.il/~freiman
>>>>
>>>> _______________________________________________
>>>> Insight-users mailing list
>>>> Insight-users at itk.org
>>>> http://www.itk.org/mailman/listinfo/insight-users
>>>>
>>>>
>>>>
>>>
>>> --
>>> Karthik Krishnan
>>> R&D Engineer,
>>> Kitware Inc.
>>> Ph: 518 371 3971 x119
>>> Fax: 518 371 3971
>>> _______________________________________________
>>> Insight-users mailing list
>>> Insight-users at itk.org
>>> http://www.itk.org/mailman/listinfo/insight-users
>>>
>>>
>
> --
> Karthik Krishnan
> R & D Engineer,
> Kitware Inc,
> Ph: 518 371 3971 x119
> Fax: 518 371 3971
>
>
More information about the Insight-users
mailing list