text: /*=========================================================================
  
    Program:   Visualization Toolkit
    Module:    vtkVoxelModeller.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 "vtkVoxelModeller.h"
  
  #include "vtkBitArray.h"
  #include "vtkCell.h"
  #include "vtkImageData.h"
  #include "vtkInformation.h"
  #include "vtkInformationVector.h"
  #include "vtkObjectFactory.h"
  #include "vtkPointData.h"
  #include "vtkStreamingDemandDrivenPipeline.h"
  
  #include <cmath>
  
  vtkStandardNewMacro(vtkVoxelModeller);
  
  // Construct an instance of vtkVoxelModeller with its sample dimensions
  // set to (50,50,50), and so that the model bounds are
  // automatically computed from its input. The maximum distance is set to
  // examine the whole grid. This could be made much faster, and probably
  // will be in the future.
  vtkVoxelModeller::vtkVoxelModeller()
  {
    this->MaximumDistance = 1.0;
  
    this->ModelBounds[0] = 0.0;
    this->ModelBounds[1] = 0.0;
    this->ModelBounds[2] = 0.0;
    this->ModelBounds[3] = 0.0;
    this->ModelBounds[4] = 0.0;
    this->ModelBounds[5] = 0.0;
  
    this->SampleDimensions[0] = 50;
    this->SampleDimensions[1] = 50;
    this->SampleDimensions[2] = 50;
  
    this->ScalarType = VTK_BIT;
    this->ForegroundValue = 1.0;
    this->BackgroundValue = 0.0;
  }
  
  // Specify the position in space to perform the voxelization.
  void vtkVoxelModeller::SetModelBounds(const double bounds[6])
  {
    vtkVoxelModeller::SetModelBounds(
      bounds[0], bounds[1], bounds[2], bounds[3], bounds[4], bounds[5]);
  }
  
  void vtkVoxelModeller::SetModelBounds(
    double xmin, double xmax, double ymin, double ymax, double zmin, double zmax)
  {
    if (this->ModelBounds[0] != xmin || this->ModelBounds[1] != xmax ||
      this->ModelBounds[2] != ymin || this->ModelBounds[3] != ymax || this->ModelBounds[4] != zmin ||
      this->ModelBounds[5] != zmax)
    {
      this->Modified();
      this->ModelBounds[0] = xmin;
      this->ModelBounds[1] = xmax;
      this->ModelBounds[2] = ymin;
      this->ModelBounds[3] = ymax;
      this->ModelBounds[4] = zmin;
      this->ModelBounds[5] = zmax;
    }
  }
  
  int vtkVoxelModeller::RequestInformation(vtkInformation* vtkNotUsed(request),
    vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector)
  {
    // get the info objects
    vtkInformation* outInfo = outputVector->GetInformationObject(0);
  
    int i;
    double ar[3], origin[3];
  
    outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), 0, this->SampleDimensions[0] - 1,
      0, this->SampleDimensions[1] - 1, 0, this->SampleDimensions[2] - 1);
  
    for (i = 0; i < 3; i++)
    {
      origin[i] = this->ModelBounds[2 * i];
      if (this->SampleDimensions[i] <= 1)
      {
        ar[i] = 1;
      }
      else
      {
        ar[i] =
          (this->ModelBounds[2 * i + 1] - this->ModelBounds[2 * i]) / (this->SampleDimensions[i] - 1);
      }
    }
    outInfo->Set(vtkDataObject::ORIGIN(), origin, 3);
    outInfo->Set(vtkDataObject::SPACING(), ar, 3);
  
    vtkDataObject::SetPointDataActiveScalarInfo(outInfo, this->ScalarType, 1);
    return 1;
  }
  
  int vtkVoxelModeller::RequestData(vtkInformation* vtkNotUsed(request),
    vtkInformationVector** inputVector, vtkInformationVector* outputVector)
  {
    // get the input
    vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
    vtkDataSet* input = vtkDataSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
  
    // get the output
    vtkInformation* outInfo = outputVector->GetInformationObject(0);
    vtkImageData* output = vtkImageData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
  
    // We need to allocate our own scalars since we are overriding
    // the superclasses "Execute()" method.
    output->SetExtent(outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()));
    output->AllocateScalars(outInfo);
  
    vtkIdType cellNum, i;
    int j, k;
    const double* bounds;
    double adjBounds[6];
    vtkCell* cell;
    double maxDistance, pcoords[3];
    vtkIdType numPts, idx, numCells;
    int subId;
    int min[3], max[3];
    double x[3], distance2;
    int jkFactor;
    double* weights = new double[input->GetMaxCellSize()];
    double closestPoint[3];
    double voxelHalfWidth[3], origin[3], spacing[3];
    vtkDataArray* newScalars = output->GetPointData()->GetScalars();
  
    //
    // Initialize self; create output objects
    //
    vtkDebugMacro(<< "Executing Voxel model");
  
    numPts = this->SampleDimensions[0] * this->SampleDimensions[1] * this->SampleDimensions[2];
    for (i = 0; i < numPts; i++)
    {
      newScalars->SetComponent(i, 0, this->BackgroundValue);
    }
  
    maxDistance = this->ComputeModelBounds(origin, spacing);
    outInfo->Set(vtkDataObject::SPACING(), spacing, 3);
    outInfo->Set(vtkDataObject::ORIGIN(), origin, 3);
    //
    // Voxel widths are 1/2 the height, width, length of a voxel
    //
    for (i = 0; i < 3; i++)
    {
      voxelHalfWidth[i] = spacing[i] / 2.0;
    }
    //
    // Traverse all cells; computing distance function on volume points.
    //
    numCells = input->GetNumberOfCells();
    for (cellNum = 0; cellNum < numCells; cellNum++)
    {
      cell = input->GetCell(cellNum);
      bounds = cell->GetBounds();
      for (i = 0; i < 3; i++)
      {
        adjBounds[2 * i] = bounds[2 * i] - maxDistance;
        adjBounds[2 * i + 1] = bounds[2 * i + 1] + maxDistance;
      }
  
      // compute dimensional bounds in data set
      for (i = 0; i < 3; i++)
      {
        min[i] = static_cast<int>(static_cast<double>(adjBounds[2 * i] - origin[i]) / spacing[i]);
        max[i] = static_cast<int>(static_cast<double>(adjBounds[2 * i + 1] - origin[i]) / spacing[i]);
        if (min[i] < 0)
        {
          min[i] = 0;
        }
        if (max[i] >= this->SampleDimensions[i])
        {
          max[i] = this->SampleDimensions[i] - 1;
        }
      }
  
      jkFactor = this->SampleDimensions[0] * this->SampleDimensions[1];
      for (k = min[2]; k <= max[2]; k++)
      {
        x[2] = spacing[2] * k + origin[2];
        for (j = min[1]; j <= max[1]; j++)
        {
          x[1] = spacing[1] * j + origin[1];
          for (i = min[0]; i <= max[0]; i++)
          {
            idx = jkFactor * k + this->SampleDimensions[0] * j + i;
            if (!(newScalars->GetComponent(idx, 0)))
            {
              x[0] = spacing[0] * i + origin[0];
  
              if (cell->EvaluatePosition(x, closestPoint, subId, pcoords, distance2, weights) != -1 &&
                ((fabs(closestPoint[0] - x[0]) <= voxelHalfWidth[0]) &&
                  (fabs(closestPoint[1] - x[1]) <= voxelHalfWidth[1]) &&
                  (fabs(closestPoint[2] - x[2]) <= voxelHalfWidth[2])))
              {
                newScalars->SetComponent(idx, 0, this->ForegroundValue);
              }
            }
          }
        }
      }
    }
    delete[] weights;
  
    return 1;
  }
  
  // Compute the ModelBounds based on the input geometry.
  double vtkVoxelModeller::ComputeModelBounds(double origin[3], double spacing[3])
  {
    const double* bounds;
    double maxDist;
    int i, adjustBounds = 0;
  
    // compute model bounds if not set previously
    if (this->ModelBounds[0] >= this->ModelBounds[1] ||
      this->ModelBounds[2] >= this->ModelBounds[3] || this->ModelBounds[4] >= this->ModelBounds[5])
    {
      adjustBounds = 1;
      vtkDataSet* ds = vtkDataSet::SafeDownCast(this->GetInput());
      // ds better be non null otherwise something is very wrong here
      bounds = ds->GetBounds();
    }
    else
    {
      bounds = this->ModelBounds;
    }
  
    for (maxDist = 0.0, i = 0; i < 3; i++)
    {
      if ((bounds[2 * i + 1] - bounds[2 * i]) > maxDist)
      {
        maxDist = bounds[2 * i + 1] - bounds[2 * i];
      }
    }
    maxDist *= this->MaximumDistance;
  
    // adjust bounds so model fits strictly inside (only if not set previously)
    if (adjustBounds)
    {
      for (i = 0; i < 3; i++)
      {
        this->ModelBounds[2 * i] = bounds[2 * i] - maxDist;
        this->ModelBounds[2 * i + 1] = bounds[2 * i + 1] + maxDist;
      }
    }
  
    // Set volume origin and data spacing
    for (i = 0; i < 3; i++)
    {
      origin[i] = this->ModelBounds[2 * i];
      spacing[i] =
        (this->ModelBounds[2 * i + 1] - this->ModelBounds[2 * i]) / (this->SampleDimensions[i] - 1);
    }
  
    return maxDist;
  }
  
  // Set the i-j-k dimensions on which to sample the distance function.
  void vtkVoxelModeller::SetSampleDimensions(int i, int j, int k)
  {
    int dim[3];
  
    dim[0] = i;
    dim[1] = j;
    dim[2] = k;
  
    this->SetSampleDimensions(dim);
  }
  
  void vtkVoxelModeller::SetSampleDimensions(int dim[3])
  {
    int dataDim, i;
  
    vtkDebugMacro(<< " setting SampleDimensions to (" << dim[0] << "," << dim[1] << "," << dim[2]
                  << ")");
  
    if (dim[0] != this->SampleDimensions[0] || dim[1] != this->SampleDimensions[1] ||
      dim[2] != this->SampleDimensions[2])
    {
      if (dim[0] < 1 || dim[1] < 1 || dim[2] < 1)
      {
        vtkErrorMacro(<< "Bad Sample Dimensions, retaining previous values");
        return;
      }
      for (dataDim = 0, i = 0; i < 3; i++)
      {
        if (dim[i] > 1)
        {
          dataDim++;
        }
      }
      if (dataDim < 3)
      {
        vtkErrorMacro(<< "Sample dimensions must define a volume!");
        return;
      }
  
      for (i = 0; i < 3; i++)
      {
        this->SampleDimensions[i] = dim[i];
      }
      this->Modified();
    }
  }
  
  int vtkVoxelModeller::FillInputPortInformation(int vtkNotUsed(port), vtkInformation* info)
  {
    info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
    return 1;
  }
  
  void vtkVoxelModeller::PrintSelf(ostream& os, vtkIndent indent)
  {
    this->Superclass::PrintSelf(os, indent);
  
    os << indent << "Maximum Distance: " << this->MaximumDistance << "\n";
    os << indent << "Sample Dimensions: (" << this->SampleDimensions[0] << ", "
       << this->SampleDimensions[1] << ", " << this->SampleDimensions[2] << ")\n";
    os << indent << "Model Bounds: \n";
    os << indent << "  Xmin,Xmax: (" << this->ModelBounds[0] << ", " << this->ModelBounds[1] << ")\n";
    os << indent << "  Ymin,Ymax: (" << this->ModelBounds[2] << ", " << this->ModelBounds[3] << ")\n";
    os << indent << "  Zmin,Zmax: (" << this->ModelBounds[4] << ", " << this->ModelBounds[5] << ")\n";
    os << indent << "ScalarType: " << this->ScalarType << endl;
    os << indent << "ForegroundValue: " << this->ForegroundValue << endl;
    os << indent << "BackgroundValue: " << this->BackgroundValue << endl;
  }
  