gbezier.cpp

一个非常有用的开源代码

/*
	Copyright (C) 2006, Mike Gashler

	This library is free software; you can redistribute it and/or
	modify it under the terms of the GNU Lesser General Public
	License as published by the Free Software Foundation; either
	version 2.1 of the License, or (at your option) any later version.

	see http://www.gnu.org/copyleft/lesser.html
*/

#include "GBezier.h"
#include "GRayTrace.h"
#include "GXML.h"
#include "GMacros.h"
#include <stdlib.h>

void Point3D::Transform(const struct Transform* pTransform)
{
	double c, s, t;

	// Size
	m_vals[0] *= pTransform->dScale;
	m_vals[1] *= pTransform->dScale;
	m_vals[2] *= pTransform->dScale;

	// Roll
//	c = pTransform->GetCosRoll();
//	s = pTransform->GetSinRoll();
//	t = x;
//	x = x * c + y * s;
//	y = y * c - t * s;

	// Pitch;
	c = pTransform->GetCosPitch();
	s = pTransform->GetSinPitch();
	t = m_vals[1];
	m_vals[1] = m_vals[1] * c + m_vals[2] * s;
	m_vals[2] = m_vals[2] * c - t * s;

	// Yaw
	c = pTransform->GetCosYaw();
	s = pTransform->GetSinYaw();
	t = m_vals[2];
	m_vals[2] = m_vals[2] * c + m_vals[0] * s;
	m_vals[0] = m_vals[0] * c - t * s;

	// Offset
	Add(&pTransform->offset);
}

void Point3D::Untransform(const struct Transform* pTransform)
{
	double c, s, t;

	// Offset
	Subtract(&pTransform->offset);

	// Yaw
	c = pTransform->GetCosYaw();
	s = pTransform->GetSinYaw();
	t = m_vals[2];
	m_vals[2] = m_vals[2] * c - m_vals[0] * s;
	m_vals[0] = m_vals[0] * c + t * s;

	// Pitch;
	c = pTransform->GetCosPitch();
	s = pTransform->GetSinPitch();
	t = m_vals[1];
	m_vals[1] = m_vals[1] * c - m_vals[2] * s;
	m_vals[2] = m_vals[2] * c + t * s;
	// Roll
//	c = pTransform->GetCosRoll();
//	s = pTransform->GetSinRoll();
//	t = x;
//	x = x * c - y * s;
//	y = y * c + t * s;

	// Size
	m_vals[0] /= pTransform->dScale;
	m_vals[1] /= pTransform->dScale;
	m_vals[2] /= pTransform->dScale;
}

void Point3D::FromXML(GXMLTag* pTag)
{
	GXMLAttribute* pAttr;
	pAttr = pTag->GetAttribute("x");
	if(pAttr)
		m_vals[0] = atof(pAttr->GetValue());
	pAttr = pTag->GetAttribute("y");
	if(pAttr)
		m_vals[1] = atof(pAttr->GetValue());
	pAttr = pTag->GetAttribute("z");
	if(pAttr)
		m_vals[2] = atof(pAttr->GetValue());
}

char* dtoa(double d, char* szBuff)
{
	sprintf(szBuff, "%f", d);
	
	// Remove trailing zeros
	int n;
	for(n = strlen(szBuff) - 1; n > 0 && szBuff[n] != '.'; n--)
	{
		if(szBuff[n] == '0')
			szBuff[n] = '\0';
		else
			break;
	}
	if(szBuff[n] == '.')
		szBuff[n] = '\0';

	return szBuff;
}

void Point3D::ToXML(GXMLTag* pTag)
{
	char szBuff[256];
	pTag->AddAttribute(new GXMLAttribute("x", dtoa(m_vals[0], szBuff)));
	pTag->AddAttribute(new GXMLAttribute("y", dtoa(m_vals[1], szBuff)));
	pTag->AddAttribute(new GXMLAttribute("z", dtoa(m_vals[2], szBuff)));
}


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



