[Insight-users] call metric in Observer code...bug?

Luis Ibanez luis.ibanez at kitware.com
Wed Dec 2 20:12:05 EST 2009


Hi Serena,


Thanks for the update.


The new error message indicates that the number of samples
landing  in the overlapping region between the FixedImageRegion
and the moving image is too small.

This is usually an  indication of:


            A)   A poorly initialized Transform,   or

            B)   A transform that went too far
                  usually due to a misconfigured optimizer

            C)   Setting the number of samples
                  to a value that is too small.


Let's start with (A),.... how are you initializing
the Transform used to compute the value of MI
in the local small regions ?


     Please let us know,


          Thanks


                Luis


-----------------------------------------------------
On Wed, Dec 2, 2009 at 5:52 PM, Serena Fabbri <fabbri at u.washington.edu> wrote:
>
> 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