[vtkusers] My investigation results on the three vtkOpenGLGPUVolumeRayCastMapper bugs. Please read!

[AGPX]

Hi to all,

I have the time to address the 3 ancients bugs that affect the mapper 
vtkOpenGLGPUVolumeRayCastMapper, discovered and signaled on mantis bugtracer by 
me on February 2009 (bugs never solved). These bugs make the mapper quite 
unusable and this is a serious issue because that mapper is the best compromise 
between quality and speed for volume rendering.
OK, we have 3 different issue. Let me talk about the three.

1) Plane clipping bugs (to better figure out the problem, here a video of the 
problem: http://www.youtube.com/watch?v=NjUTTA86yTM&feature=related). I have 
solved this bug. The problem is in the method: 
vtkOpenGLGPUVolumeRayCastMapper::RenderClippedBoundingBox. That method check the 
correct winding order of the polygon and fails its job. Actually the algorithm 
check the relationship between the polygon plane and the SPATIAL center of the 
volume. Note that the volume is always convex. The problem is that you haven't 
to use spatial center, but center of mass. If you use spatial center (center of 
bounding box) you wrong because it not always lie inside the convex polyhedra: 
it could lie on a polygon boundary plane. Let me show this in 2D:

 __________
|         /
|        /
|       /
|      /
|     x
|    /
|   /
|  /
| /
|/

If you have a triangle, the spatial center lie on an edge and due to floating 
point inaccuracy this point could lie in front or behind the plane. When point 
lie behind the plane the algorithm flip the vertex order resulting in back-faced 
polygon. The solution is to compute center of mass (average of vertices):

 __________
|         /
|        /
|       /
|      /
|  x  /
|    /
|   /
|  /
| /
|/

The center of mass always lie inside the convex polyhedra and this solve the 
problem.

Here my fix (in vtkOpenGLGPUVolumeRayCastMapper::RenderClippedBoundingBox):

  double center[3] = {0,0,0};
  double min[3] = {VTK_DOUBLE_MAX, VTK_DOUBLE_MAX, VTK_DOUBLE_MAX};
  double max[3] = {VTK_DOUBLE_MIN, VTK_DOUBLE_MIN, VTK_DOUBLE_MIN};

  // First compute center point
  npts = points->GetNumberOfPoints();
  for ( i = 0; i < npts; i++ )
    {
    double pt[3];
    points->GetPoint( i, pt );
    for ( j = 0; j < 3; j++ )
      {
          // AGPX MODIFIED
//      min[j] = (pt[j]<min[j])?(pt[j]):(min[j]);
//      max[j] = (pt[j]>max[j])?(pt[j]):(max[j]);
          center[j] += pt[j];
          // AGPX MODIFIED
      }
    }

 // AGPX MODIFIED
  //center[0] = 0.5*(min[0]+max[0]);
  //center[1] = 0.5*(min[1]+max[1]);
  //center[2] = 0.5*(min[2]+max[2]);
  center[0] /= ((double)npts);
  center[1] /= ((double)npts);
  center[2] /= ((double)npts);
 // AGPX MODIFIED

This completely solve this issue. Now go to the next bug.

2) "Volume disappearing". Suddenly, if you fly-through inside a volume, it can 
disappear. The problem is related to the near clipping plane. Actually the 
method vtkOpenGLGPUVolumeRayCastMapper::ClipBoundingBox, clip the volume box 
with the near clipping plane (and eventually also with custom clip planes) in 
order to avoid that the volume results opened. The problem is that due to 
floating point inaccuracy the resulting clipped convex polyhedra could be still 
opened! The problem could be solved pushing the near clipping plane slightly 
ahead. That is:

if(this->NearPlane==0)
    {
    this->NearPlane= vtkPlane::New();
    }

  // AGPX ADDED
  const double offset = 0.01;
  camNearPoint[0] += camPlaneNormal[0] * offset;
  camNearPoint[1] += camPlaneNormal[1] * offset;
  camNearPoint[2] += camPlaneNormal[2] * offset;
  // AGPX ADDED

  this->NearPlane->SetOrigin( camNearPoint );
  this->NearPlane->SetNormal( camPlaneNormal );
  this->Planes->AddItem(this->NearPlane);

