suzanne.soy/FreeCAD
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

+ optimize B-spline approximation

Browse Source
This commit is contained in:
wmayer 2015-11-03 00:55:15 +01:00
parent 2d4d60d2ab
commit 466d1a60b1
4 changed files with 160 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/Mod/ReverseEngineering/App/AppReverseEngineering.cpp
View File

@ -160,6 +160,10 @@ private:
throw Py::RuntimeError("Computation of B-Spline surface failed");
}
catch (const Py::Exception&) {
// re-throw
throw;
}
catch (Standard_Failure &e) {
std::string str;
Standard_CString msg = e.GetMessageString();

150
src/Mod/ReverseEngineering/App/ApproxSurface.cpp
View File

@ -26,9 +26,15 @@
#include <math_Householder.hxx>
#include <Geom_BSplineSurface.hxx>
#include <QFuture>
#include <QFutureWatcher>
#include <QtConcurrentMap>
#include <boost/bind.hpp>
#include <Mod/Mesh/App/Core/Approximation.h>
#include <Base/Sequencer.h>
#include <Base/Tools2D.h>
#include <Base/Tools.h>
#include "ApproxSurface.h"
@ -174,6 +180,23 @@ void BSplineBasis::AllBasisFunctions(double fParam, TColStd_Array1OfReal& vFuncV
}
}
BSplineBasis::ValueT BSplineBasis::LocalSupport(int iIndex, double fParam)
{
int m = _vKnotVector.Length()-1;
int p = _iOrder-1;
if ((iIndex == 0 && fParam == _vKnotVector(0)) ||
(iIndex == m-p-1 && fParam == _vKnotVector(m))) {
return BSplineBasis::Full;
}
if (fParam < _vKnotVector(iIndex) || fParam >= _vKnotVector(iIndex+p+1)) {
return BSplineBasis::Zero;
}
return BSplineBasis::Other;
}
double BSplineBasis::BasisFunction(int iIndex, double fParam)
{
int m = _vKnotVector.Length()-1;
@ -889,10 +912,29 @@ bool BSplineParameterCorrection::SolveWithoutSmoothing()
double fV = (*_pvcUVParam)(i).Y();
unsigned long ulIdx=0;
// Vorberechnung der Werte der Basis-Funktionen
std::vector<double> basisU(_usUCtrlpoints);
for (unsigned short j=0; j<_usUCtrlpoints; j++) {
for (unsigned short k=0; k<_usVCtrlpoints; k++) {
M(i,ulIdx) = _clUSpline.BasisFunction(j,fU)*_clVSpline.BasisFunction(k,fV);
ulIdx++;
basisU[j] = _clUSpline.BasisFunction(j,fU);
}
std::vector<double> basisV(_usVCtrlpoints);
for (unsigned short k=0; k<_usVCtrlpoints; k++) {
basisV[k] = _clVSpline.BasisFunction(k,fV);
}
for (unsigned short j=0; j<_usUCtrlpoints; j++) {
double valueU = basisU[j];
if (valueU == 0.0) {
for (unsigned short k=0; k<_usVCtrlpoints; k++) {
M(i,ulIdx) = 0.0;
ulIdx++;
}
}
else {
for (unsigned short k=0; k<_usVCtrlpoints; k++) {
M(i,ulIdx) = valueU * basisV[k];
ulIdx++;
}
}
}
}
@ -925,52 +967,116 @@ bool BSplineParameterCorrection::SolveWithoutSmoothing()
return true;
}
namespace Reen {
class ScalarProduct
{
public:
ScalarProduct(const math_Matrix& mat) : mat(mat)
{
}
std::vector<double> multiply(int col) const
{
math_Vector vec = mat.Col(col);
std::vector<double> out(mat.ColNumber());
for (int n=mat.LowerCol(); n<=mat.UpperCol(); n++) {
out[n] = vec * mat.Col(n);
}
return out;
}
private:
const math_Matrix& mat;
};
}
bool BSplineParameterCorrection::SolveWithSmoothing(double fWeight)
{
unsigned long ulSize = _pvcPoints->Length();
unsigned long ulDim = _usUCtrlpoints*_usVCtrlpoints;
math_Matrix M (0, ulSize-1, 0,_usUCtrlpoints*_usVCtrlpoints-1);
math_Matrix MTM(0, _usUCtrlpoints*_usVCtrlpoints-1,
0, _usUCtrlpoints*_usVCtrlpoints-1);
math_Vector Xx (0, _usUCtrlpoints*_usVCtrlpoints-1);
math_Vector Xy (0, _usUCtrlpoints*_usVCtrlpoints-1);
math_Vector Xz (0, _usUCtrlpoints*_usVCtrlpoints-1);
math_Matrix M (0, ulSize-1, 0, ulDim-1);
math_Vector Xx (0, ulDim-1);
math_Vector Xy (0, ulDim-1);
math_Vector Xz (0, ulDim-1);
math_Vector bx (0, ulSize-1);
math_Vector by (0, ulSize-1);
math_Vector bz (0, ulSize-1);
math_Vector Mbx(0, _usUCtrlpoints*_usVCtrlpoints-1);
math_Vector Mby(0, _usUCtrlpoints*_usVCtrlpoints-1);
math_Vector Mbz(0, _usUCtrlpoints*_usVCtrlpoints-1);
math_Vector Mbx(0, ulDim-1);
math_Vector Mby(0, ulDim-1);
math_Vector Mbz(0, ulDim-1);
//Bestimmung der Koeffizientenmatrix des ueberbestimmten LGS
for (unsigned long i=0; i<ulSize; i++) {
double fU = (*_pvcUVParam)(i).X();
double fV = (*_pvcUVParam)(i).Y();
unsigned long ulIdx=0;
int ulIdx=0;
// Vorberechnung der Werte der Basis-Funktionen
std::vector<double> basisU(_usUCtrlpoints);
for (unsigned short j=0; j<_usUCtrlpoints; j++) {
basisU[j] = _clUSpline.BasisFunction(j,fU);
}
std::vector<double> basisV(_usVCtrlpoints);
for (unsigned short k=0; k<_usVCtrlpoints; k++) {
basisV[k] = _clVSpline.BasisFunction(k,fV);
}
for (unsigned short j=0; j<_usUCtrlpoints; j++) {
for (unsigned short k=0; k<_usVCtrlpoints; k++) {
M(i,ulIdx) = _clUSpline.BasisFunction(j,fU)*_clVSpline.BasisFunction(k,fV);
ulIdx++;
double valueU = basisU[j];
if (valueU == 0.0) {
for (unsigned short k=0; k<_usVCtrlpoints; k++) {
M(i,ulIdx) = 0.0;
ulIdx++;
}
}
else {
for (unsigned short k=0; k<_usVCtrlpoints; k++) {
M(i,ulIdx) = valueU * basisV[k];
ulIdx++;
}
}
}
}
//Das Produkt aus ihrer Transformierten und ihr selbst ergibt die quadratische Systemmatrix
//Das Produkt aus ihrer Transformierten und ihr selbst ergibt die quadratische Systemmatrix (langsam)
#if 0
math_Matrix MTM = M.TMultiply(M);
#elif 0
math_Matrix MTM(0, ulDim-1, 0, ulDim-1);
for (unsigned long m=0; m<ulDim; m++) {
math_Vector Mm = M.Col(m);
for (unsigned long n=m; n<ulDim; n++) {
MTM(m,n)=MTM(n,m)=M.Col(m)*M.Col(n);
MTM(m,n) = MTM(n,m) = Mm * M.Col(n);
}
}
#else // multi-threaded
std::vector<int> columns(ulDim);
std::generate(columns.begin(), columns.end(), Base::iotaGen<int>(0));
ScalarProduct scalar(M);
QFuture< std::vector<double> > future = QtConcurrent::mapped
(columns, boost::bind(&ScalarProduct::multiply, &scalar, _1));
QFutureWatcher< std::vector<double> > watcher;
watcher.setFuture(future);
watcher.waitForFinished();
math_Matrix MTM(0, ulDim-1, 0, ulDim-1);
int rowIndex=0;
for (QFuture< std::vector<double> >::const_iterator it = future.begin(); it != future.end(); ++it) {
int colIndex=0;
for (std::vector<double>::const_iterator jt = it->begin(); jt != it->end(); ++jt, colIndex++)
MTM(rowIndex, colIndex) = *jt;
rowIndex++;
}
#endif
//Bestimmen der rechten Seite
for (int ii=_pvcPoints->Lower(); ii<=_pvcPoints->Upper(); ii++) {
bx(ii) = (*_pvcPoints)(ii).X(); by(ii) = (*_pvcPoints)(ii).Y(); bz(ii) = (*_pvcPoints)(ii).Z();
}
for (unsigned long i=0; i<ulDim; i++) {
Mbx(i) = M.Col(i) * bx;
Mby(i) = M.Col(i) * by;
Mbz(i) = M.Col(i) * bz;
math_Vector Mi = M.Col(i);
Mbx(i) = Mi * bx;
Mby(i) = Mi * by;
Mbz(i) = Mi * bz;
}
// Loese das LGS mit der LU-Zerlegung
@ -1048,7 +1154,7 @@ void BSplineParameterCorrection::CalcSecondSmoothMatrix(Base::SequencerLauncher&
for (unsigned short j=0; j<_usVCtrlpoints; j++) {
_clSecondMatrix(m,n) = _clUSpline.GetIntegralOfProductOfBSplines(i,k,2,2) *
_clVSpline.GetIntegralOfProductOfBSplines(j,l,0,0) +
2*_clUSpline.GetIntegralOfProductOfBSplines(i,k,1,1) *
2*_clUSpline.GetIntegralOfProductOfBSplines(i,k,1,1) *
_clVSpline.GetIntegralOfProductOfBSplines(j,l,1,1) +
_clUSpline.GetIntegralOfProductOfBSplines(i,k,0,0) *
_clVSpline.GetIntegralOfProductOfBSplines(j,l,2,2);

26
src/Mod/ReverseEngineering/App/ApproxSurface.h
View File

@ -45,6 +45,11 @@ namespace Reen {
class ReenExport SplineBasisfunction
{
public:
enum ValueT {
Zero = 0,
Full,
Other
};
/**
* Konstruktor
* @param iSize Length of Knots vector
@ -71,6 +76,16 @@ public:
virtual ~SplineBasisfunction();
/**
* Gibt an, ob der Funktionswert Nik(t) an der Stelle fParam
* 0, 1 oder ein Wert dazwischen ergibt.
* Dies dient dazu, um die Berechnung zu u.U. zu beschleunigen.
*
* @param iIndex Index
* @param fParam Parameterwert
* @return ValueT
*/
virtual ValueT LocalSupport(int iIndex, double fParam)=0;
/**
* Berechnet den Funktionswert Nik(t) an der Stelle fParam
* (aus: Piegl/Tiller 96 The NURBS-Book)
@ -167,6 +182,17 @@ public:
*/
virtual void AllBasisFunctions(double fParam, TColStd_Array1OfReal& vFuncVals);
/**
* Gibt an, ob der Funktionswert Nik(t) an der Stelle fParam
* 0, 1 oder ein Wert dazwischen ergibt.
* Dies dient dazu, um die Berechnung zu u.U. zu beschleunigen.
*
* @param iIndex Index
* @param fParam Parameterwert
* @return ValueT
*/
virtual ValueT LocalSupport(int iIndex, double fParam);
/**
* Berechnet den Funktionswert Nik(t) an der Stelle fParam
* (aus: Piegl/Tiller 96 The NURBS-Book)

2
src/Mod/ReverseEngineering/App/CMakeLists.txt
View File

@ -18,6 +18,7 @@ include_directories(
${EIGEN3_INCLUDE_DIR}
${PCL_INCLUDE_DIRS}
${FLANN_INCLUDE_DIRS}
${QT_QTCORE_INCLUDE_DIR}
)
link_directories(${OCC_LIBRARY_DIR})
@ -32,6 +33,7 @@ set(Reen_LIBS
${PCL_FEATURES_LIBRARIES}
${PCL_SEARCH_LIBRARIES}
${PCL_SURFACE_LIBRARIES}
${QT_QTCORE_LIBRARY}
)
SET(Reen_SRCS