[Insight-users] call metric in Observer code...bug?
Luis Ibanez
luis.ibanez at kitware.com
Thu Dec 3 22:08:06 EST 2009
Hi Serena,
Thanks for the details....
but... I have now more questions here
In (A), I guess that you meant that you pass to the local
BSpline transform the same AFFINE transform that you
used for the Bulk Transform of the Global BSpline Transform.
Is that right ?
--
The Error:
BSpline_MMI_RSG(4731) malloc: *** vm_allocate(size=4294926336) failed
(error code=3)
Seems to be the result of storing a negative number into an
unsigned type, in which case you get the modulo 2 of that
number under the 32 bits conversion.
The number "4294926336" is what you will get if you
assign -40960 to an unsigned type.
One of your computations may have slipped an incorrect sign....
B) When you set the optimizer for the local region you probably should
start with a very small step length.
What value of "maximum step length" are you currently using ?
C) Using All pixels in the small regions
sounds like the right thing to do in this case.
Regards,
Luis
--------------------------------------------------------------------------------------------
On Thu, Dec 3, 2009 at 2:05 PM, Serena Fabbri <fabbri at u.washington.edu> wrote:
> Hi Luis,
>
> A)
> I am using BSpline Transform.
> I pass the Trasform used for Global MI to the metric2.
>
> If I define a Transform for the small region I obtain this error:
>
> BSpline_MMI_RSG(4731) malloc: *** vm_allocate(size=4294926336) failed (error
> code=3)
> BSpline_MMI_RSG(4731) malloc: *** error: can't allocate region
> BSpline_MMI_RSG_1(4731) malloc: *** set a breakpoint in szone_error to debug
> Abort trap
>
>
>
> This is the code for the Transform
>
> FixedImageType::IndexType start;
> FixedImageType::SizeType size;
>
> start[0]=43;
> start[1]=17;
> start[2]=0;
>
> size[0]=10;
> size[1]=13;
> size[2]=20;
>
>
> FixedImageType::RegionType fixedRegionMetric;
> fixedRegionMetric.SetIndex(start);
> fixedRegionMetric.SetSize(size);
>
>
> typedef DeformableTransformType::RegionType RegionTypeM;
> RegionTypeM bsplineRegionM;
> RegionTypeM::SizeType gridSizeOnImageM;
> RegionTypeM::SizeType gridBorderSizeM;
> RegionTypeM::SizeType totalGridSizeM;
>
> unsigned int numberOfGridNodesInOneDimensionCoarseM = 3;
> numberOfGridNodesInOneDimensionCoarseM = 10;
> gridSizeOnImageM.Fill( numberOfGridNodesInOneDimensionCoarseM );
> gridBorderSizeM.Fill( SplineOrder ); // Border for spline order = 3 ( 1
> lower, 2 upper )
> totalGridSizeM = gridSizeOnImageM + gridBorderSizeM;
>
> bsplineRegionM.SetSize( totalGridSizeM );
>
>
> SpacingType spacingM = fixedImage->GetSpacing();
> OriginType originM = fixedImage->GetOrigin();
>
> FixedImageType::SizeType fixedImageSizeM = fixedRegionMetric.GetSize();
> std::cout<<fixedRegionMetric.GetSize()<<std::endl;
>
> for(unsigned int r=0; r<ImageDimension; r++)
> {
> spacingM[r] *= static_cast<double>(fixedImageSizeM[r] - 1) /
> static_cast<double>(gridSizeOnImageM[r] - 1);
> }
>
> FixedImageType::DirectionType gridDirectionM = fixedImage->GetDirection();
> SpacingType gridOriginOffsetM = gridDirectionM * spacingM;
>
> OriginType gridOriginM = originM - gridOriginOffsetM;
>
> bsplineTransformCoarseM->SetGridSpacing( spacingM );
> bsplineTransformCoarseM->SetGridOrigin( gridOriginM );
> bsplineTransformCoarseM->SetGridRegion( bsplineRegionM );
> bsplineTransformCoarseM->SetGridDirection( gridDirectionM );
> bsplineTransformCoarseM->SetParameters(
> initialDeformableTransformParameters );
>
>
>
>
> B)
> local MI uses the same optimizer than global MI (Regular Step Gradient
> Descent).
> It works fine in the registration task.
>
>
> C)
> I used SetUseAllPixels for the small region too. The region is small in fact
> the number of pixel is 2600.
>
>
>
>
>
>
>
> On Wed, 2 Dec 2009, Luis Ibanez wrote:
>
>> 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