void Transform::FromXML(GXMLTag* pTag)
{
	offset.FromXML(pTag);
	GXMLAttribute* pAttr;
	pAttr = pTag->GetAttribute("Scale");
	if(pAttr)
		dScale = atof(pAttr->GetValue());
//	pAttr = pTag->GetAttribute("Roll");
//	if(pAttr)
//		dRoll = atof(pAttr->GetValue());
	pAttr = pTag->GetAttribute("Pitch");
	if(pAttr)
		dPitch = atof(pAttr->GetValue());
	pAttr = pTag->GetAttribute("Yaw");
	if(pAttr)
		dYaw = atof(pAttr->GetValue());
}

void Transform::ToXML(GXMLTag* pTag)
{
	char szBuff[256];
	offset.ToXML(pTag);
	if(dScale != 0)
		pTag->AddAttribute(new GXMLAttribute("Scale", dtoa(dScale, szBuff)));
//	if(dRoll != 0)
//		pTag->AddAttribute(new GXMLAttribute("Roll", dtoa(dScale, szBuff)));
	if(dPitch != 0)
		pTag->AddAttribute(new GXMLAttribute("Pitch", dtoa(dScale, szBuff)));
	if(dYaw != 0)
		pTag->AddAttribute(new GXMLAttribute("Yaw", dtoa(dScale, szBuff)));
}


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


struct GBezierPoint
{
	struct Point3D point;
	double weight;

	inline void Multiply(double dMag)
	{
		point.Multiply(dMag);
		weight *= dMag;
	}

	inline void Add(GBezierPoint* pPoint)
	{
		point.Add(&pPoint->point);
		weight += pPoint->weight;
	}
};

GBezier::GBezier(int nControlPoints)
{
	m_nControlPoints = nControlPoints;
	m_pPoints = new GBezierPoint[nControlPoints];
}

GBezier::GBezier(GBezier* pThat)
{
	m_nControlPoints = pThat->m_nControlPoints;
	m_pPoints = new GBezierPoint[m_nControlPoints];
	struct Point3D point;
	double weight;
	int n;
	for(n = 0; n < m_nControlPoints; n++)
	{
		pThat->GetControlPoint(&point, &weight, n);
		SetControlPoint(n, &point, weight);
	}
}

GBezier::~GBezier()
{
	delete(m_pPoints);
}

GBezier* GBezier::Copy()
{
	GBezier* pNewBezier = new GBezier(m_nControlPoints);
	struct Point3D point;
	double weight;
	int n;
	for(n = 0; n < m_nControlPoints; n++)
	{
		GetControlPoint(&point, &weight, n);
		pNewBezier->SetControlPoint(n, &point, weight);
	}
	return pNewBezier;
}

// Determine the point on the rational Bezier curve using the Horner algorithm
void GBezier::GetPoint(double t, struct Point3D* pOutPoint)
{
	int n = m_nControlPoints - 1;
	if(t == 1.0)
	{
		// Special case because the calcualtion for "u" below can't handle the case where t is 1
		*pOutPoint = m_pPoints[n].point;
		return;
	}
	double fw = m_pPoints[n].weight;
	pOutPoint->m_vals[0] = m_pPoints[n].point.m_vals[0] * fw;
	pOutPoint->m_vals[1] = m_pPoints[n].point.m_vals[1] * fw;
	pOutPoint->m_vals[2] = m_pPoints[n].point.m_vals[2] * fw;
	double u = t / (1.0 - t); // t/(1-t) in "monomialization trick" 
	double bc = 1.0; // binomial coefficient
	double den = 1.0; // denominator in P(t)/(1-t)^n
	double aux; // auxiliary variable - efficiency

	int i;
	for (i = n - 1; i >= 0; i--)
	{
		bc = (bc * (i + 1)) / (n - i); // new binomial coefficient 
		aux = bc * m_pPoints[i].weight;
		pOutPoint->m_vals[0] = pOutPoint->m_vals[0] * u + m_pPoints[i].point.m_vals[0] * aux;
		pOutPoint->m_vals[1] = pOutPoint->m_vals[1] * u + m_pPoints[i].point.m_vals[1] * aux;
		pOutPoint->m_vals[2] = pOutPoint->m_vals[2] * u + m_pPoints[i].point.m_vals[2] * aux;
		fw = fw * u + aux;
		den *= (1.0 - t);
	}

	pOutPoint->m_vals[0] /= fw;
	pOutPoint->m_vals[1] /= fw;
	pOutPoint->m_vals[2] /= fw;
}

