[Insight-users] call metric in Observer code...bug?
Serena Fabbri
fabbri at u.washington.edu
Thu Dec 3 14:05:53 EST 2009
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