[Insight-users] Fwd: Bug in itk::kdTree ?

Bill Lorensen bill.lorensen at gmail.com
Tue Dec 1 19:40:33 EST 2009 

Luis,

I think the problem has been isolated to the fact that itkKdTree is
not thread-safe. Unfortunately, this issue has been active in multiple
threads.

Motes can verify?

Bill

On Tue, Dec 1, 2009 at 6:02 PM, Luis Ibanez <luis.ibanez at kitware.com> wrote:
> ---------- 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
>>>>>>
>>>>>
>>>>
>>>
>>
> _____________________________________
> 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