Lack of precision ImageData finding procedures

[David Gobbi]

Hi Christophe,

One problem with your proposed code is that the compiler must
convert all the extents to float before doing the comparisons,
which will be quite a performance hit.

A better solution would be to change

 ijk[i] = (int) floatLoc;

to

 ijk[i] = (int)floor(floatLoc);


For my own VTK classes, I have a simple macro to replace the
expensive call to floor:

static inline int vtkFloor(float x)
{
  int ix = int(x);
  if (x-ix < 0) { ix--; }

  return ix;
}

And another macro specifically for parametric coords:

static inline int vtkFloor(float x, float &p)
{
  int ix = int(x);
  p = x-ix;
  if (p < 0) { p = x - (--ix); }

  return ix;
}

which is used as follows:

  ijk[i] = vtkFloor(floatLoc,pcoords[i])

This guarantees round-to-neg-infinity behaviour, and the
parametric coords are always in the range [0,1).

 - David

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


On Fri, 19 May 2000, Odet Christophe wrote:

> I think that there is a problem with the function vtkImageData::ComputeStructuredCoordinates
> because the 
> 
> ijk[i] = (int) floatLoc 
> 
> statement rounds towards zero the floatLoc value. When the floatLoc equals -0.5, for example,
> it is rounded to 0 and the x point seems to belongs to the cell beginning at point 0,0,0. The 
> problem arises when dealing with x values greater or lower than than the extent.
> 
> I suggest to change some reference to ijk[i] with floatLoc in the code lines:
> 
> ######### OLD CODE
>         if ( ijk[i] >= this->Extent[i*2] && ijk[i] < this->Extent[i*2 + 1] )
>       {
>       pcoords[i] = floatLoc - (float)ijk[i];
>       }
> 
>     else if ( ijk[i] < this->Extent[i*2] || ijk[i] > this->Extent[i*2+1] ) 
>       {
>       return 0;
>       } 
> ########## END OF OLD CODE
> 
> to obtain 
> 
> ######### NEW PROPOSED CODE
>         if ( floatLoc >= this->Extent[i*2] && floatLoc < this->Extent[i*2 + 1] )
>       {
>       pcoords[i] = floatLoc - (float)ijk[i];
>       }
> 
>     else if ( floatLoc < this->Extent[i*2] || floatLoc > this->Extent[i*2+1] ) 
>       {
>       return 0;
>       } 
> ########## END OF NEW CODE
> 
> It test the patch with the following little tcl/tk script and it seems to work well. The dataset implicit 
> clipping is more precise with the new code (this is visible with small objects as used in the script).
> In the original version a too large piece of the sphere is selected. In the new version, only some
> part of the cells clipped by the dataset are selected. I hope this will help, but I think that the same problem
> could appear in some other functions of the VTK library (FindCells, FindPoints....).
> 
> ##########
> 
> # Object to be cut. A simple little sphere
> 
> vtkSphereSource sphe
> sphe SetRadius 2
> 
> vtkPolyDataMapper spheMapper
>     spheMapper SetInput [sphe GetOutput]
>     spheMapper ScalarVisibilityOff
> vtkActor spheActor
>     spheActor SetMapper spheMapper
> 
> # The cutter... a simple parallelogram
> 
> vtkStructuredPoints sp
> sp SetDimensions 2 2 2
> sp SetSpacing 5 0.1 2
> sp AllocateScalars
> set spsc [[sp GetPointData] GetScalars]
> for {set i 0} {$i<8} { incr i } {
>  $spsc SetScalar $i 1
> }
> vtkDataSetMapper spMapper
>   spMapper SetInput sp
> vtkActor spActor
>   spActor SetMapper spMapper 
> 
> vtkImplicitVolume spds
> spds SetVolume sp
> spds SetOutValue 0
> spds SetOutGradient 1 0 0
> 
> # Cutting...
> 
> vtkClipPolyData clipper
>  clipper SetInput [sphe GetOutput]
>  clipper SetClipFunction spds
>  clipper GenerateClipScalarsOn
>  clipper SetValue 0.5
> 
> vtkPolyDataMapper clipMapper
>     clipMapper SetInput [clipper GetOutput]
> vtkActor clipActor
>     clipActor SetMapper clipMapper
> 
> 
> # Create graphics stuff
> #
> vtkRenderer ren1
> vtkRenderWindow renWin
>     renWin AddRenderer ren1
> vtkRenderWindowInteractor iren
>     iren SetRenderWindow renWin
> 
> # Add the actors to the renderer, set the background and size
> ren1 AddActor spheActor
> [spheActor GetProperty] SetRepresentationToWireframe
> ren1 AddActor spActor 
> #spActor AddPosition -3 0 0
> ren1 AddActor clipActor
> #clipActor AddPosition 0 3 0
> ren1 SetBackground 0 0 0
> [ren1 GetActiveCamera] Azimuth 30
> [ren1 GetActiveCamera] Elevation 30
> [ren1 GetActiveCamera] Dolly 1.5
> ren1 ResetCameraClippingRange
> 
> renWin SetSize 400 400
> iren Initialize
> 
> # render the image
> #
> iren SetUserMethod {wm deiconify .vtkInteract}
> 
> 
> scale .c
> 
> proc toto { val } { clipper SetValue $val; renWin Render }
> pack .c
> .c configure -from 0 -to 1 -resolution 0.1 -command toto
> .c set 0.5
>  
> 
> --
> ********************************************************************
> ODET Christophe
> CREATIS, INSA 502, F-69621 Villeurbanne cedex (France)
> tel 04 72 43 85 62
> Fax 04 72 43 85 26
> Email: christophe.odet at creatis.insa-lyon.fr
> URL    :http://www.creatis.insa-lyon.fr
> ********************************************************************

--------------------------------------------------------------------
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.
--------------------------------------------------------------------