Wednesday, May 1, 2013

Generate a Torus in VTK without using the superquadric function

Older versions of VTK only provide the superquadric function which is also able to generate a Torus. This post provides a vtkTorusSource which directly generates a vtkPolyData which is Torus shaped.
Phi is the angle of the Torus' body, theta the angle of the circle in the centre of the Torus. This means that the whole extent is [-Radius - BodyRadius, Radius + BodyRadius] and the diameter of the inner whole is 2 * (Radius - BodyRadius). The other parameters should be clear by their names.
vtkTorusSource.h
#ifndef VTKTORUSSOURCE_H
#define VTKTORUSSOURCE_H

#include "vtkPolyDataAlgorithm.h"

#define VTK_MAX_TORUS_RESOLUTION 1024

class VTK_GRAPHICS_EXPORT vtkTorusSource : public vtkPolyDataAlgorithm
{
public:
  vtkTypeMacro(vtkTorusSource, vtkPolyDataAlgorithm);
  void PrintSelf(ostream& os, vtkIndent indent);

  // Description:
  // Construct torus with default resolution 8 in both Phi
  // and Theta directions. Theta ranges from (0,360) and phi (0,360) degrees.
  static vtkTorusSource *New();

  // Description: Radius of the center circle of the torus
  // Set radius of sphere. Default is 2. => together with BodyRadius this gives an
  // outer radius of 3
  vtkSetClampMacro(Radius,double,0.0,VTK_DOUBLE_MAX);
  vtkGetMacro(Radius,double);

  // Description: Radius of the body of the torus
  // Set radius of sphere. Default is 1.
  vtkSetClampMacro(BodyRadius,double,0.0,VTK_DOUBLE_MAX);
  vtkGetMacro(BodyRadius,double);

  // Description:
  // Set the center of the sphere. Default is 0,0,0.
  vtkSetVector3Macro(Center,double);
  vtkGetVectorMacro(Center,double,3);

  // Description:
  // Set the number of points in the longitude direction (ranging from
  // StartTheta to EndTheta).
  vtkSetClampMacro(ThetaResolution,int,3,VTK_MAX_TORUS_RESOLUTION);
  vtkGetMacro(ThetaResolution,int);

  // Description:
  // Set the number of points in the latitude direction (ranging
  // from StartPhi to EndPhi).
  vtkSetClampMacro(PhiResolution,int,3,VTK_MAX_TORUS_RESOLUTION);
  vtkGetMacro(PhiResolution,int);

  // Description:
  // Set the starting longitude angle. By default StartTheta=0 degrees.
  vtkSetClampMacro(StartTheta,double,0.0,360.0);
  vtkGetMacro(StartTheta,double);

  // Description:
  // Set the ending longitude angle. By default EndTheta=360 degrees.
  vtkSetClampMacro(EndTheta,double,0.0,360.0);
  vtkGetMacro(EndTheta,double);

  // Description:
  // Set the starting latitude angle (0 is at north pole). By default
  // StartPhi=0 degrees.
  vtkSetClampMacro(StartPhi,double,0.0,360.0);
  vtkGetMacro(StartPhi,double);

  // Description:
  // Set the ending latitude angle. By default EndPhi=360 degrees.
  vtkSetClampMacro(EndPhi,double,0.0,360.0);
  vtkGetMacro(EndPhi,double);

protected:
  vtkTorusSource();

  int RequestData(vtkInformation *, vtkInformationVector **,
    vtkInformationVector *);
  int RequestInformation(vtkInformation *, vtkInformationVector **,
   vtkInformationVector *);

  double Radius;
  double BodyRadius;
  double Center[3];
  int ThetaResolution;
  int PhiResolution;
  double StartTheta;
  double EndTheta;
  double StartPhi;
  double EndPhi;

private:
  vtkTorusSource(const vtkTorusSource&);  // Not implemented.
  void operator=(const vtkTorusSource&);  // Not implemented.
};

#endif // VTKTORUSSOURCE_H
vtkTorusSource.cxx
#include "vtkTorusSource.h"

