[vtkusers] question about vtkMeshQuality.cxx

Tue Jul 20 09:49:45 EDT 2004

Indeed , vtkMeshQuality was not enough for us too.

The problem with the "quality" of a cell is that many measures exist in the 
litterature depending on the application.
As philippe seems interested in making something additions, I implemented 
one of the triangle shape measures existing :
1/2 perimeter * longest edge / area
normalized by the value for an equilateral triangle.
It measure the "stretching factor" of a triangle.

please find the code attached.
it's not perfect, some checking steps were never implemented, but in those 
cases, a comment line exists to remind that a check should be made.
If it can be of any help, use it.

alex.

A 11:05 19/07/2004 -0700, Philippe P. Pebay a ecrit :
>Hi
>
>I was looking at using vtkMeshQuality for some paper results, but it seems 
>that the name of this class is too generic, because:
>1. it only addresses tetrahedral meshes
>2. in the case of tetrahedral meshes, there are many quality measures that 
>can be used, but this class only provides the "radius-ratio" (normalized 
>circumradius / inradius).
>
>I haven't seen any other quality metrics in VTK; would anyone mind if we 
>revise vtkMeshQuality to compute various quality measures on several 
>different types of cells ?  We would like to compute the following quality 
>measures:
>1. tetrahedral meshes: radius-ratio, aspect-ratio, shape measures 
>(Frobenius-like), etc.
>2. hexahedral meshes
>3. triangular and quadrilateral elements: various measures in the literature.
>
>Also, I noticed that the Execute method in vtkMeshQuality.cxx declares p# 
>and dp# (#= 1,2,3,4) as double arrays, and then copies p's into dp's. But 
>it appears that dp's are never altered, so why do 2 copies of the points 
>exist ?
>
>Cheers
>--
>===============================================================
>Philippe P. Pebay             |  pppebay at ca.sandia.gov
>Sandia National Laboratories  |  Phone: +1 925 294 2024
>PO Box 969, MS 9051           |  Cell : +1 925 784 3255
>Livermore, CA 94550 U.S.A.    |  Fax  : +1 925 294 2595
>            www.ca.sandia.gov/CRF/staff/pebay.html
>===============================================================
>"Success is the ability to go from one failure
>to another with no loss of enthusiasm."
>-- Sir Winston Churchill
>
>
>_______________________________________________
>This is the private VTK discussion list. Please keep messages on-topic. 
>Check the FAQ at: <http://public.kitware.com/cgi-bin/vtkfaq>
>Follow this link to subscribe/unsubscribe:
>http://www.vtk.org/mailman/listinfo/vtkusers
-------------- next part --------------
/*=========================================================================

  Program:   Mailleur 3D multi-r?solution (Creatis PFE  2001~)
  Module:    vtkTriangleShapeFactor.h
  Language:  C++
  Date:      2003/03
  Auteurs:   ALex. GOUAILLARD

=========================================================================*/
// .NAME vtkTriangleShapeFactor 
// .SECTION Description
// This class compute the mean shape factor 
// of all the triangles of a given PolyData

#ifndef __vtkTriangleShapeFactor_h
#define __vtkTriangleShapeFactor_h

#include <vtkObject.h>
#include <vtkPolyData.h>
#include <vtkMath.h>

class VTK_EXPORT vtkTriangleShapeFactor : public vtkObject 
{

public:

  // VTK Constructor 
  static vtkTriangleShapeFactor *New();

  // VTK inheritance scheme
  vtkTypeMacro(vtkTriangleShapeFactor,vtkObject);

  // To get the result as a float
  float GetMeanShapeFactor() {return this->MeanShapeFactor;};

  //To get thie min angle of 
  static double GetMinAngle(vtkPolyData *input,int IdCell);

  // Set a PolyData as input
  void SetInput(vtkPolyData *UserSetInput) {this->Input = UserSetInput;};  

  // compute the value and put it in this->MeanShapeFactor
  void Update();

protected:
  vtkTriangleShapeFactor(); //constructeur
  ~vtkTriangleShapeFactor() {}; // destructeur
  vtkTriangleShapeFactor(const vtkTriangleShapeFactor&) {};
  void operator=(const vtkTriangleShapeFactor&) {};

  float MeanShapeFactor;

  vtkPolyData *Input;

};

#endif
-------------- next part --------------
/*=========================================================================

  Program:   Mailleur 3D multi-r?solution (Creatis PFE  2001~)
  Module:    vtkTriangleShapeFactor.cxx
  Language:  C++
  Date:      2003/03
  Auteurs:   Alex. GOUAILLARD

=========================================================================*/
#include <vtkObjectFactory.h>
#include <vtkMath.h>
#include <vtkTriangle.h>

