/*========================================================================= 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"; }