Rendering/Volume/vtkEncodedGradientEstimator.cxx

/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkEncodedGradientEstimator.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
#include "vtkEncodedGradientEstimator.h"

#include "vtkGarbageCollector.h"
#include "vtkImageData.h"
#include "vtkMultiThreader.h"
#include "vtkRecursiveSphereDirectionEncoder.h"
#include "vtkTimerLog.h"

#include <cmath>


vtkCxxSetObjectMacro(vtkEncodedGradientEstimator, InputData, vtkImageData );

// Construct a vtkEncodedGradientEstimator with initial values of nullptr for
// the Input, EncodedNormal, and GradientMagnitude. Also,
// indicate that the IndexTable has not yet been initialized. The
// GradientMagnitudeRange and the GradientMangitudeTable are
// initialized to default values - these will change in the future
// when magnitude of gradient opacities are included
vtkEncodedGradientEstimator::vtkEncodedGradientEstimator()
{
  this->InputData                      = nullptr;
  this->EncodedNormals             = nullptr;
  this->EncodedNormalsSize[0]      = 0;
  this->EncodedNormalsSize[1]      = 0;
  this->EncodedNormalsSize[2]      = 0;
  this->GradientMagnitudes         = nullptr;
  this->GradientMagnitudeScale     = 1.0;
  this->GradientMagnitudeBias      = 0.0;
  this->Threader                   = vtkMultiThreader::New();
  this->NumberOfThreads            = this->Threader->GetNumberOfThreads();
  this->DirectionEncoder           = vtkRecursiveSphereDirectionEncoder::New();
  this->ComputeGradientMagnitudes  = 1;
  this->CylinderClip               = 0;
  this->CircleLimits               = nullptr;
  this->CircleLimitsSize           = -1;
  this->UseCylinderClip            = 0;
  this->LastUpdateTimeInSeconds    = -1.0;
  this->LastUpdateTimeInCPUSeconds = -1.0;
  this->ZeroNormalThreshold        = 0.0;
  this->ZeroPad                    = 1;
  this->BoundsClip                 = 0;
  this->Bounds[0] =
    this->Bounds[1] =
    this->Bounds[2] =
    this->Bounds[3] =
    this->Bounds[4] =
    this->Bounds[5] = 0;

}

// Destruct a vtkEncodedGradientEstimator - free up any memory used
vtkEncodedGradientEstimator::~vtkEncodedGradientEstimator()
{
  this->SetInputData(nullptr);
  this->Threader->Delete();
  this->Threader = nullptr;

  delete [] this->EncodedNormals;

  delete [] this->GradientMagnitudes;

  if ( this->DirectionEncoder )
  {
    this->DirectionEncoder->UnRegister( this );
  }

  delete [] this->CircleLimits;
}

void vtkEncodedGradientEstimator::SetZeroNormalThreshold( float v )
{
  if ( this->ZeroNormalThreshold != v )
  {
    if ( v < 0.0 )
    {
      vtkErrorMacro( << "The ZeroNormalThreshold must be a value >= 0.0" );
      return;
    }

    this->ZeroNormalThreshold = v;
    this->Modified();
  }
}

void
vtkEncodedGradientEstimator::SetDirectionEncoder(vtkDirectionEncoder *direnc)
{
  // If we are setting it to its current value, don't do anything
  if ( this->DirectionEncoder == direnc )
  {
    return;
  }

  // If we already have a direction encoder, unregister it.
  if ( this->DirectionEncoder )
  {
    this->DirectionEncoder->UnRegister(this);
    this->DirectionEncoder = nullptr;
  }

  // If we are passing in a non-nullptr encoder, register it
  if ( direnc )
  {
    direnc->Register( this );
  }

  // Actually set the encoder, and consider the object Modified
  this->DirectionEncoder = direnc;
  this->Modified();
}