#include <vtkObjectFactory.h>
#include <vtkPointData.h>
#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkDataObject.h>
#include <vtkStreamingDemandDrivenPipeline.h>
#include <vtkInformation.h>
#include <vtkInformationVector.h>
#include <vtkFloatArray.h>
#include <vtkCellArray.h>
#include <vtkMath.h>

vtkStandardNewMacro(vtkTorusSource);

vtkTorusSource::vtkTorusSource() : vtkPolyDataAlgorithm()
{
  Radius = 2.;
  BodyRadius = 1.;
  Center[0] = 0.0;
  Center[1] = 0.0;
  Center[2] = 0.0;

  ThetaResolution = 8;
  PhiResolution = 8;
  StartTheta = 0.0;
  EndTheta = 360.0;
  StartPhi = 0.0;
  EndPhi = 360.0;

  this->SetNumberOfInputPorts(0);
}

int vtkTorusSource::RequestData(
    vtkInformation *vtkNotUsed(request),
    vtkInformationVector **vtkNotUsed(inputVector),
    vtkInformationVector *outputVector)
{
  // get the info object
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

  // get the ouptut
  vtkPolyData *output = vtkPolyData::SafeDownCast(
        outInfo->Get(vtkDataObject::DATA_OBJECT()));

  int i, j;
  int jStart, jEnd;
  int numPts, numPolys;
  vtkPoints *newPoints;
  vtkFloatArray *newNormals;
  vtkCellArray *newPolys;
  double x[3], n[3], deltaPhi, deltaTheta, phi, theta, radius, norm;
  double startTheta, endTheta, startPhi, endPhi;
  int base, thetaResolution, phiResolution;
  vtkIdType pts[4];
  int piece =
      outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
  int numPieces =
      outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());

  if (numPieces > ThetaResolution)
  {
    numPieces = ThetaResolution;
  }
  if (piece >= numPieces)
  {
    // Although the super class should take care of this,
    // it cannot hurt to check here.
    return 1;
  }

  // I want to modify the ivars resoultion start theta and end theta,
  // so I will make local copies of them.  THese might be able to be merged
  // with the other copies of them, ...
  int localThetaResolution = ThetaResolution;
  double localStartTheta = StartTheta;
  double localEndTheta = EndTheta;

  while (localEndTheta < localStartTheta)
  {
    localEndTheta += 360.0;
  }
  deltaTheta = (localEndTheta - localStartTheta) / localThetaResolution;

  // Change the ivars based on pieces.
  int start, end;
  start = piece * localThetaResolution / numPieces;
  end = (piece+1) * localThetaResolution / numPieces;
  localEndTheta = localStartTheta + (double)(end) * deltaTheta;
  localStartTheta = localStartTheta + (double)(start) * deltaTheta;
  localThetaResolution = end - start;

  // Set things up; allocate memory
  //
  vtkDebugMacro("TorusSource Executing piece index " << piece
                << " of " << numPieces << " pieces.");

  numPts = this->PhiResolution * localThetaResolution + 2;
  numPolys = this->PhiResolution * 2 * localThetaResolution;

  newPoints = vtkPoints::New();
  newPoints->Allocate(numPts);
  newNormals = vtkFloatArray::New();
  newNormals->SetNumberOfComponents(3);
  newNormals->Allocate(3*numPts);
  newNormals->SetName("Normals");

  newPolys = vtkCellArray::New();
  newPolys->Allocate(newPolys->EstimateSize(numPolys, 4));

  // Check data, determine increments, and convert to radians
  startTheta = (localStartTheta < localEndTheta ? localStartTheta : localEndTheta);
  startTheta *= vtkMath::Pi() / 180.0;
  endTheta = (localEndTheta > localStartTheta ? localEndTheta : localStartTheta);
  endTheta *= vtkMath::Pi() / 180.0;

  startPhi = (this->StartPhi < this->EndPhi ? this->StartPhi : this->EndPhi);
  startPhi *= vtkMath::Pi() / 180.0;
  endPhi = (this->EndPhi > this->StartPhi ? this->EndPhi : this->StartPhi);
  endPhi *= vtkMath::Pi() / 180.0;

  phiResolution = this->PhiResolution;
  deltaPhi = (endPhi - startPhi) / (this->PhiResolution - 1);
  thetaResolution = localThetaResolution;
  if (fabs(localStartTheta - localEndTheta) < 360.0)
  {
    ++localThetaResolution;
  }
  deltaTheta = (endTheta - startTheta) / thetaResolution;

  jStart = 0;
  jEnd = PhiResolution;

  this->UpdateProgress(0.1);

  // Create intermediate points
  for (i = 0; i < localThetaResolution; ++i)
  {
    theta = localStartTheta * vtkMath::Pi() / 180.0 + i * deltaTheta;

    for (j = jStart; j < jEnd; ++j)
    {
      phi = startPhi + j * deltaPhi;
      n[0] = BodyRadius * cos(phi) * cos(theta);
      n[1] = BodyRadius * sin(phi) * sin(theta);
      n[2] = BodyRadius * sin(phi);
      x[0] = (Radius + BodyRadius * cos(phi)) * cos(theta);
      x[1] = (Radius + BodyRadius * cos(phi)) * sin(theta);
      x[2] = BodyRadius * sin(phi);
      newPoints->InsertNextPoint(x);

      if ( (norm = vtkMath::Norm(n)) == 0.0 )
      {
        norm = 1.0;
      }
      n[0] /= norm; n[1] /= norm; n[2] /= norm;
      newNormals->InsertNextTuple(n);
    }
    this->UpdateProgress (0.10 + 0.50 * i / static_cast<float>(localThetaResolution));
  }

  // Generate mesh connectivity
  base = phiResolution * localThetaResolution;

  if (fabs(localStartTheta - localEndTheta) < 360.0)
  {
    --localThetaResolution;
  }

  this->UpdateProgress (0.70);

  for (i = 0; i < localThetaResolution; ++i)
  {
    for (j = 0; j < (phiResolution - 1); ++j)
    {
      pts[1] = phiResolution * i + j;
      pts[0] = pts[1] + 1;
      pts[3] = ((phiResolution * (i + 1) + j) % base) + 1;
      pts[2] = pts[3] - 1;
      newPolys->InsertNextCell(4, pts);
    }
    this->UpdateProgress (0.70 + 0.30 * i / static_cast<double>(localThetaResolution));
  }

  // Update ourselves and release memory
  //
  newPoints->Squeeze();
  output->SetPoints(newPoints);
  newPoints->Delete();

  newNormals->Squeeze();
  output->GetPointData()->SetNormals(newNormals);
  newNormals->Delete();

  newPolys->Squeeze();
  output->SetPolys(newPolys);
  newPolys->Delete();

  return 1;
}

