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

  Program:   Visualization Toolkit
  Module:    vtkEllipsoidalGaussianKernel.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 "vtkEllipsoidalGaussianKernel.h"
#include "vtkAbstractPointLocator.h"
#include "vtkObjectFactory.h"
#include "vtkIdList.h"
#include "vtkDoubleArray.h"
#include "vtkDataSet.h"
#include "vtkPointData.h"
#include "vtkMath.h"

vtkStandardNewMacro(vtkEllipsoidalGaussianKernel);

//----------------------------------------------------------------------------
vtkEllipsoidalGaussianKernel::vtkEllipsoidalGaussianKernel()
{
  this->Radius = 1.0;
  this->Sharpness = 2.0;
  this->Eccentricity = 2.0;

  this->Normals = NULL;
  this->Scalars = NULL;
  this->UseNormals = true;
  this->UseScalars = false;

  this->F2 = this->Sharpness / this->Radius;
  this->E2 = this->Eccentricity * this->Eccentricity;
}


//----------------------------------------------------------------------------
vtkEllipsoidalGaussianKernel::~vtkEllipsoidalGaussianKernel()
{
  this->FreeStructures();
}


//----------------------------------------------------------------------------
void vtkEllipsoidalGaussianKernel::
FreeStructures()
{
  this->Superclass::FreeStructures();

  if ( this->Normals )
    {
    this->Normals->Delete();
    this->Normals = NULL;
    }

  if ( this->Scalars )
    {
    this->Scalars->Delete();
    this->Scalars = NULL;
    }
}

//----------------------------------------------------------------------------
void vtkEllipsoidalGaussianKernel::
Initialize(vtkAbstractPointLocator *loc, vtkDataSet *ds, vtkPointData *pd)
{
  this->Superclass::Initialize(loc, ds, pd);

  this->Scalars = pd->GetScalars();
  if ( this->UseScalars && this->Scalars &&
       this->Scalars->GetNumberOfComponents() == 1 )
    {
    this->Scalars->Register(this);
    }
  else
    {
    this->Scalars = NULL;
    }

  this->Normals = pd->GetNormals();
  if ( this->UseNormals && this->Normals )
    {
    this->Normals->Register(this);
    }
  else
    {
    this->Normals = NULL;
    }

  this->F2 = this->Sharpness / this->Radius;
  this->F2 = this->F2 * this->F2;
  this->E2 = this->Eccentricity * this->Eccentricity;
}

//----------------------------------------------------------------------------
vtkIdType vtkEllipsoidalGaussianKernel::
ComputeBasis(double x[3], vtkIdList *pIds)
{
  this->Locator->FindPointsWithinRadius(this->Radius, x, pIds);
  return pIds->GetNumberOfIds();
}

//----------------------------------------------------------------------------
vtkIdType vtkEllipsoidalGaussianKernel::
ComputeWeights(double x[3], vtkIdList *pIds, vtkDoubleArray *weights)
{
  vtkIdType numPts = pIds->GetNumberOfIds();
  int i;
  vtkIdType id;
  double sum = 0.0;
  weights->SetNumberOfTuples(numPts);
  double *w = weights->GetPointer(0);
  double y[3], v[3], r2, z2, rxy2, mag;
  double n[3], s;
  double f2=this->F2, e2=this->E2;

  for (i=0; i<numPts; ++i)
    {
    id = pIds->GetId(i);
    this->DataSet->GetPoint(id,y);

    v[0] = x[0] - y[0];
    v[1] = x[1] - y[1];
    v[2] = x[2] - y[2];
    r2 = vtkMath::Dot(v,v);

    if ( r2 == 0.0 ) //precise hit on existing point
      {
      pIds->SetNumberOfIds(1);
      pIds->SetId(0,id);
      weights->SetNumberOfTuples(1);
      weights->SetValue(0,1.0);
      return 1;
      }
    else // continue computing weights
      {
      // Normal affect
      if ( this->Normals )
        {
        this->Normals->GetTuple(id,n);
        mag = vtkMath::Dot(n,n);
        mag = ( mag == 0.0 ? 1.0 : sqrt(mag) );
        }
      else
        {
        mag = 1.0;
        }

      // Scalar scaling
      if ( this->Scalars )
        {
        this->Scalars->GetTuple(id,&s);
        }
      else
        {
        s = 1.0;
        }

      z2 = vtkMath::Dot(v,n) / mag;
      z2 = z2*z2;
      rxy2 = r2 - z2;

      w[i] = s * exp(-f2 * (rxy2/e2 + z2));
      sum += w[i];
      }//computing weights
    }//over all points

  // Normalize
  for (i=0; i<numPts; ++i)
    {
    w[i] /= sum;
    }

  return numPts;
}

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

  os << indent << "Radius: " << this->Radius << endl;
  os << indent << "Sharpness: " << this->Sharpness << endl;
  os << indent << "Eccentricity: " << this->Eccentricity << endl;
  os << indent << "Use Normals: "
     << (this->UseNormals? "On" : " Off") << "\n";
  os << indent << "Use Scalars: "
     << (this->UseScalars? "On" : " Off") << "\n";
}
