[Insight-users] Fwd: Bug in itk::kdTree ?
Luis Ibanez
luis.ibanez at kitware.com
Tue Dec 1 18:02:25 EST 2009
---------- Forwarded message ----------
From: Luis Ibanez <luis.ibanez at kitware.com>
Date: Tue, Dec 1, 2009 at 6:02 PM
Subject: Re: [Insight-users] Bug in itk::kdTree ?
To: motes motes <mort.motes at gmail.com>
Hi Motes,
Here is a suggestion for tracking this problem.
If we take the test:
Insight/Testing/Code/Numerics/Statistics/
itkKdTreeGeneratorTest.cxx
and replace the initialization of the array in lines: 39-45
with loading the list of points that you are generating
for your test,
then we should be able to reproduce the problem
that you are seeing.
Could you describe how are you generating the set
of points ?
Or send us a file with the list of points ?
(assuming that is less than 1Mb in size...)
In this way we can track is this is actually a bug in
the KdTree (or its generator)
Thanks
Luis
------------------------------------------------------------------------------
On Thu, Nov 26, 2009 at 2:59 PM, motes motes <mort.motes at gmail.com> wrote:
> I have written my own transform where no of the regular grid
> structures (bspline grid, coeff_image, etc) are used.
>
> Since the transform gets points in physical coordinates from the
> metric the radius and other measures used in the transform must also
> be in converted to physical coordinates.
>
>
>
>
>
> Here is an example:
>
> I define the radius as the basisfunction doman divided by 2. I still
> use the BSpline kernel which has a domain over 4 "blocks" and the
> radius is therefore 2 "blocks".
>
>
>
> One "block" in pixel coordinates is:
>
> blockInPixels = ImageDimension/(numberOfNodes-1)
>
> So if we have a 128*128*128 volume and define 4*4*4 nodes each block is:
>
> blockInPixel = 128/3 = 43 pixels
>
> Since the radius is two blocks the radius becomes 2*43 = 85 pixels.
>
> Now we need to convert this to physical space which is done with the formula:
>
> physical = I*S+O
>
> where I is index, S is scaling and O is origin of the fixed image.
> Assume that the fixed image has the following spacing and origin
> respectivily:
>
> [3.125, 3.125, 4] and [0,0,0]
>
> the radius in each dimension therefore becomes:
>
> [3.125 * 85, 3.125 * 85, 4 * 85] = [266, 266, 340]
>
>
> Currently I choose the first element as radius, 266 for all dimensions
> and feed that to the kd-tree.
>
> But maybe I am doing something wrong in this custom transform?
>
> But in principle it should be possible to specify a radius of any size
> to the kd-tree without getting a segmentation fault.
>
>
>
>
>
>
>
>
>
>
>
>
>
> On Thu, Nov 26, 2009 at 7:59 PM, Luis Ibanez <luis.ibanez at kitware.com> wrote:
>> Hi Motes,
>>
>> Why are you using such a large radius ?
>>
>> The method will return ALL the points that are inside
>> a hyper-shpere of radius 262.
>>
>> Do you really anticipate that to be the region of the
>> space that you need is that large ?
>>
>>
>> Luis
>>
>>
>> --------------------------------------------------------
>> On Thu, Nov 26, 2009 at 1:55 PM, motes motes <mort.motes at gmail.com> wrote:
>>> The value is:
>>>
>>> m_radius = 262
>>>
>>> when I do std::cout << "m_radius" << std::endl; I have also tried:
>>>
>>> static_cast<double>(m_kernel_radius)/1.0;
>>>
>>> But I still get the error.
>>>
>>>
>>>
>>>
>>> When I look in the implementation (itkKdTree.txx) in the ::SearchLoop
>>> function (a recursive function) I see that the radius is used:
>>>
>>>
>>> // search the other node, if necessary
>>> tempValue = lowerBound[partitionDimension];
>>> lowerBound[partitionDimension] = partitionValue;
>>> if ( this->BoundsOverlapBall(query, lowerBound, upperBound,
>>> m_SearchRadius) )
>>> {
>>> SearchLoop(node->Right(), query, lowerBound, upperBound);
>>> }
>>> lowerBound[partitionDimension] = tempValue;
>>>
>>>
>>> My thought was therefore (under the assumption that a "bad" radius is
>>> supplied) that the function never terminates since its a recursive
>>> function.
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>> On Thu, Nov 26, 2009 at 7:46 PM, Luis Ibanez <luis.ibanez at kitware.com> wrote:
>>>> Hi Motes,
>>>>
>>>> What is the value of m_radius ?
>>>>
>>>> Please report this to the mailing list.
>>>>
>>>>
>>>> Thanks
>>>>
>>>>
>>>> Luis
>>>>
>>>>
>>>> -----------------------------------------
>>>> On Thu, Nov 26, 2009 at 1:32 PM, motes motes <mort.motes at gmail.com> wrote:
>>>>> For the last couple of days I have tried to find a solution for a
>>>>> segmentation fault that happens when I use the itk::kdTree.
>>>>>
>>>>> I am using the version of search where the radius is specified:
>>>>>
>>>>> /** Searches the neighbors fallen into a hypersphere */
>>>>> void Search(const MeasurementVectorType &query,
>>>>> double radius,
>>>>> InstanceIdentifierVectorType& result) const;
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> From my function I call it with:
>>>>>
>>>>> VectorType vectorPoint;
>>>>> for (int i=0; i<NDimensions; i++) {
>>>>> vectorPoint[i] = point[i];
>>>>> }
>>>>>
>>>>> NeighborsType neighbors;
>>>>> tree->Search(vectorPoint, m_radius, neighbors);
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> Where radius is computed like:
>>>>>
>>>>> unsigned int FVOrder = 3;
>>>>> double Order = (double)(FVOrder+1.0);
>>>>> double m_radius = (1.0*maxValue * Order)/2.0;
>>>>>
>>>>>
>>>>> Now if I do:
>>>>>
>>>>> tree->Search(vectorPoint, 5.0, neighbors);
>>>>>
>>>>> it works fine without any segmentation errors. But if I do:
>>>>>
>>>>> tree->Search(vectorPoint, m_radius, neighbors);
>>>>>
>>>>> where m_radius is compted as above I get the segmentation error (after
>>>>> a few thousand calls).
>>>>>
>>>>> Could this be a bug or am I doing something wrong here?
>>>>> _____________________________________
>>>>> 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