void vtkTorusSource::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);

  os << indent << "Theta Resolution: " << ThetaResolution << "\n";
  os << indent << "Phi Resolution: " << PhiResolution << "\n";
  os << indent << "Theta Start: " << StartTheta << "\n";
  os << indent << "Phi Start: " << StartPhi << "\n";
  os << indent << "Theta End: " << EndTheta << "\n";
  os << indent << "Phi End: " << EndPhi << "\n";
  os << indent << "Radius: " << Radius << "\n";
  os << indent << "BodyRadius: " << BodyRadius << "\n";
  os << indent << "Center: (" << Center[0] << ", "
     << Center[1] << ", " << Center[2] << ")\n";
}

int vtkTorusSource::RequestInformation(
    vtkInformation *vtkNotUsed(request),
    vtkInformationVector **vtkNotUsed(inputVector),
    vtkInformationVector *outputVector)
{
  // get the info object
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

  outInfo->Set(vtkStreamingDemandDrivenPipeline::MAXIMUM_NUMBER_OF_PIECES(),
               -1);

  outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_BOUNDING_BOX(),
               Center[0] - Radius - BodyRadius,
      Center[0] + Radius + BodyRadius,
      Center[1] - Radius - BodyRadius,
      Center[1] + Radius + BodyRadius,
      Center[2] - BodyRadius,
      Center[2] + BodyRadius);

  return 1;
}
UPDATE: Fixed calculation of normals.

No comments:

Post a Comment