BUG: Lack of precision in ImageData finding procedures

[David Gobbi]

Hi Christophe,

I see what you mean, but I think there is a fix that doesn't require
floating-point comparisons, only a check whether  pcoords[i] == 0.

This should do the trick:

//--------------------------------------------------------------------
// Convenience function computes the structured coordinates for
// a point x[3].  The voxel is specified by the array ijk[3],
// and the parametric coordinates in the cell are specified with
// pcoords[3]. The function returns a 0 if the point x is outside
// of the volume, and a 1 if inside the volume.
int vtkImageData::ComputeStructuredCoordinates(float x[3], int ijk[3], 
					       float pcoords[3])
{
  int i;
  float d, floatLoc;
  float *origin = this->GetOrigin();
  float *spacing = this->GetSpacing();
  
  //
  //  Compute the ijk location
  //
  for (i=0; i<3; i++) 
    {
    d = x[i] - origin[i];
    floatLoc = d / spacing[i];
    ijk[i] = (int)floor(floatLoc);
    pcoords[i] = floatLoc-ijk[i];

    if ( ijk[i] >= this->Extent[i*2] && ijk[i] < this->Extent[i*2+1] )
      {
      continue;
      }

    else if ( ijk[i] < this->Extent[i*2] || ijk[i] > this->Extent[i*2+1] 
	      || pcoords[i] != 0) 
      {
      return 0;
      } 
 
    else if (this->Dimensions[i] != 1) //&& (ijk[i] == this->Extent[i*2+1]
      {                                //    && pcoords[i] == 0)
      ijk[i] -= 1;
      pcoords[i] = 1.0;
      }
    }
  return 1;
}

 - David

--
  David Gobbi, MSc                    dgobbi at irus.rri.on.ca
  Advanced Imaging Research Group
  Robarts Research Institute, University of Western Ontario

On Tue, 23 May 2000, Odet Christophe wrote:

> Hi David,
> 
> Thanks for your help.
> 
> I agree with you about the floating overhead but I think that your solution does not solve 
> completely the problem (may be I don't understand your solution !)
> 
> in vtkImageData::ComputeStructuredCoordinates the value -0.5 is converted to 0. Using a floor
> statement (or vtkFloor) will correctly convert it to -1 and this effectively solves the problem for this
> kind of values. Nevertheless, considering the case where your extent is (0,1), a value of 1.5 is converted to
> 1 (with or without the floor function) and the lines
>  
> else // if ( ijk[i] == this->Extent[i*2+1] )
>       {
>       if (this->Dimensions[i] == 1)
>         {
>         pcoords[i] = 0.0;
>         }
>       else
>         {
>         ijk[i] -= 1;
>         pcoords[i] = 1.0;
>         }
>       }
> 
> will consider the 1.5 point to belong to the extent. This is false and the extent is artificially extended.
> if the value x is exactly equal to 1 the fucntion is ok as it interpolates inside the 0,1 cell with pcoords=1.0.
> But if x is greater than 1.0 (i.e 1.2) then answer is wrong giving a pcoords of 1.0 . 
> 
> So this is not only a problem of floor'ed int values.
> 
> The vtkImageData::ComputeStructuredCoordinates needs to be modified in a much more complicated
> way. 
> 
> Can we really do something without some foat comparisons ?
> May be with comparison based on some floor'ed and ceil'ed values....
> Do you agree with that ?
> 
> Christophe
> 
> --------------------------------------------------------------------
> This is the private VTK discussion list. Please keep messages on-topic.
> Check the FAQ at: <http://public.kitware.com/cgi-bin/vtkfaq>
> To UNSUBSCRIBE, send message body containing "unsubscribe vtkusers" to
> <majordomo at public.kitware.com>. For help, send message body containing
> "info vtkusers" to the same address.
> --------------------------------------------------------------------
> 

--------------------------------------------------------------------
This is the private VTK discussion list. Please keep messages on-topic.
Check the FAQ at: <http://public.kitware.com/cgi-bin/vtkfaq>
To UNSUBSCRIBE, send message body containing "unsubscribe vtkusers" to
<majordomo at public.kitware.com>. For help, send message body containing
"info vtkusers" to the same address.
--------------------------------------------------------------------