int GBezier::GetControlPointCount()
{
	return m_nControlPoints;
}

void GBezier::GetControlPoint(struct Point3D* pOutPoint, double* pOutWeight, int n)
{
	GAssert(n >= 0 && n <= m_nControlPoints, "out of range");
	*pOutPoint = m_pPoints[n].point;
	*pOutWeight = m_pPoints[n].weight;
}

void GBezier::SetControlPoint(int n, struct Point3D* pPoint, double weight)
{
	GAssert(n >= 0 && n <= m_nControlPoints, "out of range");
	m_pPoints[n].point = *pPoint;
	m_pPoints[n].weight = weight;
}

void GBezier::ElevateDegree()
{
	// Allocate the new control points and set the first and last one
	GBezierPoint* pNewControlPoints = new GBezierPoint[m_nControlPoints + 1];
	pNewControlPoints[0] = m_pPoints[0];
	pNewControlPoints[m_nControlPoints] = m_pPoints[m_nControlPoints - 1];

	// Convert from weighted (rational) 3D control points to integral (non-rational) 4D control points
	int n;
	for(n = 0; n < m_nControlPoints; n++)
		m_pPoints[n].point.Multiply(m_pPoints[n].weight);

	// Determine the new intermediate control points
	struct GBezierPoint p1;
	struct GBezierPoint p2;
	for(n = 1; n < m_nControlPoints; n++)
	{
		// Find the new 4D control point
		double dFac = (double)n / (m_nControlPoints);
		p1 = m_pPoints[n - 1];
		p1.Multiply(dFac);
		p2 = m_pPoints[n];
		p2.Multiply(1.0 - dFac);
		pNewControlPoints[n] = p1;
		pNewControlPoints[n].Add(&p2);

		// Convert back to a weighted (rational) 3D point
		pNewControlPoints[n].point.Multiply(1.0 / pNewControlPoints[n].weight);
	}

	// Swap in the new control points
	delete(m_pPoints);
	m_pPoints = pNewControlPoints;
	m_nControlPoints++;
}

// Returns a 2D triangular array of 4D control points determined by the deCasteljau algorithm
struct GBezierPoint* GBezier::deCasteljau(double t, bool bTail)
{
	// Convert from weighted (rational) 3D control points to integral (non-rational) 4D control points
	int n;
	for(n = 0; n < m_nControlPoints; n++)
		m_pPoints[n].point.Multiply(m_pPoints[n].weight);

	// Populate the biggest row with the current control points
	GBezierPoint* pDeCasteljauPoints = (GBezierPoint*)alloca(sizeof(GBezierPoint) * m_nControlPoints * (m_nControlPoints + 1) / 2);
	int y = m_nControlPoints - 1;
	n = (y * (y + 1)) / 2;
	int x;
	for(x = 0; x < m_nControlPoints; x++)
		pDeCasteljauPoints[n + x] = m_pPoints[x];

	// Calculate the rest of the points
	struct GBezierPoint p1;
	struct GBezierPoint p2;
	for(y--; y >= 0; y--)
	{
		n = (y * (y + 1)) / 2;
		for(x = 0; x <= y; x++)
		{
			p1 = pDeCasteljauPoints[n + x + y + 1];
			p2 = pDeCasteljauPoints[n + x + y + 2];
			p1.Multiply(1.0 - t);
			p2.Multiply(t);
			pDeCasteljauPoints[n + x] = p1;
			pDeCasteljauPoints[n + x].Add(&p2);
		}
	}

	// Extract the segment control points and convert back to weighted (rational) 3D points
	GBezierPoint* pSegment = new GBezierPoint[m_nControlPoints];
	for(y = 0; y < m_nControlPoints; y++)
	{
		if(bTail)
			n = (y * (y + 1)) / 2 + y;
		else
			n = ((m_nControlPoints - 1 - y) * (m_nControlPoints - y)) / 2;
		pSegment[y] = pDeCasteljauPoints[n];
		pSegment[y].point.Multiply(1 / pSegment[y].weight);
	}
	return pSegment;