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

Copyright (c) 1998-2010 Kitware Inc. 28 Corporate Drive, Suite 204,
Clifton Park, NY, 12065, USA.

All rights reserved. No part of this software may be reproduced,
distributed,
or modified, in any form or by any means, without permission in writing from
Kitware Inc.

IN NO EVENT SHALL THE AUTHORS OR DISTRIBUTORS BE LIABLE TO ANY PARTY FOR
DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT
OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY DERIVATIVES THEREOF,
EVEN IF THE AUTHORS HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES,
INCLUDING,
BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE, AND NON-INFRINGEMENT.  THIS SOFTWARE IS PROVIDED ON AN
"AS IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE NO OBLIGATION TO
PROVIDE
MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
=========================================================================*/
#include "vtkCMBSphericalPointSource.h"

#include "vtkDoubleArray.h"
#include "vtkIdTypeArray.h"
#include "vtkCellArray.h"
#include "vtkCellLocator.h"
#include "vtkMath.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkPolyData.h"
#include "vtkTransform.h"
#include "vtkPointData.h"
#include "vtkSmartPointer.h"
#include <math.h>

vtkStandardNewMacro(vtkCMBSphericalPointSource);

//----------------------------------------------------------------------------
// Construct with defaults
vtkCMBSphericalPointSource::vtkCMBSphericalPointSource()
{
  this->Radius = 0.0;
  this->RResolution = 1;
  this->ThetaResolution = 1;
  this->PhiResolution = 1;
  this->IgnorePhi = false;
  this->Center[0] = this->Center[1] = this->Center[2] = 0.0;
  this->SetNumberOfInputPorts(0);
}

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

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

  // Determine the number of points we will be evaluating as well as the deltas
  int size;
  double deltaR, deltaTheta, deltaPhi;
  const double PI = 3.141592654;
  deltaTheta = 2.0 * PI / static_cast<double>(this->ThetaResolution);
  deltaPhi = PI / static_cast<double>(this->PhiResolution - 1);
  deltaR = this->Radius;

  if (this->IgnorePhi) // Evaluating a disk
    {
    if ((this->Radius == 0.0) ||
        ((this->RResolution == 1) && (this->ThetaResolution == 1)))
      {
      size = 1;
      }
    else
      {
      if (this->RResolution > 1)
        {
        size = ((this->RResolution - 1) * this->ThetaResolution) + 1;
        deltaR = this->Radius / static_cast<double>(this->RResolution-1);
        }
      else
        {
        size = this->ThetaResolution;
        }
      }
    }
  else // Evaluating a Sphere
     {
     if (this->RResolution > 1)
       {
       // We need to add points along the the z axis which is why there is a + 2 as well as the
       // single point at the origin (the reason for the +1)
       size = (this->RResolution - 1) * (this->ThetaResolution * (this->PhiResolution-2) + 2) + 1;
       deltaR = this->Radius / static_cast<double>(this->RResolution-1);
       }
    else
      {
      size = (this->ThetaResolution * (this->PhiResolution-2)) + 2;
      }
     }

  // Lets allocate the arrays
  vtkSmartPointer<vtkPoints> newPoints;
  newPoints = vtkSmartPointer<vtkPoints>::New();
  newPoints->SetDataTypeToDouble(); //used later during transformation
  newPoints->Allocate(size);
  output->SetPoints(newPoints);
  vtkIdType nextPt;
  vtkSmartPointer<vtkCellArray> outVerts = vtkSmartPointer<vtkCellArray>::New();
  output->SetVerts(outVerts);

  // Populate the data
  if (this->Radius == 0.0)
    {
    // Just test the center Point
    nextPt = newPoints->InsertNextPoint(this->Center);
    outVerts->InsertNextCell(1, &nextPt);
    return 1;
    }

  double theta, r, phi, p[3];
  int iTheta, iPhi, iR=0;

  // Insert Center if needed
  if (this->RResolution > 1)
    {
    nextPt = newPoints->InsertNextPoint(this->Center);
    outVerts->InsertNextCell(1, &nextPt);
    iR = 1;
    }
  if (this->IgnorePhi)
    {
    p[2] = this->Center[2];
    for (r = deltaR;
         iR < this->RResolution;
         iR++, r+=deltaR)
      {
      for (theta = 0.0, iTheta = 0;
           iTheta < this->ThetaResolution;
           iTheta++, theta+=deltaTheta)
        {
        p[0] = this->Center[0] + (r*cos(theta));
        p[1] = this->Center[1] + (r*sin(theta));
        nextPt = newPoints->InsertNextPoint(p);
        outVerts->InsertNextCell(1, &nextPt);
        }
      }
    }
  else
    {
    // First insert the points along the Phi = 0 axis
    p[0] = this->Center[0];
    p[1] = this->Center[1];
    //Are we inserting points in the r direction?
    if (iR == 1)
      {
      for (r = deltaR;
           iR < this->RResolution;
           iR++, r+=deltaR)
        {
        p[2] = this->Center[2] + r;
        nextPt = newPoints->InsertNextPoint(p);
        outVerts->InsertNextCell(1, &nextPt);
        p[2] = this->Center[2] - r;
        nextPt = newPoints->InsertNextPoint(p);
        outVerts->InsertNextCell(1, &nextPt);
        }
      iR = 1;
      }
    else
      {
      // Just insert points at the poles
        p[2] = this->Center[2] + this->Radius;
        nextPt = newPoints->InsertNextPoint(p);
        outVerts->InsertNextCell(1, &nextPt);
        p[2] = this->Center[2] - this->Radius;
        nextPt = newPoints->InsertNextPoint(p);
        outVerts->InsertNextCell(1, &nextPt);
      }
    double sphi;
    for (r = deltaR;
         iR < this->RResolution;
         iR++, r+=deltaR)
      {
      for (phi = deltaPhi, iPhi = 2;
           iPhi < this->PhiResolution;
           iPhi++, phi+=deltaPhi)
        {
        p[2] = this->Center[2] +  (r*cos(phi));
        sphi =  (r*sin(phi));

        for (theta = 0.0, iTheta = 0;
             iTheta < this->ThetaResolution;
             iTheta++, theta+=deltaTheta)
          {
          p[0] = this->Center[0] + (sphi*cos(theta));
          p[1] = this->Center[1] + (sphi*sin(theta));
          nextPt = newPoints->InsertNextPoint(p);
          outVerts->InsertNextCell(1, &nextPt);
          }
        }
      }
    }
  return 1;
}

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

  os << indent << "RResolution: " << this->RResolution << "\n";
  os << indent << "ThetaResolution: " << this->ThetaResolution << "\n";
  os << indent << "PhiResolution: " << this->PhiResolution << "\n";
  os << indent << "IgnorePhi: " << (this->IgnorePhi ? "On\n" : "Off\n");
  os << indent << "Center: (" << this->Center[0] << ", "
     << this->Center[1] << ", " << this->Center[2] << ")\n";
}
