Permalink
Fetching contributors…
Cannot retrieve contributors at this time
171 lines (141 sloc) 4.97 KB
/* ********************************************************************
* PointAreaComputer
*
* Copyright 2012 Joachim Giard
* Université catholique de Louvain, Belgium
* Apache License, Version 2.0
*
* Computes the voronoi area of each point
*
* Arno Klein added lines to avoid NaN errors in 2017.
*
* *******************************************************************/
#include "PointAreaComputer.h"
#include <vtkCellArray.h>
#include <vtkDoubleArray.h>
#include <vtkMath.h>
#include <vtkPolyDataWriter.h>
#include <vtkPointData.h>
PointAreaComputer::PointAreaComputer(vtkPolyData *mesh)
{
m_mesh = mesh;
m_pointsArea = vtkDoubleArray::New();
}
PointAreaComputer::~PointAreaComputer()
{
m_pointsArea->Delete();
}
void PointAreaComputer::ComputeArea()
{
m_pointsArea = ComputeVoronoiArea();
}
vtkDoubleArray* PointAreaComputer::GetArea()
{
return m_pointsArea;
}
void PointAreaComputer::WriteIntoFile(char *fileName)
{
vtkPolyDataWriter* writer=vtkPolyDataWriter::New();
writer->SetFileName(fileName);
m_mesh->GetPointData()->SetScalars(m_pointsArea);
writer->SetInputData(m_mesh);
writer->Write();
writer->Delete();
}
vtkDoubleArray* PointAreaComputer::ComputeVoronoiArea()
{
// initialisation
vtkCellArray* cells=m_mesh->GetPolys();
int nbPolys = cells->GetNumberOfCells();
vtkDoubleArray* pointsArea = vtkDoubleArray::New();
for(int i=0;i<m_mesh->GetNumberOfPoints();i++)
{
pointsArea->InsertNextValue(0);
}
int cellIds[3];
double pos1[3];
double pos2[3];
double pos3[3];
float a;
float b;
float c;
float alpha;
float beta;
float gamma;
float area, areaA, areaB, areaC;
bool obtuse = false;
for (int i=0;i<nbPolys;i++)
{
cellIds[0]=m_mesh->GetCell(i)->GetPointId(0);
cellIds[1]=m_mesh->GetCell(i)->GetPointId(1);
cellIds[2]=m_mesh->GetCell(i)->GetPointId(2);
if(cellIds[0]>m_mesh->GetNumberOfPoints()||cellIds[1]>m_mesh->GetNumberOfPoints()||cellIds[2]>m_mesh->GetNumberOfPoints())
{
cout<<"point surface cell id problem: "<<i<<" "<<cellIds[0]<<" "<<cellIds[1]<<" "<<cellIds[2]<<endl;
}
else
{
m_mesh->GetPoint(cellIds[0],pos1);
m_mesh->GetPoint(cellIds[1],pos2);
m_mesh->GetPoint(cellIds[2],pos3);
// distances between points of each triangle
a=sqrt(vtkMath::Distance2BetweenPoints(pos3,pos2));
b=sqrt(vtkMath::Distance2BetweenPoints(pos1,pos3));
c=sqrt(vtkMath::Distance2BetweenPoints(pos2,pos1));
// set 1/3 areas to zero by default
areaA = 0;
areaB = 0;
areaC = 0;
// if distances greater than zero
if(a>0 && b>0 && c>0)
{
// angles
alpha = acos( (pow(b,2) + pow(c,2) - pow(a,2)) / (2*b*c) );
beta = acos( (pow(a,2) + pow(c,2) - pow(b,2)) / (2*a*c) );
gamma = acos( (pow(b,2) + pow(a,2) - pow(c,2)) / (2*b*a) );
// don't set 1/3 areas if angles are not numbers
if( vtkMath::IsNan(alpha) || vtkMath::IsNan(beta) || vtkMath::IsNan(gamma) )
{
}
// else set 1/3 areas
else
{
if( alpha > vtkMath::Pi()/2 || beta > vtkMath::Pi()/2 || gamma > vtkMath::Pi()/2 )
{
// Herons's formula
area = 0.25 * sqrt( (a+b+c)*(b+c-a)*(a-b+c)*(a+b-c) );
if(alpha > vtkMath::Pi()/2)
{
areaA = area/2;
areaB = area/4;
areaC = area/4;
}
else if (beta > vtkMath::Pi()/2)
{
areaA = area/4;
areaB = area/2;
areaC = area/4;
}
else if (gamma > vtkMath::Pi()/2)
{
areaA = area/4;
areaB = area/4;
areaC = area/2;
}
}
else
{
areaA = (pow(b,2) / tan(beta) + pow(c,2) / tan(gamma)) /8;
areaB = (pow(a,2) / tan(alpha) + pow(c,2) / tan(gamma)) /8;
areaC = (pow(b,2) / tan(beta) + pow(a,2) / tan(alpha)) /8;
}
}
}
// add 1/3 of the face area to each point of the triangle
pointsArea->SetValue(cellIds[0],pointsArea->GetValue(cellIds[0])+ areaA);
pointsArea->SetValue(cellIds[1],pointsArea->GetValue(cellIds[1])+ areaB);
pointsArea->SetValue(cellIds[2],pointsArea->GetValue(cellIds[2])+ areaC);
}
}
return pointsArea;
}