[Insight-users] Re: Antw: (Addition) MattesMutualInformation Registration

Karthik Krishnan Karthik.Krishnan at kitware.com
Tue Jun 7 18:41:21 EDT 2005


Jan,

Thanks for reporting your findings

It helped us find a bug in 
itkMattesMutualInformationImageToImageMetric.txx.   
SampleFullFixedImageDomain as you thought you wasn't really sampling the 
entire fixed image domain. The bug should be fixed in a few minutes. 
This was the reason for the randomness with UseAllPixels and for your 
occassional crashes. Please update 
itkMattesMutualInformationImageToImageMetric.txx

On another note, the parameters in your code were not really 
appropriate. A maximum step size of 20 is too high, especially when you 
are using Multi-resolution registration. You want to make it say 3 and 
double it at every resolution level. [ my guess is that this drove the 
registration away ]. You were also using way too many bins.

Please find attached your code, with the parameters changed. You should 
get a reasonable registration with the following command line.
RegistrationDebug BrainT1SliceBorder20.png 
BrainProtonDensitySliceShifted13x17y.png

Thanks
Regards
Karthik

Jan Schreiber wrote:

>Hi Karthik, 
>
>could it be that there is a conflict with 
>    SetNumberOfSpatialSamples( ... ) and SetUseAllPixels( ... )
>in MattesMutualInformation?
>
>If I set SetUseAllPixels( true ), the computation time still depends on the number of spatial samples - which surprises me, becaus I thought it would just iterate though all the pixels in the fixed image region and hence should only depend on the region's size. 
>
>Furthermore I debugged MattesMutualInformationImageToImageMetric::SampleFullFixedImageDomain() and it showed that the variable regionIter reaches its end (regionIter.IsAtEnd()) before (!) being incremented. 
>
>The risk of a crashes increases with a higher number of spatial samples. 
>
>I hope, this information is valuable for you.
>
>Thanks and best regards, 
>Jan
>
>P.S.: also strange is, that the registration results are not allways the same when running the application several times. I thought with setting to use all pixels, the randomness should be gone...
>
>P.P.S.: RegistrationDebug image1 image2
>
>
>  
>
>>>>Karthik Krishnan <Karthik.Krishnan at kitware.com> 06/06/05 5:16  >>>
>>>>        
>>>>
>Oh, I'm sorry. I was assuming that you were using the class 
>ImageRegistrationMethod and manually resampling and registering 
>different resolution levels. You are right. You needn't do this stuff if 
>you are using MutiResolutionImageRegistrationMethod cause the pyramid 
>helper class recomputes the fixed image regions for you.
>
>That still doesn't explain the crash from the pixel accessor. You may 
>have to post a minimal compiling code so we can look at it.
>
>Thanks
>regards
>karthik
>
>Jan Schreiber wrote:
>
>  
>
>>Hi Karthik, 
>>
>>thank you very much for your detailed advise. 
>>But, are you sure that I have to re-compute the FixedImageRegions for every level?
>>I thought, this is automatically done in 
>>
>>MultiResolutionImageRegistrationMethod<TFixedImage,TMovingImage>::PreparePyramids( void )
>>
>>If I am wrong: where do I have to implement the code below, into the IterationEventObserver? Where do I get the variable idx[] from?
>>
>>Thank you very much for your help! 
>>Jan
>>
>>
>> 
>>
>>    
>>
>>>>>Karthik Krishnan <Karthik.Krishnan at kitware.com> 06/03/05 8:17  >>>
>>>>>       
>>>>>
>>>>>          
>>>>>
>>In case I wasn't clear enough:
>>
>>You need to be doing this at every resolution level. So you would have 
>>something like
>>
>> for(unsigned int resolution=0; resolution < 
>>numberOfResolutionLevelsToUse; resolution++)
>>   {
>>   RegisterResolutionLevel(); 
>>   }
>>
>>RegisterResolutionLevel()
>>{
>>fixedImageRegionStart[0] = static_cast<int>( idx[0] / factor );
>>fixedImageRegionStart[1] = static_cast<int>( idx[2] / factor );
>>fixedImageRegionStart[2] = static_cast<int>( idx[4] / factor );
>>
>>fixedImageRegionSize[0] =  static_cast<int>( (idx[1] - idx[0] + 1 ) / 
>>factor );
>>fixedImageRegionSize[1] =  static_cast<int>( (idx[3] - idx[2] + 1 ) / 
>>factor );
>>fixedImageRegionSize[2] =  static_cast<int>( (idx[5] - idx[4] + 1 ) / 
>>factor );
>>
>>fixedImageRegion.SetIndex( fixedImageRegionStart );
>>fixedImageRegion.SetSize( fixedImageRegionSize );
>>
>>m_RegistrationMethod->SetFixedImageRegion( fixedImageRegion ); 
>>m_RegistrationMethod->StartRegistration()
>>}
>>
>>
>>Karthik Krishnan wrote:
>>
>> 
>>
>>    
>>
>>>Jan,
>>>
>>>If you are using Multi-resolution registration and you are restricting 
>>>the fixed image region, you need to shrink it down yourself. The 
>>>ImageRegistrationMethod class will not do this for you. In other 
>>>words, you should be doing:
>>>
>>>fixedImageRegionStart[0] = static_cast<int>( idx[0] / factor );
>>>fixedImageRegionStart[1] = static_cast<int>( idx[2] / factor );
>>>fixedImageRegionStart[2] = static_cast<int>( idx[4] / factor );
>>>
>>>fixedImageRegionSize[0] =  static_cast<int>( (idx[1] - idx[0] + 1 ) / 
>>>factor );
>>>fixedImageRegionSize[1] =  static_cast<int>( (idx[3] - idx[2] + 1 ) / 
>>>factor );
>>>fixedImageRegionSize[2] =  static_cast<int>( (idx[5] - idx[4] + 1 ) / 
>>>factor );
>>>
>>>fixedImageRegion.SetIndex( fixedImageRegionStart );
>>>fixedImageRegion.SetSize( fixedImageRegionSize );
>>>
>>>m_RegistrationMethod->SetFixedImageRegion( fixedImageRegion );
>>>
>>>
>>>HTH
>>>karthik
>>>
>>>
>>>Jan Schreiber wrote:
>>>
>>>   
>>>
>>>      
>>>
>>>>Hi Karthik,
>>>>unfortunately, an update of the cvs sources did not help.
>>>>I am not using ImageAdaptors directly and a breakpoint in the 
>>>>ImageAdaptor's constructor has never been reached. So it seams also 
>>>>not to be used indirectly.
>>>>It seams that the m_Position variable gets an invalid value within 
>>>>the ImageRegionConstIteratorWithIndex::opearator++(). At least the 
>>>>instruction  InternalPixelType ipt = *m_Position;
>>>>at the beginning of the function does not make any problems,  ipt = 
>>>>*m_Position;
>>>>at the end of the function throws an exception...
>>>>
>>>>Another thing I do not really understand is the ShrinkImageFilter 
>>>>(used by the pyramid). As the shrinkfactor is 8, the pixel spacing of 
>>>>my shrinked image is also 8.0 (having a spacing of 1.0 in the 
>>>>original image). But as my image has width and height that is not 
>>>>divisible by 8, the ( shrinked image with ) * ( shrinkfactor ) != ( 
>>>>original image width)
>>>>Could that matter?
>>>>
>>>>Thanks and regards, Jan
>>>>
>>>>P.S.: thank you also for your replay on the demons questions
>>>>
>>>>
>>>>
>>>>
>>>>     
>>>>
>>>>        
>>>>
>>>>>>>Karthik Krishnan <Karthik.Krishnan at kitware.com> 06/02/05 7:22  >>>
>>>>>>>     
>>>>>>>           
>>>>>>>
>>>>>>>              
>>>>>>>
>>>>Ah,
>>>>
>>>>We ran into nearly the same problem yesterday. Are you using 
>>>>ImageAdaptors in your application. There was a problem with the 
>>>>adaptor update mechanism that caused the application to crash cause 
>>>>the m_PixelAccessor found the buffered region to be empty.
>>>>
>>>>Try updating itkImageBase.h and itkImageBase.txx and 
>>>>itkImageAdaptor.cxx all in Code/Common
>>>>
>>>>Thanks
>>>>regards
>>>>karthik
>>>>
>>>>Jan Schreiber wrote:
>>>>
>>>>
>>>>
>>>>     
>>>>
>>>>        
>>>>
>>>>>Hi,
>>>>>my registation of two x-ray images (sometimes) crashes.
>>>>>I am using the cvs version (about 1-2 weeks old) of ITK, Microsoft 
>>>>>Visual C++ .NET and the following components:
>>>>>itk::RegularStepGradientDescentOptimizer
>>>>>itk::LinearInterpolateImageFunction
>>>>>itk::MultiResolutionImageRegistrationMethod
>>>>>itk::(Recursive)MultiResolutionPyramidImageFilter
>>>>>itk::(Recursive)MultiResolutionPyramidImageFilter
>>>>>itk::MattesMutualInformationImageToImageMetric
>>>>>
>>>>>
>>>>>The Error occures in ImageConstIteratorWithIndex.h line 234
>>>>>PixelType Get(void) const    { return 
>>>>>m_PixelAccessor.Get(*m_Position); }
>>>>>because m_Position is not valid.
>>>>>This Get() Function is called in line 532 of 
>>>>>itkMattesMutualInformationImageToImageMetric.txx
>>>>>   (*iter).FixedImageValue = regionIter.Get();
>>>>>
>>>>>Consulting the debugger gives information shown at the end of this 
>>>>>mail ( have a look at the line with <<========== )
>>>>>
>>>>>The Spacing of my image is usually 1.0 (as I am reading PNG Images). 
>>>>>Within the pyramid (the output below is from level 4) the spacing is 
>>>>>8.0. Is this realted to the downsampling?
>>>>>Originally, my Images have a size of 2828x2320 Pixels, 8 bit stored 
>>>>>as float and the FixedImageRegion ist set to x:1800, y:1000, w:940, 
>>>>>h:1290.
>>>>>
>>>>>m_Position gets invalid if the level of the 
>>>>>MultiResolutionImageRegistrationMethod is >= 3.
>>>>>The chance for this crah is almost 100% when SetUseAllPixels(true) 
>>>>>is used and the FixedImageRegion is not the LergestPossibleRegion. 
>>>>>Taking random Pixels reduces the risk of a crah but does not 
>>>>>eliminate it. Setting the FixedImageRegion to something smaller than 
>>>>>the LargestPossibleRegion increases the crash-frequency .
>>>>>
>>>>>I would be very glad if you had any suggestion...
>>>>>
>>>>>Thanks in advance, Jan
>>>>>
>>>>>
>>>>>
>>>>>Debugger output:
>>>>>
>>>>>-    regionIter    {...}    
>>>>>itk::ImageRegionConstIteratorWithIndex<ibia::Image2DFloat>
>>>>>-    itk::ImageConstIteratorWithIndex<ibia::Image2DFloat>    
>>>>>{ImageDimension=??? m_Image={m_Pointer=0x0a088e98 } 
>>>>>m_PositionIndex={m_Index=0x011edc1c } ...}    
>>>>>itk::ImageConstIteratorWithIndex<ibia::Image2DFloat>
>>>>>-    m_Image    {m_Pointer=0x0a088e98 }    
>>>>>itk::WeakPointer<ibia::Image2DFloat const >
>>>>>-    m_Pointer    0x0a088e98    const ibia::Image2DFloat *
>>>>>-    ibia::Image2D<float>    {...}    const ibia::Image2D<float>
>>>>>-    ibia::ImageBase<float,2>    {...}    const 
>>>>>ibia::ImageBase<float,2>
>>>>>-    itk::Image<float,2>    {ImageDimension=??? }    const 
>>>>>itk::Image<float,2>
>>>>>-    itk::ImageBase<2>    {ImageDimension=??? m_Spacing={...} 
>>>>>m_Origin={...} ...}    const itk::ImageBase<2>
>>>>>+    itk::DataObject    {m_Source={m_Pointer=0x0a088db8 
>>>>>{m_Updating=??? m_OutputInformationMTime=??? m_Inputs=??? ...} } 
>>>>>m_SourceOutputIndex=0 m_UpdateMTime={m_ModifiedTime=2988 } ...}    
>>>>>itk::DataObject
>>>>>-    m_Spacing    {...}    itk::Vector<double,2>
>>>>>-    itk::FixedArray<double,2>    {Length=??? }    
>>>>>itk::FixedArray<double,2>
>>>>>-    m_InternalArray    0x0a088ee0    double [2]
>>>>>   [0]    8.0000000000000000    double
>>>>>   [1]    8.0000000000000000    double
>>>>>-    m_Origin    {...}    itk::Point<double,2>
>>>>>-    itk::FixedArray<double,2>    {Length=??? }    
>>>>>itk::FixedArray<double,2>
>>>>>-    m_InternalArray    0x0a088ef0    double [2]
>>>>>   [0]    0.00000000000000000    double
>>>>>   [1]    0.00000000000000000    double
>>>>>+    m_Direction    {RowDimensions=??? }    itk::Matrix<double,2,2>
>>>>>+    m_OffsetTable    0x0a088f20    long [3]
>>>>>-    m_LargestPossibleRegion    {ImageDimension=??? 
>>>>>SliceDimension=??? }    itk::ImageRegion<2>
>>>>>+    itk::Region    {...}    const itk::Region
>>>>>-    m_Index    {m_Index=0x0a088f30 }    itk::Index<2>
>>>>>-    m_Index    0x0a088f30    long [2]
>>>>>   [0]    0    long
>>>>>   [1]    0    long
>>>>>-    m_Size    {m_Size=0x0a088f38 }    itk::Size<2>
>>>>>-    m_Size    0x0a088f38    unsigned long [2]
>>>>>   [0]    353    unsigned long
>>>>>   [1]    290    unsigned long
>>>>>-    m_RequestedRegion    {ImageDimension=??? SliceDimension=??? 
>>>>>}    itk::ImageRegion<2>
>>>>>+    itk::Region    {...}    const itk::Region
>>>>>-    m_Index    {m_Index=0x0a088f44 }    itk::Index<2>
>>>>>-    m_Index    0x0a088f44    long [2]
>>>>>   [0]    0    long
>>>>>   [1]    0    long
>>>>>-    m_Size    {m_Size=0x0a088f4c }    itk::Size<2>
>>>>>-    m_Size    0x0a088f4c    unsigned long [2]
>>>>>   [0]    353    unsigned long
>>>>>   [1]    290    unsigned long
>>>>>-    m_BufferedRegion    {ImageDimension=??? SliceDimension=??? }    
>>>>>itk::ImageRegion<2>
>>>>>+    itk::Region    {...}    const itk::Region
>>>>>-    m_Index    {m_Index=0x0a088f58 }    itk::Index<2>
>>>>>-    m_Index    0x0a088f58    long [2]
>>>>>   [0]    0    long
>>>>>   [1]    0    long
>>>>>-    m_Size    {m_Size=0x0a088f60 }    itk::Size<2>
>>>>>-    m_Size    0x0a088f60    unsigned long [2]
>>>>>   [0]    353    unsigned long
>>>>>   [1]    290    unsigned long
>>>>>+    m_Buffer    {m_Pointer=0x0a088fb0 {m_ImportPointer=??? 
>>>>>m_Size=??? m_Capacity=??? ...} }    
>>>>>itk::SmartPointer<itk::ImportImageContainer<unsigned long,float> >
>>>>>-    m_PositionIndex    {m_Index=0x011edc1c }    itk::Index<2>
>>>>>-    m_Index    0x011edc1c    long [2]
>>>>>   [0]    257    long
>>>>>   [1]    129    long
>>>>>-    m_BeginIndex    {m_Index=0x011edc24 }    itk::Index<2>
>>>>>-    m_Index    0x011edc24    long [2]
>>>>>   [0]    225    long
>>>>>   [1]    125    long
>>>>>-    m_EndIndex    {m_Index=0x011edc2c }    itk::Index<2>
>>>>>-    m_Index    0x011edc2c    long [2]
>>>>>   [0]    342    long
>>>>>   [1]    286    long
>>>>>-    m_Region    {ImageDimension=??? SliceDimension=??? }    
>>>>>itk::ImageRegion<2>
>>>>>+    itk::Region    {...}    const itk::Region
>>>>>+    m_Index    {m_Index=0x011edc38 }    itk::Index<2>
>>>>>+    m_Size    {m_Size=0x011edc40 }    itk::Size<2>
>>>>>-    m_OffsetTable    0x011edc48    unsigned long [3]
>>>>>   [0]    1    unsigned long
>>>>>   [1]    353    unsigned long
>>>>>   [2]    102370    unsigned long
>>>>>-    m_Position    0x0a178000    const float *
>>>>>       CXX0030: Fehler: Ausdruck kann nicht ausgewertet werden    
>>>>>const float       <<================
>>>>>-    m_Begin    0x0a13f520    const float *
>>>>>       86.568977    const float
>>>>>-    m_End    0x0a176970    const float *
>>>>>       8.7256823    const float
>>>>>   m_Remaining    true    bool
>>>>>   m_PixelAccessor    {...}    itk::DefaultPixelAccessor<float>
>>>>>
>>>>>_______________________________________________
>>>>>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 
>>>
>>>   
>>>
>>>      
>>>
>> 
>>    
>>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: main.cpp
Type: text/x-c++src
Size: 12972 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20050607/c0313610/main-0001.cpp


More information about the Insight-users mailing list