The problem here is to establish how much you have to push plane ahead. 0.01 is 
working for me. Ok, following me to the last (hard) bug.

3) Volume Rendering Deformation. This bug show a non correct rendering of the 
volume when you perform a fly-through inside a volume (here 
http://www.youtube.com/watch?gl=IT&hl=it&v=C_GA6UHbUGo you can see a video of 
the problem). Actually this bug is unsolved and here I need some help from 
authors. Look like a bad entry-ray and/or ray-direction calculation. The problem 
is showed only when the near clipping plane clips the volume box. The entry 
point is given by the texture coordinate of the point being rasterized (a 3D 
texture coordinate is applied to the vertices of the rendered volume box). This 
look like correct. The ray direction is given by the difference between the 
entry point and the camera position in texture space. This should be always ok. 
The method vtkOpenGLGPUVolumeRayCastMapper::LoadProjectionParameters computed 
the camera position in texture space and looks correct (notice that I'm always 
talking about a perspective camera):

double *bounds=this->CurrentScalar->GetLoadedBounds();

  double dx=bounds[1]-bounds[0];
  double dy=bounds[3]-bounds[2];
  double dz=bounds[5]-bounds[4];


 // Perspective projection

    // Compute camera position in texture coordinates
    // Position of the center of the camera in world frame
    double cameraPosWorld[4];
    // Position of the center of the camera in the dataset frame
    // (the transform of the volume is taken into account)
    double cameraPosDataset[4];
   // Position of the center of the camera in the texture frame
    // the coordinates are translated and rescaled
    double cameraPosTexture[4];

    ren->GetActiveCamera()->GetPosition(cameraPosWorld);
    cameraPosWorld[3]=1.0; // we use homogeneous coordinates.

    datasetToWorld->MultiplyPoint(cameraPosWorld,cameraPosDataset);

   // From homogeneous to cartesian coordinates.
    if(cameraPosDataset[3]!=1.0)
      {
      double ratio=1/cameraPosDataset[3];
      cameraPosDataset[0]*=ratio;
      cameraPosDataset[1]*=ratio;
      cameraPosDataset[2]*=ratio;
      }

    cameraPosTexture[0] = (cameraPosDataset[0]-bounds[0])/dx;
    cameraPosTexture[1] = (cameraPosDataset[1]-bounds[2])/dy;
    cameraPosTexture[2] = (cameraPosDataset[2]-bounds[4])/dz;

Ray direction is computed by the shader 
(vtkGPUVolumeRayCastMapper_PerspectiveProjectionFS.glsl):

// Entry position (global scope)
vec3 pos;
// Incremental vector in texture space (global scope)
vec3 rayDir;
// Camera position in texture space
uniform vec3 cameraPosition;

// Defined in the right projection method.
void incrementalRayDirection()
{
  // Direction of the ray in texture space, not normalized.
  rayDir=pos-cameraPosition;
  ...
}

and looks correct.

The entry point seems correct, because if I use a totally opaque volume (in 
order to avoid ray traversal, that is the algorithm stop immediately on the 
first hitted voxel) and clipped with a plane to enter inside it, the rendering 
is always ok. That is entry point looks correct. When the volume isn't fully 
opaque the result is wrong. So, I suspect that the ray direction is wrong. Note 
that the problem disappear (or at least is not much noticeable) IF I try push 
the near clipping plane ahead of a greater quantity (>> 0.01). If I push it 
ahead of 16.0 instead of 0.01 (as stated in the previous bug), the deformation 
is not very noticeable. Note that 16.0 is an empirically determinated value for 
my tested datasets and reasonably it's not a valid choice for every dataset (a 
way to determine a good value?). It's not clear why the near clipping plane 
cause this problem. The texture coordinate seems to be interpolated correctly 
and apparently there are no reason why the method should not works. Here I'm 
stuck and I need some help from some CG Guru! Please support!

Thanks in advance for your attention,

- Gianluca Arcidiacono (a.k.a. AGPX)



      
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20100812/eb5e4d68/attachment.htm>