int vtkEncodedGradientEstimator::GetEncodedNormalIndex( vtkIdType xyzIndex )
{
  this->Update();
  return *(this->EncodedNormals + xyzIndex);
}

int vtkEncodedGradientEstimator::GetEncodedNormalIndex( int xIndex,
                                                        int yIndex,
                                                        int zIndex )
{
  vtkIdType ystep, zstep;

  this->Update();

  // Compute steps through the volume in x, y, and z
  ystep = this->InputSize[0];
  zstep = ystep * this->InputSize[1];

  return *(this->EncodedNormals + zIndex * zstep + yIndex * ystep + xIndex);
}

unsigned short *vtkEncodedGradientEstimator::GetEncodedNormals()
{
  this->Update();

  return this->EncodedNormals;
}

unsigned char *vtkEncodedGradientEstimator::GetGradientMagnitudes()
{
  this->Update();

  return this->GradientMagnitudes;
}

void vtkEncodedGradientEstimator::Update( )
{
  int                scalarInputSize[3];
  double             scalarInputAspect[3];
  double             startSeconds, endSeconds;
  double             startCPUSeconds, endCPUSeconds;

  if ( !this->InputData )
  {
    vtkErrorMacro(<< "No input in gradient estimator.");
    return;
  }

  if ( this->GetMTime() > this->BuildTime ||
       this->DirectionEncoder->GetMTime() > this->BuildTime ||
       this->InputData->GetMTime() > this->BuildTime ||
       !this->EncodedNormals )
  {

    startSeconds = vtkTimerLog::GetUniversalTime();
    startCPUSeconds = vtkTimerLog::GetCPUTime();

    // Get the dimensions of the data and its aspect ratio
    this->InputData->GetDimensions( scalarInputSize );
    this->InputData->GetSpacing( scalarInputAspect );

    // If we previously have allocated space for the encoded normals,
    // and this space is no longer the right size, delete it
    if ( this->EncodedNormalsSize[0] != scalarInputSize[0] ||
         this->EncodedNormalsSize[1] != scalarInputSize[1] ||
         this->EncodedNormalsSize[2] != scalarInputSize[2] )
    {
      delete [] this->EncodedNormals;
      this->EncodedNormals = nullptr;

      delete [] this->GradientMagnitudes;
      this->GradientMagnitudes = nullptr;
    }

    // Compute the number of encoded voxels
    vtkIdType encodedSize = scalarInputSize[0];
    encodedSize *= scalarInputSize[1];
    encodedSize *= scalarInputSize[2];

    // Allocate space for the encoded normals if necessary
    if ( !this->EncodedNormals )
    {
      this->EncodedNormals = new unsigned short[ encodedSize ];
      this->EncodedNormalsSize[0] = scalarInputSize[0];
      this->EncodedNormalsSize[1] = scalarInputSize[1];
      this->EncodedNormalsSize[2] = scalarInputSize[2];
    }

    if ( !this->GradientMagnitudes && this->ComputeGradientMagnitudes )
    {
      this->GradientMagnitudes = new unsigned char[ encodedSize ];
    }

    // Copy info that multi threaded function will need into temp variables
    memcpy( this->InputSize, scalarInputSize, 3 * sizeof(int) );
    // TODO cleanup when double changes are further along
    this->InputAspect[0] = static_cast<float>(scalarInputAspect[0]);
    this->InputAspect[1] = static_cast<float>(scalarInputAspect[1]);
    this->InputAspect[2] = static_cast<float>(scalarInputAspect[2]);
    // memcpy( this->InputAspect, scalarInputAspect, 3 * sizeof(float) );

    if ( this->CylinderClip &&
         (this->InputSize[0] == this->InputSize[1]) )
    {
      this->UseCylinderClip = 1;
      this->ComputeCircleLimits( this->InputSize[0] );
    }
    else
    {
      this->UseCylinderClip = 0;
    }
    this->UpdateNormals();

    this->BuildTime.Modified();

    endSeconds = vtkTimerLog::GetUniversalTime();
    endCPUSeconds = vtkTimerLog::GetCPUTime();

    this->LastUpdateTimeInSeconds    = static_cast<float>(endSeconds    - startSeconds);
    this->LastUpdateTimeInCPUSeconds = static_cast<float>(endCPUSeconds - startCPUSeconds);
  }
}