#include "vtkTriangleShapeFactor.h"
#include "vtkParameterization.h"
//--------------------------------------------------------------------------

vtkTriangleShapeFactor* vtkTriangleShapeFactor::New()
{
  // First try to create the object from the vtkObjectFactory
  vtkObject* ret = vtkObjectFactory::CreateInstance("vtkTriangleShapeFactor");
  if(ret)
    {
    return (vtkTriangleShapeFactor*)ret;
    }
  // If the factory was unable to create the object, then create it here.
  return new vtkTriangleShapeFactor;
}

vtkTriangleShapeFactor::vtkTriangleShapeFactor()
{
	this->Input = NULL;
	this->MeanShapeFactor = 0.0;
}

void vtkTriangleShapeFactor::Update() 
{
	// check inputs
	if (!this->Input) return;

	// check if it's a triangular mesh

	// name variables
	double QualitySum;
	double RegularTriangleShapeFactor;
	double LongestEdgeLength;
	double TempLength;
	double Perimeter;

	int NbCells;
	int IdVertex1;
	int IdVertex2;
	int IdVertex3;

	register int i;

	double Vertex1[3];
	double Vertex2[3];
	double Vertex3[3];

	// initialization
	NbCells = this->Input->GetNumberOfCells();
	RegularTriangleShapeFactor = 2.0 * sqrt(3.0);
	QualitySum = 0.0;

	// loop on all the cells and compute the unnormalyzed shape factor
	for (i=0;i<NbCells;i++)
	{
		// Initialyze
		LongestEdgeLength = 0.0;
		Perimeter = 0.0;

		// get the points coordinates
		vtkParameterization::GetFaceVertices(this->Input,i,IdVertex1,IdVertex2,IdVertex3);
		this->Input->GetPoint(IdVertex1,Vertex1);
		this->Input->GetPoint(IdVertex2,Vertex2);
		this->Input->GetPoint(IdVertex3,Vertex3);

		// compute the length of each edges
		TempLength = sqrt(vtkMath::Distance2BetweenPoints(Vertex1,Vertex2));
		if(TempLength > LongestEdgeLength) LongestEdgeLength = TempLength;
		Perimeter += TempLength;

		TempLength = sqrt(vtkMath::Distance2BetweenPoints(Vertex1,Vertex3));
		if(TempLength > LongestEdgeLength) LongestEdgeLength = TempLength;
		Perimeter += TempLength;

		TempLength = sqrt(vtkMath::Distance2BetweenPoints(Vertex3,Vertex2));
		if(TempLength > LongestEdgeLength) LongestEdgeLength = TempLength;
		Perimeter += TempLength;

		// shape factor = half perimeter * longest edge / area.
		QualitySum += 0.5 * Perimeter * LongestEdgeLength / vtkTriangle::TriangleArea(Vertex1,Vertex2,Vertex3);

	}

	// normalize
	this->MeanShapeFactor = NbCells * RegularTriangleShapeFactor / QualitySum;

}

double vtkTriangleShapeFactor::GetMinAngle(vtkPolyData *input,int IdCell)
{

	register int i;

	int Id[3];

	double A[3];
	double B[3];
	double C[3];

	double BA[3];
	double BC[3];
	double AC[3];

	double dot[3];
	double cross[3][3];
	double tan_teta[3];

	double min;

	for(i=0;i<3;i++)
		tan_teta[i] = -1;

	vtkParameterization::GetFaceVertices(input,IdCell,Id[0],Id[1],Id[2]);
	input->GetPoint(Id[0],A);
	input->GetPoint(Id[1],B);
	input->GetPoint(Id[2],C);

	for(i=0;i<3;i++)
	{
		BA[i] = A[i] - B[i];
		BC[i] = C[i] - B[i];
		AC[i] = C[i] - A[i];
	}

	dot[0] = -1.0 * vtkMath::Dot(BA,AC);
	dot[1] = vtkMath::Dot(BA,BC);
	dot[2] = vtkMath::Dot(AC,BC);

	vtkMath::Cross(BA,AC,cross[0]);
	vtkMath::Cross(BA,BC,cross[1]);
	vtkMath::Cross(AC,BC,cross[2]);

	for(i=0;i<3;i++)
	{
		if(dot[i] < 0)
			tan_teta[i] = -1.0 * vtkMath::Norm(cross[i]) / dot[i];
		else if(dot[i] > 0)
			tan_teta[i] = vtkMath::Norm(cross[i]) / dot[i];
	}

	if(tan_teta[0] != -1) //if(teta[0] != pi/2)
		min = tan_teta[0];

	else 
		min = tan_teta[1];

	for(i=0;i<3;i++)
		if((tan_teta[i] < min) && (tan_teta[i] > 0))
			min = tan_teta[i];

	return(min);
}