[Insight-users] MattesMutualInformation Registration
UseAllPixels() -> solved!
Jan Schreiber
Jan.Schreiber at umit.at
Wed Jun 8 09:25:58 EDT 2005
Hi Karthik,
thank you very much for your help and the fast solving of the bug - that is really phantastic!
My registration now works without any problems, even with the high step size that is needed for the 2828x2320 x-ray images. ;-)
Thanks and best regards,
Jan
>>> Karthik Krishnan <Karthik.Krishnan at kitware.com> 06/08/05 12:41 >>>
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
>>>
>>>
>>>
>>>
>>>
>>
>>
>>
More information about the Insight-users
mailing list