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

  Program:   Visualization Toolkit
  Module:    $RCSfile: vtkGaussPointCloudSource.cxx,v $
  Language:  C++
  Date:      $Date: 2002/04/18 13:52:49 $
  Version:   $Revision: 1.4 $
  
  Made by Rasmus Paulsen
  email:  rrp@imm.dtu.dk
  web:    www.imm.dtu.dk/~rrp/VTK

  This class is not mature enough to enter the official VTK release.
=========================================================================*/
#include "vtkGaussPointCloudSource.h"
#include "vtkPoints.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"

vtkCxxRevisionMacro(vtkGaussPointCloudSource, "$Revision: 1.4 $");
vtkStandardNewMacro(vtkGaussPointCloudSource);

//----------------------------------------------------------------------------
vtkGaussPointCloudSource::vtkGaussPointCloudSource(float xL, float yL, float zL, int n)
{
	this->XSdev = fabs(xL);
	this->YSdev = fabs(yL);
	this->ZSdev = fabs(zL);

	this->Center[0] = 0.0;
	this->Center[1] = 0.0;
	this->Center[2] = 0.0;
	
	this->NPoint  = n;
	this->Seed  = 17;
}

//----------------------------------------------------------------------------

// The following routine is adapted from Numerical Recipes in C
double vtkGaussPointCloudSource::GaussRandomNumber()
{
	double rsq = 0;
	double v1 = 0;
	double v2 = 0;
	do
	{
		// Pick two uniform numbers in the square extending from -1 to +1 in each direction
		v1 = vtkMath::Random(-1, 1);;
		v2 = vtkMath::Random(-1, 1);;
		// See if they are in the unit circle, and if they are not, try again.
		rsq = v1*v1 + v2*v2;
	} 
	while (rsq >= 1.0 || rsq == 0);
	double fac = sqrt(-2.0 * log(rsq)/rsq);
	// Make the Box-Muller transform to get a normal deviate.

	return v1 * fac;
}

void vtkGaussPointCloudSource::Execute()
{
	vtkDebugMacro(<<"Creating Gaussian point cloud");
	
	vtkPolyData *output = this->GetOutput();
	
	vtkPoints *newPoints = vtkPoints::New();
	newPoints->Allocate(this->NPoint);
	
	vtkPoints *pts = vtkPoints::New();
	pts->Allocate(this->NPoint);
	
	vtkCellArray *verts = vtkCellArray::New();

	vtkMath::RandomSeed(this->Seed);

	for (int i = 0; i < this->NPoint; i++)
	{
		double x = GaussRandomNumber() * XSdev + Center[0];
		double y = GaussRandomNumber() * YSdev + Center[1];
		double z = GaussRandomNumber() * ZSdev + Center[2];

		vtkIdType id = pts->InsertNextPoint(x, y, z);
		verts->InsertNextCell(1);
		verts->InsertCellPoint(id);
	}
	
	output->SetPoints(pts);
	pts->Delete();
	output->SetVerts(verts);
	verts->Delete();
}

	
void vtkGaussPointCloudSource::PrintSelf(ostream& os, vtkIndent indent)
{
	vtkPolyDataSource::PrintSelf(os,indent);
	
	os << indent << "X Sdev: " << this->XSdev << "\n";
	os << indent << "Y Sdev: " << this->YSdev << "\n";
	os << indent << "Z Sdev: " << this->ZSdev << "\n";
	os << indent << "Center: (" << this->Center[0] << ", " 
		<< this->Center[1] << ", " << this->Center[2] << ")\n";
	os << indent << "NPoints: " << this->NPoint << "\n";
}