Skip to content

Commit

Permalink
+ optimize B-spline approximation
Browse files Browse the repository at this point in the history
  • Loading branch information
@wwmayer
wwmayer committed Nov 2, 2015
1 parent 2d4d60d commit 466d1a6
Show file tree
Hide file tree
Showing 4 changed files with 160 additions and 22 deletions.
4 changes: 4 additions & 0 deletions src/Mod/ReverseEngineering/App/AppReverseEngineering.cpp
View file
Expand Up @@ -160,6 +160,10 @@ class Module : public Py::ExtensionModule<Module>

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();
Expand Down
150 changes: 128 additions & 22 deletions src/Mod/ReverseEngineering/App/ApproxSurface.cpp
View file
Expand Up @@ -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"

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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++) {
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++;
}
}
}
}
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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);
Expand Down
26 changes: 26 additions & 0 deletions src/Mod/ReverseEngineering/App/ApproxSurface.h
View file
Expand Up @@ -45,6 +45,11 @@ namespace Reen {
class ReenExport SplineBasisfunction
{
public:
enum ValueT {
Zero = 0,
Full,
Other
};
/**
* Konstruktor
* @param iSize Length of Knots vector
Expand All @@ -71,6 +76,16 @@ class ReenExport SplineBasisfunction

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)
Expand Down Expand Up @@ -167,6 +182,17 @@ class ReenExport BSplineBasis : public SplineBasisfunction
*/
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)
Expand Down
2 changes: 2 additions & 0 deletions src/Mod/ReverseEngineering/App/CMakeLists.txt
View file
Expand Up @@ -18,6 +18,7 @@ include_directories(
${EIGEN3_INCLUDE_DIR}
${PCL_INCLUDE_DIRS}
${FLANN_INCLUDE_DIRS}
${QT_QTCORE_INCLUDE_DIR}
)

link_directories(${OCC_LIBRARY_DIR})
Expand All @@ -32,6 +33,7 @@ set(Reen_LIBS
${PCL_FEATURES_LIBRARIES}
${PCL_SEARCH_LIBRARIES}
${PCL_SURFACE_LIBRARIES}
${QT_QTCORE_LIBRARY}
)

SET(Reen_SRCS
Expand Down

0 comments on commit 466d1a6

Please sign in to comment.