void vtkEncodedGradientEstimator::ComputeCircleLimits( int size )
{
  int     *ptr, y;
  double  w, halfsize, length, start, end;

  if ( this->CircleLimitsSize != size )
  {
    delete [] this->CircleLimits;
    this->CircleLimits = new int[2*size];
    this->CircleLimitsSize = size;
  }

  ptr = this->CircleLimits;

  halfsize = (size-1)/2.0;

  for ( y = 0; y < size; y++ )
  {
    w = halfsize - y;
    length = static_cast<int>( sqrt( (halfsize*halfsize) - (w*w) ) + 0.5 );
    start = halfsize - length - 1;
    end   = halfsize + length + 1;
    start = (start<0)?(0):(start);
    end   = (end>(size-1))?(size-1):(end);

    *(ptr++) = static_cast<int>(start);
    *(ptr++) = static_cast<int>(end);
  }
}

// Print the vtkEncodedGradientEstimator
void vtkEncodedGradientEstimator::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  if ( this->InputData )
  {
    os << indent << "InputData: (" << this->InputData << ")\n";
  }
  else
  {
    os << indent << "Input: (none)\n";
  }

  if ( this->DirectionEncoder )
  {
    os << indent << "DirectionEncoder: (" << this->DirectionEncoder << ")\n";
  }
  else
  {
    os << indent << "DirectionEncoder: (none)\n";
  }

  os << indent << "Build Time: "
     << this->BuildTime.GetMTime() << endl;

  os << indent << "Gradient Magnitude Scale: "
     << this->GradientMagnitudeScale << endl;

  os << indent << "Gradient Magnitude Bias: "
     << this->GradientMagnitudeBias << endl;

  os << indent << "Zero Pad: "
     << ((this->ZeroPad)?"On":"Off") << endl;

  os << indent << "Bounds Clip: "
     << ((this->BoundsClip)?"On":"Off") << endl;

  os << indent << "Bounds: ("
     << this->Bounds[0] << ", " << this->Bounds[1] << ", "
     << this->Bounds[2] << ", " << this->Bounds[3] << ", "
     << this->Bounds[4] << ", " << this->Bounds[5] << ")\n";

  os << indent << "Zero Normal Threshold: "
     << this->ZeroNormalThreshold << endl;

  os << indent << "Compute Gradient Magnitudes: "
     << ((this->ComputeGradientMagnitudes)?"On":"Off") << endl;

  os << indent << "Cylinder Clip: "
     << ((this->CylinderClip)?"On":"Off") << endl;

  os << indent << "Number Of Threads: "
     << this->NumberOfThreads << endl;

  os << indent << "Last Update Time In Seconds: "
     << this->LastUpdateTimeInSeconds << endl;

  os << indent << "Last Update Time In CPU Seconds: "
     << this->LastUpdateTimeInCPUSeconds << endl;

  // I don't want to print out these variables - they are
  // internal and the get methods are included only for access
  // within the threaded function
  // os << indent << "Use Cylinder Clip: "
  //    << this->UseCylinderClip << endl;
  // os << indent << " Input Size: "
  //    << this->InputSize << endl;
  // os << indent << " Input Aspect Clip: "
  //    << this->InputAspect << endl;

}

//----------------------------------------------------------------------------
void
vtkEncodedGradientEstimator::ReportReferences(vtkGarbageCollector* collector)
{
  this->Superclass::ReportReferences(collector);
  vtkGarbageCollectorReport(collector, this->InputData, "Input");
}