[Insight-users] call metric in Observer code...bug?
Serena Fabbri
fabbri at u.washington.edu
Wed Dec 2 17:52:52 EST 2009
Hi Luis,
Thank you for your replay.
I think it was my oversight because I printed the globalMI twice......sorry about that!!
I fixed the code and now I am obtaining this error:
ExceptionObject caught !
itk::ExceptionObject (0xd01060)
Location: "typename itk::ImageToImageMetric<TFixedImage, TMovingImage>::MeasureType
itk::MattesMutualInformationImageToImageMetric<TFixedImage, TMovingImage>::GetValue(typename
itk::ImageToImageMetric<TFixedImage, TMovingImage>::ParametersType&) const [with TFixedImage =
itk::OrientedImage<float, 3>, TMovingImage = itk::OrientedImage<float, 3>]"
File: /Users/fabbri/InsightToolkit-
3.10.0/Code/Algorithms/itkMattesMutualInformationImageToImageMetric.txx
Line: 830
Description: itk::ERROR: MattesMutualInformationImageToImageMetric(0xd01b00): Too many samples
map outside moving image buffer: 0 / 1000
I have taken a look to MattesMutualInformation code and I have found that:
Initialize() calls PreComputeTransformValues()
In PreComputeTransformValues() weights and indices are calculated.
I have found that weights are 0 or e-306 and the indices are always 0.
These values are used in TransformPoint (GetValue() code) to calculate the mappedPoints.
The calculated mappedPoints are outside my movingImage.
I'd like to ask you if it is a bug or I get this error because the region I used for MI local is really little
compared to fixed and moving image.
Any suggestion will be very appreciate.
Thank you.
Serena.
On Tue, 1 Dec 2009, Luis Ibanez wrote:
> Hi Serena,
>
> Thanks for letting us know of your progress.
>
> --
>
> Regarding your question,
>
> Please note that Mutual Information is NOT and additive measure.
>
>
> That is, if you computed a value of Mutual Information, in a region R,
> and then you partition that Region R into two subregions R1 and R2,
> and compute Mutual information in R1 and then in R2, the sum of
> these two last numbers will NOT be equal to the mutual information
> of R.
>
> In short:
>
> MI( R ) != MI( R1 ) + MI( R2 )
>
> This is because MI is related to the logs of probabilities of intensity
> pairs (pixel counts).
>
> So... it still may be ok that you get similar values of Mutual Information
> from MI(R) and MI( R1). In principle, it just tells you that the statistical
> distribution of pixels intensity pairs in the Region R is similar to the one
> in the subregion R1.
>
> I don't know if that helps you with the original goal that you had for
> computing MI in subregions of the image.....
>
>
> Regards,
>
>
> Luis
>
>
> --------------------------------------------------------------------------------
---
> On Mon, Nov 30, 2009 at 8:10 PM, Serena Fabbri <fabbri at u.washington.edu> wrote:
>>
>> Hi Luis,
>>
>> thank you very much to explain clearly the mistake.
>> I fixed the code and now I am able to calculate local MI.
>>
>>
>> Now I am a bit surprised because the SetFixedImageRegion is (10,17,20)
>> pixelsize and the FixedImage is (80,80,101) pixelsize and
>> I am finding that the values of local MI are pretty close to the global
>> values of MI. I use all pixel for local MI and 10% of statistic for global
>> MI.
>>
>> Do you think this is reasonable?
>>
>> Thanks again for any suggestion.
>>
>> Serena.
>>
>>
>> These are the results of LOCAL MI:
>>
>> 0 -0.420331 1 -0.437285 2 -0.452543 3 -0.46165 4 -0.468013
>> 5 -0.473045 6 -0.476698 7 -0.479594 8 -0.482887 9
>> -0.484799 10 -0.48746 11 -0.489388 12 -0.490882 13 -0.492886
>> 14 -0.494813 15 -0.495938 16 -0.497812 17 -0.498596 18
>> -0.500178 19 -0.501496 20 -0.502828 21 -0.503783 22
>> -0.504879 23 -0.506142 24 -0.507478 25 -0.508179 26
>> -0.509647 27 -0.51022 28 -0.511936 29 -0.51237 30 -0.51411
>> 31 -0.514496 32 -0.516189 33 -0.516839 34 -0.518303 35
>> -0.518713 36 -0.520112 37 -0.520491 38 -0.521601 39
>> -0.522019 40 -0.523052 41 -0.523425 42 -0.524304 43
>> -0.524585 44 -0.525526 45 -0.52578 46 -0.526725 47 -0.526987
>> 48 -0.528026 49 -0.528595 50 -0.529578
>>
>>
>> -GLOBAL MI
>> 0 -0.420331 1 -0.437285 2 -0.452543 3 -0.46165 4 -0.468013
>> 5 -0.473045 6 -0.476698 7 -0.479594 8 -0.482887 9
>> -0.484799 10 -0.48746 11 -0.489388 12 -0.490882 13 -0.492886
>> 14 -0.494813 15 -0.495938 16 -0.497812 17 -0.498596 18
>> -0.500178 19 -0.501496 20 -0.502828 21 -0.503783 22
>> -0.504879 23 -0.506142 24 -0.507478 25 -0.508179 26
>> -0.509647 27 -0.51022 28 -0.511936 29 -0.51237 30 -0.51411
>> 31 -0.514496 32 -0.516189 33 -0.516839 34 -0.518303 35
>> -0.518713 36 -0.520112 37 -0.520491 38 -0.521601 39
>> -0.522019 40 -0.523052 41 -0.523425 42 -0.524304 43
>> -0.524585 44 -0.525526 45 -0.52578 46 -0.526725 47 -0.526987
>> 48 -0.528026 49 -0.528595
>>
>>
>> On Thu, 26 Nov 2009, Luis Ibanez wrote:
>>
>>> Serena,
>>>
>>> The "caller" argument of the Execute method is an Observer
>>> is the pointer to the object that invoked the event.
>>>
>>> In your particular case, this is the Optimizer.
>>>
>>> As Bill explained,
>>> you can only cast "caller" to the OptimizerType.
>>>
>>>
>>> For all the other elements of the registration framework, you have
>>> to pass them through other means to the command.
>>>
>>> One way of doing this is to have them as Member variables,
>>> and to have "Set" methods for them.
>>>
>>> Like
>>>
>>>> class CommandIterationUpdate : public itk::Command {
>>>> public:
>>>> typedef CommandIterationUpdate Self;
>>>> typedef itk::Command Superclass;
>>>> typedef itk::SmartPointer<Self> Pointer;
>>>> itkNewMacro( Self );
>>>>
>>>> protected:
>>>> CommandIterationUpdate() {};
>>>>
>>>> public:
>>>> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
>>>> typedef const OptimizerType * OptimizerPointer;
>>>>
>>>> typedef itk::OrientedImage< float, 3 > FixedImageType;
>>>> typedef const FixedImageType* FixedImagePointer;
>>>>
>>>> typedef itk::OrientedImage< float, 3 > MovingImageType;
>>>>
>>>> typedef itk::MattesMutualInformationImageToImageMetric<
>>>> FixedImageType, MovingImageType
>>>>>
>>>>> MetricType;
>>>>
>>>> typedef const MetricType* MetricPointer;
>>>>
>>>> typedef itk::BSplineDeformableTransform< double,3,3 >
>>>> DeformableTransformType;
>>>> typedef const DeformableTransformType* TransformPointer;
>>>>
>>> MetricPointer m_Metric;
>>> TransformPointer m_Transform;
>>> FixedImagePointer m_FixedImage;
>>> MovingImagePointer m_MovingImage;
>>>
>>>
>>> public:
>>>
>>> SetTransform( TransformType * transform ) { m_Transform = transform; }
>>> SetFixedImage( FixedImageType * fixed ) { m_FixedImage = fixed; }
>>> SetMovingImage( MovingImageType * moving ) { m_MovingImage = moving; }
>>>
>>>
>>> and then you can use them in your Execute method.
>>>
>>> Also,
>>> when you connect the Command to the optimizer,
>>> you should also do:
>>>
>>>
>>> observer->SetFixedImage( fixedImage );
>>> observer->SetMovingImage( movingImage );
>>> observer->SetTransform( transform );
>>>
>>> optimizer->AddObserver( itk::IterationEvent(), observer );
>>>
>>> In this way, the components of the registration framework
>>> will be available to the internals of your observer class.
>>>
>>>
>>> Regards,
>>>
>>>
>>> Luis
>>>
>>>
>>>
>>> -----------------------------------------------------------------------------
--
>>> On Wed, Nov 25, 2009 at 8:51 PM, Serena Fabbri <fabbri at u.washington.edu>
>>> wrote:
>>>>
>>>> Hi All,
>>>>
>>>> I am registering CT image and MRI image with Mattes MI and Bspline
>>>> Transformation.
>>>> I'd like to know the value of MI of a little area of the fixed image
>>>> during
>>>> the registration process.
>>>> I call the metric in the observer code but i get a bus error.
>>>>
>>>> could anybody suggest me where the error is?
>>>>
>>>> Thank you very much.
>>>>
>>>> Serena
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> class CommandIterationUpdate : public itk::Command {
>>>> public:
>>>> typedef CommandIterationUpdate Self;
>>>> typedef itk::Command Superclass;
>>>> typedef itk::SmartPointer<Self> Pointer;
>>>> itkNewMacro( Self );
>>>>
>>>> protected:
>>>> CommandIterationUpdate() {};
>>>>
>>>> public:
>>>> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
>>>> typedef const OptimizerType * OptimizerPointer;
>>>>
>>>> typedef itk::OrientedImage< float, 3 > FixedImageType;
>>>> typedef const FixedImageType* FixedImagePointer;
>>>>
>>>> typedef itk::OrientedImage< float, 3 > MovingImageType;
>>>>
>>>> typedef itk::MattesMutualInformationImageToImageMetric<
>>>> FixedImageType, MovingImageType
>>>>>
>>>>> MetricType;
>>>>
>>>> typedef const MetricType* MetricPointer;
>>>>
>>>> typedef itk::BSplineDeformableTransform< double,3,3 >
>>>> DeformableTransformType;
>>>> typedef const DeformableTransformType* TransformPointer;
>>>>
>>>>
>>>>
>>>>
>>>> void Execute(itk::Object *caller, const itk::EventObject & event)
>>>> {
>>>> Execute( (const itk::Object *)caller, event);
>>>> }
>>>>
>>>> void Execute(const itk::Object * object, const itk::EventObject & event)
>>>> {
>>>>
>>>> OptimizerPointer optimizer = dynamic_cast<
>>>> OptimizerPointer
>>>>>
>>>>> ( object );
>>>>
>>>> MetricPointer metric = dynamic_cast< MetricPointer >(
>>>> object
>>>> );
>>>> TransformPointer transform = dynamic_cast<
>>>> TransformPointer
>>>>>
>>>>> ( object );
>>>>
>>>> FixedImagePointer fixedImage = dynamic_cast<
>>>> FixedImagePointer>(
>>>> object );
>>>> FixedImageType::IndexType start;
>>>> FixedImageType::SizeType size;
>>>>
>>>> start[0]=43;
>>>> start[1]=30;
>>>> start[2]=0;
>>>>
>>>> size[0]=10;
>>>> size[1]=17;
>>>> size[2]=20;
>>>>
>>>>
>>>> FixedImageType::RegionType fixedRegionMetric ;
>>>> fixedRegionMetric.SetIndex(start);
>>>> fixedRegionMetric.SetSize(size);
>>>>
>>>> OptimizerType::ParametersType parameters;
>>>>
>>>>
>>>> if( !(itk::IterationEvent().CheckEvent( &event )) )
>>>> {
>>>> return;
>>>> }
>>>> metric->SetFixedImageRegion(fixedRegionMetric);
>>>> metric->SetUseAllPixels( true );
>>>> parameters = transform->GetParameters();
>>>>
>>>>
>>>> std::cout << optimizer->GetCurrentIteration() << "MIglobal "<<
>>>> optimizer->GetValue() << " MIlocal "<<metric->GetValue(parameters);
>>>> std::cout << std::endl;
>>>>
>>>>
>>>>
>>>> }
>>>> };
>>>>
>>>>
>>>> _____________________________________
>>>> Powered by www.kitware.com
>>>>
>>>> Visit other Kitware open-source projects at
>>>> http://www.kitware.com/opensource/opensource.html
>>>>
>>>> Kitware offers ITK Training Courses, for more information visit:
>>>> http://www.kitware.com/products/protraining.html
>>>>
>>>> Please keep messages on-topic and check the ITK FAQ at:
>>>> http://www.itk.org/Wiki/ITK_FAQ
>>>>
>>>> Follow this link to subscribe/unsubscribe:
>>>> http://www.itk.org/mailman/listinfo/insight-users
>>>>
>>>
>>
>>
>>
>
More information about the Insight-users
mailing list