f523ab6ee3
"Professional CMake" book suggest the following: "Targets should build successfully with or without compiler support for precompiled headers. It should be considered an optimization, not a requirement. In particular, do not explicitly include a precompile header (e.g. stdafx.h) in the source code, let CMake force-include an automatically generated precompile header on the compiler command line instead. This is more portable across the major compilers and is likely to be easier to maintain. It will also avoid warnings being generated from certain code checking tools like iwyu (include what you use)." Therefore, removed the "#include <PreCompiled.h>" from sources, also there is no need for the "#ifdef _PreComp_" anymore
1417 lines
45 KiB
C++
1417 lines
45 KiB
C++
/***************************************************************************
|
|
* Copyright (c) 2008 Werner Mayer <wmayer[at]users.sourceforge.net> *
|
|
* *
|
|
* This file is part of the FreeCAD CAx development system. *
|
|
* *
|
|
* This library is free software; you can redistribute it and/or *
|
|
* modify it under the terms of the GNU Library General Public *
|
|
* License as published by the Free Software Foundation; either *
|
|
* version 2 of the License, or (at your option) any later version. *
|
|
* *
|
|
* This library is distributed in the hope that it will be useful, *
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
|
* GNU Library General Public License for more details. *
|
|
* *
|
|
* You should have received a copy of the GNU Library General Public *
|
|
* License along with this library; see the file COPYING.LIB. If not, *
|
|
* write to the Free Software Foundation, Inc., 59 Temple Place, *
|
|
* Suite 330, Boston, MA 02111-1307, USA *
|
|
* *
|
|
***************************************************************************/
|
|
|
|
#include <QFuture>
|
|
#include <QFutureWatcher>
|
|
#include <QtConcurrentMap>
|
|
|
|
#include <Geom_BSplineSurface.hxx>
|
|
#include <Precision.hxx>
|
|
#include <math_Gauss.hxx>
|
|
#include <math_Householder.hxx>
|
|
|
|
#include <Base/Sequencer.h>
|
|
#include <Base/Tools.h>
|
|
#include <Mod/Mesh/App/Core/Approximation.h>
|
|
|
|
#include "ApproxSurface.h"
|
|
|
|
|
|
using namespace Reen;
|
|
namespace sp = std::placeholders;
|
|
|
|
// SplineBasisfunction
|
|
|
|
SplineBasisfunction::SplineBasisfunction(int iSize)
|
|
: _vKnotVector(0, iSize - 1)
|
|
, _iOrder(1)
|
|
{}
|
|
|
|
SplineBasisfunction::SplineBasisfunction(TColStd_Array1OfReal& vKnots,
|
|
TColStd_Array1OfInteger& vMults,
|
|
int iSize,
|
|
int iOrder)
|
|
: _vKnotVector(0, iSize - 1)
|
|
{
|
|
int sum = 0;
|
|
for (int h = vMults.Lower(); h <= vMults.Upper(); h++) {
|
|
sum += vMults(h);
|
|
}
|
|
|
|
if (vKnots.Length() != vMults.Length() || iSize != sum) {
|
|
// Throw exception
|
|
Standard_ConstructionError::Raise("BSplineBasis");
|
|
}
|
|
|
|
int k = 0;
|
|
for (int i = vMults.Lower(); i <= vMults.Upper(); i++) {
|
|
for (int j = 0; j < vMults(i); j++) {
|
|
_vKnotVector(k) = vKnots(i);
|
|
k++;
|
|
}
|
|
}
|
|
|
|
_iOrder = iOrder;
|
|
}
|
|
|
|
SplineBasisfunction::SplineBasisfunction(TColStd_Array1OfReal& vKnots, int iOrder)
|
|
: _vKnotVector(0, vKnots.Length() - 1)
|
|
{
|
|
_vKnotVector = vKnots;
|
|
_iOrder = iOrder;
|
|
}
|
|
|
|
SplineBasisfunction::~SplineBasisfunction() = default;
|
|
|
|
void SplineBasisfunction::SetKnots(TColStd_Array1OfReal& vKnots, int iOrder)
|
|
{
|
|
if (_vKnotVector.Length() != vKnots.Length()) {
|
|
Standard_RangeError::Raise("BSplineBasis");
|
|
}
|
|
|
|
_vKnotVector = vKnots;
|
|
_iOrder = iOrder;
|
|
}
|
|
|
|
void SplineBasisfunction::SetKnots(TColStd_Array1OfReal& vKnots,
|
|
TColStd_Array1OfInteger& vMults,
|
|
int iOrder)
|
|
{
|
|
int sum = 0;
|
|
for (int h = vMults.Lower(); h <= vMults.Upper(); h++) {
|
|
sum += vMults(h);
|
|
}
|
|
|
|
if (vKnots.Length() != vMults.Length() || _vKnotVector.Length() != sum) {
|
|
// Throw exception
|
|
Standard_RangeError::Raise("BSplineBasis");
|
|
}
|
|
int k = 0;
|
|
for (int i = vMults.Lower(); i <= vMults.Upper(); i++) {
|
|
for (int j = 0; j < vMults(i); j++) {
|
|
_vKnotVector(k) = vKnots(i);
|
|
k++;
|
|
}
|
|
}
|
|
|
|
_iOrder = iOrder;
|
|
}
|
|
|
|
////////////////////////////////////////// BSplineBasis
|
|
|
|
BSplineBasis::BSplineBasis(int iSize)
|
|
: SplineBasisfunction(iSize)
|
|
{}
|
|
|
|
BSplineBasis::BSplineBasis(TColStd_Array1OfReal& vKnots,
|
|
TColStd_Array1OfInteger& vMults,
|
|
int iSize,
|
|
int iOrder)
|
|
: SplineBasisfunction(vKnots, vMults, iSize, iOrder)
|
|
{}
|
|
|
|
BSplineBasis::BSplineBasis(TColStd_Array1OfReal& vKnots, int iOrder)
|
|
: SplineBasisfunction(vKnots, iOrder)
|
|
{}
|
|
|
|
BSplineBasis::~BSplineBasis() = default;
|
|
|
|
int BSplineBasis::FindSpan(double fParam)
|
|
{
|
|
int n = _vKnotVector.Length() - _iOrder - 1;
|
|
if (fParam == _vKnotVector(n + 1)) {
|
|
return n;
|
|
}
|
|
|
|
int low = _iOrder - 1;
|
|
int high = n + 1;
|
|
int mid = (low + high) / 2; // Binary search
|
|
|
|
while (fParam < _vKnotVector(mid) || fParam >= _vKnotVector(mid + 1)) {
|
|
if (fParam < _vKnotVector(mid)) {
|
|
high = mid;
|
|
}
|
|
else {
|
|
low = mid;
|
|
}
|
|
mid = (low + high) / 2;
|
|
}
|
|
|
|
return mid;
|
|
}
|
|
|
|
void BSplineBasis::AllBasisFunctions(double fParam, TColStd_Array1OfReal& vFuncVals)
|
|
{
|
|
if (vFuncVals.Length() != _iOrder) {
|
|
Standard_RangeError::Raise("BSplineBasis");
|
|
}
|
|
|
|
int iIndex = FindSpan(fParam);
|
|
|
|
TColStd_Array1OfReal zaehler_left(1, _iOrder - 1);
|
|
TColStd_Array1OfReal zaehler_right(1, _iOrder - 1);
|
|
vFuncVals(0) = 1.0;
|
|
|
|
for (int j = 1; j < _iOrder; j++) {
|
|
zaehler_left(j) = fParam - _vKnotVector(iIndex + 1 - j);
|
|
zaehler_right(j) = _vKnotVector(iIndex + j) - fParam;
|
|
double saved = 0.0;
|
|
for (int r = 0; r < j; r++) {
|
|
double tmp = vFuncVals(r) / (zaehler_right(r + 1) + zaehler_left(j - r));
|
|
vFuncVals(r) = saved + zaehler_right(r + 1) * tmp;
|
|
saved = zaehler_left(j - r) * tmp;
|
|
}
|
|
|
|
vFuncVals(j) = saved;
|
|
}
|
|
}
|
|
|
|
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;
|
|
int p = _iOrder - 1;
|
|
double saved;
|
|
TColStd_Array1OfReal N(0, p);
|
|
|
|
if ((iIndex == 0 && fParam == _vKnotVector(0))
|
|
|| (iIndex == m - p - 1 && fParam == _vKnotVector(m))) {
|
|
return 1.0;
|
|
}
|
|
|
|
if (fParam < _vKnotVector(iIndex) || fParam >= _vKnotVector(iIndex + p + 1)) {
|
|
return 0.0;
|
|
}
|
|
|
|
for (int j = 0; j <= p; j++) {
|
|
if (fParam >= _vKnotVector(iIndex + j) && fParam < _vKnotVector(iIndex + j + 1)) {
|
|
N(j) = 1.0;
|
|
}
|
|
else {
|
|
N(j) = 0.0;
|
|
}
|
|
}
|
|
|
|
for (int k = 1; k <= p; k++) {
|
|
if (N(0) == 0.0) {
|
|
saved = 0.0;
|
|
}
|
|
else {
|
|
saved = ((fParam - _vKnotVector(iIndex)) * N(0))
|
|
/ (_vKnotVector(iIndex + k) - _vKnotVector(iIndex));
|
|
}
|
|
|
|
for (int j = 0; j < p - k + 1; j++) {
|
|
double Tleft = _vKnotVector(iIndex + j + 1);
|
|
double Tright = _vKnotVector(iIndex + j + k + 1);
|
|
|
|
if (N(j + 1) == 0.0) {
|
|
N(j) = saved;
|
|
saved = 0.0;
|
|
}
|
|
else {
|
|
double tmp = N(j + 1) / (Tright - Tleft);
|
|
N(j) = saved + (Tright - fParam) * tmp;
|
|
saved = (fParam - Tleft) * tmp;
|
|
}
|
|
}
|
|
}
|
|
|
|
return N(0);
|
|
}
|
|
|
|
void BSplineBasis::DerivativesOfBasisFunction(int iIndex,
|
|
int iMaxDer,
|
|
double fParam,
|
|
TColStd_Array1OfReal& Derivat)
|
|
{
|
|
int iMax = iMaxDer;
|
|
if (Derivat.Length() != iMax + 1) {
|
|
Standard_RangeError::Raise("BSplineBasis");
|
|
}
|
|
// kth derivatives (k> degrees) are zero
|
|
if (iMax >= _iOrder) {
|
|
for (int i = _iOrder; i <= iMaxDer; i++) {
|
|
Derivat(i) = 0.0;
|
|
}
|
|
iMax = _iOrder - 1;
|
|
}
|
|
|
|
TColStd_Array1OfReal ND(0, iMax);
|
|
int p = _iOrder - 1;
|
|
math_Matrix N(0, p, 0, p);
|
|
double saved;
|
|
|
|
// if value is outside the interval, then function value and all derivatives equal null
|
|
if (fParam < _vKnotVector(iIndex) || fParam >= _vKnotVector(iIndex + p + 1)) {
|
|
for (int k = 0; k <= iMax; k++) {
|
|
Derivat(k) = 0.0;
|
|
}
|
|
return;
|
|
}
|
|
|
|
// Calculate the basis functions of Order 1
|
|
for (int j = 0; j < _iOrder; j++) {
|
|
if (fParam >= _vKnotVector(iIndex + j) && fParam < _vKnotVector(iIndex + j + 1)) {
|
|
N(j, 0) = 1.0;
|
|
}
|
|
else {
|
|
N(j, 0) = 0.0;
|
|
}
|
|
}
|
|
|
|
// Calculate a triangular table of the function values
|
|
for (int k = 1; k < _iOrder; k++) {
|
|
if (N(0, k - 1) == 0.0) {
|
|
saved = 0.0;
|
|
}
|
|
else {
|
|
saved = ((fParam - _vKnotVector(iIndex)) * N(0, k - 1))
|
|
/ (_vKnotVector(iIndex + k) - _vKnotVector(iIndex));
|
|
}
|
|
for (int j = 0; j < p - k + 1; j++) {
|
|
double Tleft = _vKnotVector(iIndex + j + 1);
|
|
double Tright = _vKnotVector(iIndex + j + k + 1);
|
|
|
|
if (N(j + 1, k - 1) == 0.0) {
|
|
N(j, k) = saved;
|
|
saved = 0.0;
|
|
}
|
|
else {
|
|
double tmp = N(j + 1, k - 1) / (Tright - Tleft);
|
|
N(j, k) = saved + (Tright - fParam) * tmp;
|
|
saved = (fParam - Tleft) * tmp;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Function value
|
|
Derivat(0) = N(0, p);
|
|
// Calculate the derivatives from the triangle table
|
|
for (int k = 1; k <= iMax; k++) {
|
|
for (int j = 0; j <= k; j++) {
|
|
// Load the (p-k)th column
|
|
ND(j) = N(j, p - k);
|
|
}
|
|
|
|
for (int jj = 1; jj <= k; jj++) {
|
|
if (ND(0) == 0.0) {
|
|
saved = 0.0;
|
|
}
|
|
else {
|
|
saved = ND(0) / (_vKnotVector(iIndex + p - k + jj) - _vKnotVector(iIndex));
|
|
}
|
|
|
|
for (int j = 0; j < k - jj + 1; j++) {
|
|
double Tleft = _vKnotVector(iIndex + j + 1);
|
|
double Tright = _vKnotVector(iIndex + j + p - k + jj + 1);
|
|
if (ND(j + 1) == 0.0) {
|
|
ND(j) = (p - k + jj) * saved;
|
|
saved = 0.0;
|
|
}
|
|
else {
|
|
double tmp = ND(j + 1) / (Tright - Tleft);
|
|
ND(j) = (p - k + jj) * (saved - tmp);
|
|
saved = tmp;
|
|
}
|
|
}
|
|
}
|
|
|
|
Derivat(k) = ND(0); // kth derivative
|
|
}
|
|
}
|
|
|
|
double BSplineBasis::DerivativeOfBasisFunction(int iIndex, int iMaxDer, double fParam)
|
|
{
|
|
int iMax = iMaxDer;
|
|
|
|
// Function value (0th derivative)
|
|
if (iMax == 0) {
|
|
return BasisFunction(iIndex, fParam);
|
|
}
|
|
|
|
// The kth derivatives (k>degrees) are null
|
|
if (iMax >= _iOrder) {
|
|
return 0.0;
|
|
}
|
|
|
|
TColStd_Array1OfReal ND(0, iMax);
|
|
int p = _iOrder - 1;
|
|
math_Matrix N(0, p, 0, p);
|
|
double saved;
|
|
|
|
// If value is outside the interval, then function value and derivatives equal null
|
|
if (fParam < _vKnotVector(iIndex) || fParam >= _vKnotVector(iIndex + p + 1)) {
|
|
return 0.0;
|
|
}
|
|
|
|
// Calculate the basis functions of Order 1
|
|
for (int j = 0; j < _iOrder; j++) {
|
|
if (fParam >= _vKnotVector(iIndex + j) && fParam < _vKnotVector(iIndex + j + 1)) {
|
|
N(j, 0) = 1.0;
|
|
}
|
|
else {
|
|
N(j, 0) = 0.0;
|
|
}
|
|
}
|
|
|
|
// Calculate triangular table of function values
|
|
for (int k = 1; k < _iOrder; k++) {
|
|
if (N(0, k - 1) == 0.0) {
|
|
saved = 0.0;
|
|
}
|
|
else {
|
|
saved = ((fParam - _vKnotVector(iIndex)) * N(0, k - 1))
|
|
/ (_vKnotVector(iIndex + k) - _vKnotVector(iIndex));
|
|
}
|
|
|
|
for (int j = 0; j < p - k + 1; j++) {
|
|
double Tleft = _vKnotVector(iIndex + j + 1);
|
|
double Tright = _vKnotVector(iIndex + j + k + 1);
|
|
|
|
if (N(j + 1, k - 1) == 0.0) {
|
|
N(j, k) = saved;
|
|
saved = 0.0;
|
|
}
|
|
else {
|
|
double tmp = N(j + 1, k - 1) / (Tright - Tleft);
|
|
N(j, k) = saved + (Tright - fParam) * tmp;
|
|
saved = (fParam - Tleft) * tmp;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Use the triangular table to calculate the derivatives
|
|
for (int j = 0; j <= iMax; j++) {
|
|
// Loading (p-iMax)th column
|
|
ND(j) = N(j, p - iMax);
|
|
}
|
|
|
|
for (int jj = 1; jj <= iMax; jj++) {
|
|
if (ND(0) == 0.0) {
|
|
saved = 0.0;
|
|
}
|
|
else {
|
|
saved = ND(0) / (_vKnotVector(iIndex + p - iMax + jj) - _vKnotVector(iIndex));
|
|
}
|
|
|
|
for (int j = 0; j < iMax - jj + 1; j++) {
|
|
double Tleft = _vKnotVector(iIndex + j + 1);
|
|
double Tright = _vKnotVector(iIndex + j + p - iMax + jj + 1);
|
|
if (ND(j + 1) == 0.0) {
|
|
ND(j) = (p - iMax + jj) * saved;
|
|
saved = 0.0;
|
|
}
|
|
else {
|
|
double tmp = ND(j + 1) / (Tright - Tleft);
|
|
ND(j) = (p - iMax + jj) * (saved - tmp);
|
|
saved = tmp;
|
|
}
|
|
}
|
|
}
|
|
|
|
return ND(0); // iMax-th derivative
|
|
}
|
|
|
|
double BSplineBasis::GetIntegralOfProductOfBSplines(int iIdx1, int iIdx2, int iOrd1, int iOrd2)
|
|
{
|
|
int iMax = CalcSize(iOrd1, iOrd2);
|
|
double dIntegral = 0.0;
|
|
double fMin, fMax;
|
|
|
|
TColStd_Array1OfReal vRoots(0, iMax), vWeights(0, iMax);
|
|
GenerateRootsAndWeights(vRoots, vWeights);
|
|
|
|
/*Calculate the integral*/
|
|
// Integration area
|
|
int iBegin = 0;
|
|
int iEnd = 0;
|
|
FindIntegrationArea(iIdx1, iIdx2, iBegin, iEnd);
|
|
|
|
for (int j = iBegin; j < iEnd; j++) {
|
|
fMax = _vKnotVector(j + 1);
|
|
fMin = _vKnotVector(j);
|
|
|
|
if (fMax > fMin) {
|
|
for (int i = 0; i <= iMax; i++) {
|
|
double fParam = 0.5 * (vRoots(i) + 1) * (fMax - fMin) + fMin;
|
|
dIntegral += 0.5 * (fMax - fMin) * vWeights(i)
|
|
* DerivativeOfBasisFunction(iIdx1, iOrd1, fParam)
|
|
* DerivativeOfBasisFunction(iIdx2, iOrd2, fParam);
|
|
}
|
|
}
|
|
}
|
|
|
|
return dIntegral;
|
|
}
|
|
|
|
void BSplineBasis::GenerateRootsAndWeights(TColStd_Array1OfReal& vRoots,
|
|
TColStd_Array1OfReal& vWeights)
|
|
{
|
|
int iSize = vRoots.Length();
|
|
|
|
// Zeroing the Legendre-Polynomials and the corresponding weights
|
|
if (iSize == 1) {
|
|
vRoots(0) = 0.0;
|
|
vWeights(0) = 2.0;
|
|
}
|
|
else if (iSize == 2) {
|
|
vRoots(0) = 0.57735;
|
|
vWeights(0) = 1.0;
|
|
vRoots(1) = -vRoots(0);
|
|
vWeights(1) = vWeights(0);
|
|
}
|
|
else if (iSize == 4) {
|
|
vRoots(0) = 0.33998;
|
|
vWeights(0) = 0.65214;
|
|
vRoots(1) = 0.86113;
|
|
vWeights(1) = 0.34785;
|
|
vRoots(2) = -vRoots(0);
|
|
vWeights(2) = vWeights(0);
|
|
vRoots(3) = -vRoots(1);
|
|
vWeights(3) = vWeights(1);
|
|
}
|
|
else if (iSize == 6) {
|
|
vRoots(0) = 0.23861;
|
|
vWeights(0) = 0.46791;
|
|
vRoots(1) = 0.66120;
|
|
vWeights(1) = 0.36076;
|
|
vRoots(2) = 0.93246;
|
|
vWeights(2) = 0.17132;
|
|
vRoots(3) = -vRoots(0);
|
|
vWeights(3) = vWeights(0);
|
|
vRoots(4) = -vRoots(1);
|
|
vWeights(4) = vWeights(1);
|
|
vRoots(5) = -vRoots(2);
|
|
vWeights(5) = vWeights(2);
|
|
}
|
|
else if (iSize == 8) {
|
|
vRoots(0) = 0.18343;
|
|
vWeights(0) = 0.36268;
|
|
vRoots(1) = 0.52553;
|
|
vWeights(1) = 0.31370;
|
|
vRoots(2) = 0.79666;
|
|
vWeights(2) = 0.22238;
|
|
vRoots(3) = 0.96028;
|
|
vWeights(3) = 0.10122;
|
|
vRoots(4) = -vRoots(0);
|
|
vWeights(4) = vWeights(0);
|
|
vRoots(5) = -vRoots(1);
|
|
vWeights(5) = vWeights(1);
|
|
vRoots(6) = -vRoots(2);
|
|
vWeights(6) = vWeights(2);
|
|
vRoots(7) = -vRoots(3);
|
|
vWeights(7) = vWeights(3);
|
|
}
|
|
else if (iSize == 10) {
|
|
vRoots(0) = 0.14887;
|
|
vWeights(0) = 0.29552;
|
|
vRoots(1) = 0.43339;
|
|
vWeights(1) = 0.26926;
|
|
vRoots(2) = 0.67940;
|
|
vWeights(2) = 0.21908;
|
|
vRoots(3) = 0.86506;
|
|
vWeights(3) = 0.14945;
|
|
vRoots(4) = 0.97390;
|
|
vWeights(4) = 0.06667;
|
|
vRoots(5) = -vRoots(0);
|
|
vWeights(5) = vWeights(0);
|
|
vRoots(6) = -vRoots(1);
|
|
vWeights(6) = vWeights(1);
|
|
vRoots(7) = -vRoots(2);
|
|
vWeights(7) = vWeights(2);
|
|
vRoots(8) = -vRoots(3);
|
|
vWeights(8) = vWeights(3);
|
|
vRoots(9) = -vRoots(4);
|
|
vWeights(9) = vWeights(4);
|
|
}
|
|
else {
|
|
vRoots(0) = 0.12523;
|
|
vWeights(0) = 0.24914;
|
|
vRoots(1) = 0.36783;
|
|
vWeights(1) = 0.23349;
|
|
vRoots(2) = 0.58731;
|
|
vWeights(2) = 0.20316;
|
|
vRoots(3) = 0.76990;
|
|
vWeights(3) = 0.16007;
|
|
vRoots(4) = 0.90411;
|
|
vWeights(4) = 0.10693;
|
|
vRoots(5) = 0.98156;
|
|
vWeights(5) = 0.04717;
|
|
vRoots(6) = -vRoots(0);
|
|
vWeights(6) = vWeights(0);
|
|
vRoots(7) = -vRoots(1);
|
|
vWeights(7) = vWeights(1);
|
|
vRoots(8) = -vRoots(2);
|
|
vWeights(8) = vWeights(2);
|
|
vRoots(9) = -vRoots(3);
|
|
vWeights(9) = vWeights(3);
|
|
vRoots(10) = -vRoots(4);
|
|
vWeights(10) = vWeights(4);
|
|
vRoots(11) = -vRoots(5);
|
|
vWeights(11) = vWeights(5);
|
|
}
|
|
}
|
|
|
|
void BSplineBasis::FindIntegrationArea(int iIdx1, int iIdx2, int& iBegin, int& iEnd)
|
|
{
|
|
// order by index
|
|
if (iIdx2 < iIdx1) {
|
|
int tmp = iIdx1;
|
|
iIdx1 = iIdx2;
|
|
iIdx2 = tmp;
|
|
}
|
|
|
|
iBegin = iIdx2;
|
|
iEnd = iIdx1 + _iOrder;
|
|
if (iEnd == _vKnotVector.Upper()) {
|
|
iEnd -= 1;
|
|
}
|
|
}
|
|
|
|
int BSplineBasis::CalcSize(int r, int s)
|
|
{
|
|
int iMaxDegree = 2 * (_iOrder - 1) - r - s;
|
|
|
|
if (iMaxDegree < 0) {
|
|
return 0;
|
|
}
|
|
else if (iMaxDegree < 4) {
|
|
return 1;
|
|
}
|
|
else if (iMaxDegree < 8) {
|
|
return 3;
|
|
}
|
|
else if (iMaxDegree < 12) {
|
|
return 5;
|
|
}
|
|
else if (iMaxDegree < 16) {
|
|
return 7;
|
|
}
|
|
else if (iMaxDegree < 20) {
|
|
return 9;
|
|
}
|
|
else {
|
|
return 11;
|
|
}
|
|
}
|
|
|
|
/////////////////// ParameterCorrection
|
|
|
|
ParameterCorrection::ParameterCorrection(unsigned usUOrder,
|
|
unsigned usVOrder,
|
|
unsigned usUCtrlpoints,
|
|
unsigned usVCtrlpoints)
|
|
: _usUOrder(usUOrder)
|
|
, _usVOrder(usVOrder)
|
|
, _usUCtrlpoints(usUCtrlpoints)
|
|
, _usVCtrlpoints(usVCtrlpoints)
|
|
, _vCtrlPntsOfSurf(0, usUCtrlpoints - 1, 0, usVCtrlpoints - 1)
|
|
, _vUKnots(0, usUCtrlpoints - usUOrder + 1)
|
|
, _vVKnots(0, usVCtrlpoints - usVOrder + 1)
|
|
, _vUMults(0, usUCtrlpoints - usUOrder + 1)
|
|
, _vVMults(0, usVCtrlpoints - usVOrder + 1)
|
|
{
|
|
_bGetUVDir = false;
|
|
_bSmoothing = false;
|
|
_fSmoothInfluence = 0.0;
|
|
}
|
|
|
|
void ParameterCorrection::CalcEigenvectors()
|
|
{
|
|
MeshCore::PlaneFit planeFit;
|
|
for (int i = _pvcPoints->Lower(); i <= _pvcPoints->Upper(); i++) {
|
|
const gp_Pnt& pnt = (*_pvcPoints)(i);
|
|
planeFit.AddPoint(Base::Vector3f((float)pnt.X(), (float)pnt.Y(), (float)pnt.Z()));
|
|
}
|
|
|
|
planeFit.Fit();
|
|
_clU = Base::toVector<double>(planeFit.GetDirU());
|
|
_clV = Base::toVector<double>(planeFit.GetDirV());
|
|
_clW = Base::toVector<double>(planeFit.GetNormal());
|
|
}
|
|
|
|
bool ParameterCorrection::DoInitialParameterCorrection(double fSizeFactor)
|
|
{
|
|
// if directions are not given, calculate yourself
|
|
if (!_bGetUVDir) {
|
|
CalcEigenvectors();
|
|
}
|
|
if (!GetUVParameters(fSizeFactor)) {
|
|
return false;
|
|
}
|
|
if (_bSmoothing) {
|
|
if (!SolveWithSmoothing(_fSmoothInfluence)) {
|
|
return false;
|
|
}
|
|
}
|
|
else {
|
|
if (!SolveWithoutSmoothing()) {
|
|
return false;
|
|
}
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
bool ParameterCorrection::GetUVParameters(double fSizeFactor)
|
|
{
|
|
// Eigenvectors as a new base
|
|
Base::Vector3d e[3];
|
|
e[0] = _clU;
|
|
e[1] = _clV;
|
|
e[2] = _clW;
|
|
|
|
// Canonical base of R^3
|
|
Base::Vector3d b[3];
|
|
b[0] = Base::Vector3d(1.0, 0.0, 0.0);
|
|
b[1] = Base::Vector3d(0.0, 1.0, 0.0);
|
|
b[2] = Base::Vector3d(0.0, 0.0, 1.0);
|
|
// Create a right system from the orthogonal eigenvectors
|
|
if ((e[0] % e[1]) * e[2] < 0) {
|
|
Base::Vector3d tmp = e[0];
|
|
e[0] = e[1];
|
|
e[1] = tmp;
|
|
}
|
|
|
|
// Now generate the transpon. Rotation matrix
|
|
Wm4::Matrix3d clRotMatTrans;
|
|
for (int i = 0; i < 3; i++) {
|
|
for (int j = 0; j < 3; j++) {
|
|
clRotMatTrans[i][j] = b[j] * e[i];
|
|
}
|
|
}
|
|
|
|
std::vector<Base::Vector2d> vcProjPts;
|
|
Base::BoundBox2d clBBox;
|
|
|
|
// Calculate the coordinates of the transf. Points and project
|
|
// these on to the x,y-plane of the new coordinate system
|
|
for (int ii = _pvcPoints->Lower(); ii <= _pvcPoints->Upper(); ii++) {
|
|
const gp_Pnt& pnt = (*_pvcPoints)(ii);
|
|
Wm4::Vector3d clProjPnt = clRotMatTrans * Wm4::Vector3d(pnt.X(), pnt.Y(), pnt.Z());
|
|
vcProjPts.emplace_back(clProjPnt.X(), clProjPnt.Y());
|
|
clBBox.Add(Base::Vector2d(clProjPnt.X(), clProjPnt.Y()));
|
|
}
|
|
|
|
if ((clBBox.MaxX == clBBox.MinX) || (clBBox.MaxY == clBBox.MinY)) {
|
|
return false;
|
|
}
|
|
double tx = fSizeFactor * clBBox.MinX - (fSizeFactor - 1.0) * clBBox.MaxX;
|
|
double ty = fSizeFactor * clBBox.MinY - (fSizeFactor - 1.0) * clBBox.MaxY;
|
|
double fDeltaX = (2 * fSizeFactor - 1.0) * (clBBox.MaxX - clBBox.MinX);
|
|
double fDeltaY = (2 * fSizeFactor - 1.0) * (clBBox.MaxY - clBBox.MinY);
|
|
|
|
// Calculate the u,v parameters with u,v from [0,1]
|
|
_pvcUVParam->Init(gp_Pnt2d(0.0, 0.0));
|
|
int ii = 0;
|
|
if (clBBox.MaxX - clBBox.MinX >= clBBox.MaxY - clBBox.MinY) {
|
|
for (const auto& pt : vcProjPts) {
|
|
(*_pvcUVParam)(ii) = gp_Pnt2d((pt.x - tx) / fDeltaX, (pt.y - ty) / fDeltaY);
|
|
ii++;
|
|
}
|
|
}
|
|
else {
|
|
for (const auto& pt : vcProjPts) {
|
|
(*_pvcUVParam)(ii) = gp_Pnt2d((pt.y - ty) / fDeltaY, (pt.x - tx) / fDeltaX);
|
|
ii++;
|
|
}
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
void ParameterCorrection::SetUV(const Base::Vector3d& clU, const Base::Vector3d& clV, bool bUseDir)
|
|
{
|
|
_bGetUVDir = bUseDir;
|
|
if (_bGetUVDir) {
|
|
_clU = clU;
|
|
_clW = clU % clV;
|
|
_clV = _clW % _clU;
|
|
}
|
|
}
|
|
|
|
void ParameterCorrection::GetUVW(Base::Vector3d& clU,
|
|
Base::Vector3d& clV,
|
|
Base::Vector3d& clW) const
|
|
{
|
|
clU = _clU;
|
|
clV = _clV;
|
|
clW = _clW;
|
|
}
|
|
|
|
Base::Vector3d ParameterCorrection::GetGravityPoint() const
|
|
{
|
|
Standard_Integer ulSize = _pvcPoints->Length();
|
|
double x = 0.0, y = 0.0, z = 0.0;
|
|
for (int i = _pvcPoints->Lower(); i <= _pvcPoints->Upper(); i++) {
|
|
const gp_Pnt& pnt = (*_pvcPoints)(i);
|
|
x += pnt.X();
|
|
y += pnt.Y();
|
|
z += pnt.Z();
|
|
}
|
|
|
|
return Base::Vector3d(x / ulSize, y / ulSize, z / ulSize);
|
|
}
|
|
|
|
void ParameterCorrection::ProjectControlPointsOnPlane()
|
|
{
|
|
Base::Vector3d base = GetGravityPoint();
|
|
for (unsigned j = 0; j < _usUCtrlpoints; j++) {
|
|
for (unsigned k = 0; k < _usVCtrlpoints; k++) {
|
|
gp_Pnt pole = _vCtrlPntsOfSurf(j, k);
|
|
Base::Vector3d pnt(pole.X(), pole.Y(), pole.Z());
|
|
pnt.ProjectToPlane(base, _clW);
|
|
pole.SetX(pnt.x);
|
|
pole.SetY(pnt.y);
|
|
pole.SetZ(pnt.z);
|
|
_vCtrlPntsOfSurf(j, k) = pole;
|
|
}
|
|
}
|
|
}
|
|
|
|
Handle(Geom_BSplineSurface) ParameterCorrection::CreateSurface(const TColgp_Array1OfPnt& points,
|
|
int iIter,
|
|
bool bParaCor,
|
|
double fSizeFactor)
|
|
{
|
|
if (_pvcPoints) {
|
|
delete _pvcPoints;
|
|
_pvcPoints = nullptr;
|
|
delete _pvcUVParam;
|
|
_pvcUVParam = nullptr;
|
|
}
|
|
|
|
_pvcPoints = new TColgp_Array1OfPnt(points.Lower(), points.Upper());
|
|
*_pvcPoints = points;
|
|
_pvcUVParam = new TColgp_Array1OfPnt2d(points.Lower(), points.Upper());
|
|
|
|
if (_usUCtrlpoints * _usVCtrlpoints > static_cast<unsigned>(_pvcPoints->Length())) {
|
|
return nullptr; // LGS under-determined
|
|
}
|
|
if (!DoInitialParameterCorrection(fSizeFactor)) {
|
|
return nullptr;
|
|
}
|
|
|
|
// Generate the approximation plane as a B-spline area
|
|
if (iIter < 0) {
|
|
bParaCor = false;
|
|
ProjectControlPointsOnPlane();
|
|
}
|
|
// No further parameter corrections
|
|
else if (iIter == 0) {
|
|
bParaCor = false;
|
|
}
|
|
|
|
if (bParaCor) {
|
|
DoParameterCorrection(iIter);
|
|
}
|
|
|
|
return new Geom_BSplineSurface(_vCtrlPntsOfSurf,
|
|
_vUKnots,
|
|
_vVKnots,
|
|
_vUMults,
|
|
_vVMults,
|
|
_usUOrder - 1,
|
|
_usVOrder - 1);
|
|
}
|
|
|
|
void ParameterCorrection::EnableSmoothing(bool bSmooth, double fSmoothInfl)
|
|
{
|
|
_bSmoothing = bSmooth;
|
|
_fSmoothInfluence = fSmoothInfl;
|
|
}
|
|
|
|
/////////////////// BSplineParameterCorrection
|
|
|
|
|
|
BSplineParameterCorrection::BSplineParameterCorrection(unsigned usUOrder,
|
|
unsigned usVOrder,
|
|
unsigned usUCtrlpoints,
|
|
unsigned usVCtrlpoints)
|
|
: ParameterCorrection(usUOrder, usVOrder, usUCtrlpoints, usVCtrlpoints)
|
|
, _clUSpline(usUCtrlpoints + usUOrder)
|
|
, _clVSpline(usVCtrlpoints + usVOrder)
|
|
, _clSmoothMatrix(0, usUCtrlpoints * usVCtrlpoints - 1, 0, usUCtrlpoints * usVCtrlpoints - 1)
|
|
, _clFirstMatrix(0, usUCtrlpoints * usVCtrlpoints - 1, 0, usUCtrlpoints * usVCtrlpoints - 1)
|
|
, _clSecondMatrix(0, usUCtrlpoints * usVCtrlpoints - 1, 0, usUCtrlpoints * usVCtrlpoints - 1)
|
|
, _clThirdMatrix(0, usUCtrlpoints * usVCtrlpoints - 1, 0, usUCtrlpoints * usVCtrlpoints - 1)
|
|
{
|
|
Init();
|
|
}
|
|
|
|
void BSplineParameterCorrection::Init()
|
|
{
|
|
// Initializations
|
|
_pvcUVParam = nullptr;
|
|
_pvcPoints = nullptr;
|
|
_clFirstMatrix.Init(0.0);
|
|
_clSecondMatrix.Init(0.0);
|
|
_clThirdMatrix.Init(0.0);
|
|
_clSmoothMatrix.Init(0.0);
|
|
|
|
/* Calculate the knot vectors */
|
|
unsigned usUMax = _usUCtrlpoints - _usUOrder + 1;
|
|
unsigned usVMax = _usVCtrlpoints - _usVOrder + 1;
|
|
|
|
// Knot vector for the CAS.CADE class
|
|
// u-direction
|
|
for (unsigned i = 0; i <= usUMax; i++) {
|
|
_vUKnots(i) = static_cast<double>(i) / static_cast<double>(usUMax);
|
|
_vUMults(i) = 1;
|
|
}
|
|
|
|
_vUMults(0) = _usUOrder;
|
|
_vUMults(usUMax) = _usUOrder;
|
|
|
|
// v-direction
|
|
for (unsigned i = 0; i <= usVMax; i++) {
|
|
_vVKnots(i) = static_cast<double>(i) / static_cast<double>(usVMax);
|
|
_vVMults(i) = 1;
|
|
}
|
|
|
|
_vVMults(0) = _usVOrder;
|
|
_vVMults(usVMax) = _usVOrder;
|
|
|
|
// Set the B-spline basic functions
|
|
_clUSpline.SetKnots(_vUKnots, _vUMults, _usUOrder);
|
|
_clVSpline.SetKnots(_vVKnots, _vVMults, _usVOrder);
|
|
}
|
|
|
|
void BSplineParameterCorrection::SetUKnots(const std::vector<double>& afKnots)
|
|
{
|
|
std::size_t numPoints = static_cast<std::size_t>(_usUCtrlpoints);
|
|
std::size_t order = static_cast<std::size_t>(_usUOrder);
|
|
if (afKnots.size() != (numPoints + order)) {
|
|
return;
|
|
}
|
|
|
|
unsigned usUMax = _usUCtrlpoints - _usUOrder + 1;
|
|
|
|
// Knot vector for the CAS.CADE class
|
|
// u-direction
|
|
for (unsigned i = 1; i < usUMax; i++) {
|
|
_vUKnots(i) = afKnots[_usUOrder + i - 1];
|
|
_vUMults(i) = 1;
|
|
}
|
|
|
|
// Set the B-spline basic functions
|
|
_clUSpline.SetKnots(_vUKnots, _vUMults, _usUOrder);
|
|
}
|
|
|
|
void BSplineParameterCorrection::SetVKnots(const std::vector<double>& afKnots)
|
|
{
|
|
std::size_t numPoints = static_cast<std::size_t>(_usVCtrlpoints);
|
|
std::size_t order = static_cast<std::size_t>(_usVOrder);
|
|
if (afKnots.size() != (numPoints + order)) {
|
|
return;
|
|
}
|
|
|
|
unsigned usVMax = _usVCtrlpoints - _usVOrder + 1;
|
|
|
|
// Knot vector for the CAS.CADE class
|
|
// v-direction
|
|
for (unsigned i = 1; i < usVMax; i++) {
|
|
_vVKnots(i) = afKnots[_usVOrder + i - 1];
|
|
_vVMults(i) = 1;
|
|
}
|
|
|
|
// Set the B-spline basic functions
|
|
_clVSpline.SetKnots(_vVKnots, _vVMults, _usVOrder);
|
|
}
|
|
|
|
void BSplineParameterCorrection::DoParameterCorrection(int iIter)
|
|
{
|
|
int i = 0;
|
|
double fMaxDiff = 0.0, fMaxScalar = 1.0;
|
|
double fWeight = _fSmoothInfluence;
|
|
|
|
Base::SequencerLauncher seq("Calc surface...",
|
|
static_cast<size_t>(iIter)
|
|
* static_cast<size_t>(_pvcPoints->Length()));
|
|
|
|
do {
|
|
fMaxScalar = 1.0;
|
|
fMaxDiff = 0.0;
|
|
|
|
Handle(Geom_BSplineSurface) pclBSplineSurf = new Geom_BSplineSurface(_vCtrlPntsOfSurf,
|
|
_vUKnots,
|
|
_vVKnots,
|
|
_vUMults,
|
|
_vVMults,
|
|
_usUOrder - 1,
|
|
_usVOrder - 1);
|
|
|
|
for (int ii = _pvcPoints->Lower(); ii <= _pvcPoints->Upper(); ii++) {
|
|
double fDeltaU, fDeltaV, fU, fV;
|
|
const gp_Pnt& pnt = (*_pvcPoints)(ii);
|
|
gp_Vec P(pnt.X(), pnt.Y(), pnt.Z());
|
|
gp_Pnt PntX;
|
|
gp_Vec Xu, Xv, Xuv, Xuu, Xvv;
|
|
// Calculate the first two derivatives and point at (u,v)
|
|
gp_Pnt2d& uvValue = (*_pvcUVParam)(ii);
|
|
pclBSplineSurf->D2(uvValue.X(), uvValue.Y(), PntX, Xu, Xv, Xuu, Xvv, Xuv);
|
|
gp_Vec X(PntX.X(), PntX.Y(), PntX.Z());
|
|
gp_Vec ErrorVec = X - P;
|
|
|
|
// Calculate Xu x Xv the normal in X(u,v)
|
|
gp_Dir clNormal = Xu ^ Xv;
|
|
|
|
// Check, if X = P
|
|
if (!(X.IsEqual(P, 0.001, 0.001))) {
|
|
ErrorVec.Normalize();
|
|
if (fabs(clNormal * ErrorVec) < fMaxScalar) {
|
|
fMaxScalar = fabs(clNormal * ErrorVec);
|
|
}
|
|
}
|
|
|
|
fDeltaU = ((P - X) * Xu) / ((P - X) * Xuu - Xu * Xu);
|
|
if (fabs(fDeltaU) < Precision::Confusion()) {
|
|
fDeltaU = 0.0;
|
|
}
|
|
fDeltaV = ((P - X) * Xv) / ((P - X) * Xvv - Xv * Xv);
|
|
if (fabs(fDeltaV) < Precision::Confusion()) {
|
|
fDeltaV = 0.0;
|
|
}
|
|
|
|
// Replace old u/v values with new ones
|
|
fU = uvValue.X() - fDeltaU;
|
|
fV = uvValue.Y() - fDeltaV;
|
|
if (fU <= 1.0 && fU >= 0.0 && fV <= 1.0 && fV >= 0.0) {
|
|
uvValue.SetX(fU);
|
|
uvValue.SetY(fV);
|
|
fMaxDiff = std::max<double>(fabs(fDeltaU), fMaxDiff);
|
|
fMaxDiff = std::max<double>(fabs(fDeltaV), fMaxDiff);
|
|
}
|
|
|
|
seq.next();
|
|
}
|
|
|
|
if (_bSmoothing) {
|
|
fWeight *= 0.5f;
|
|
SolveWithSmoothing(fWeight);
|
|
}
|
|
else {
|
|
SolveWithoutSmoothing();
|
|
}
|
|
|
|
i++;
|
|
} while (i < iIter && fMaxDiff > Precision::Confusion() && fMaxScalar < 0.99);
|
|
}
|
|
|
|
bool BSplineParameterCorrection::SolveWithoutSmoothing()
|
|
{
|
|
unsigned ulSize = _pvcPoints->Length();
|
|
unsigned ulDim = _usUCtrlpoints * _usVCtrlpoints;
|
|
math_Matrix M(0, ulSize - 1, 0, ulDim - 1);
|
|
math_Matrix Xx(0, ulDim - 1, 0, 0);
|
|
math_Matrix Xy(0, ulDim - 1, 0, 0);
|
|
math_Matrix Xz(0, ulDim - 1, 0, 0);
|
|
math_Vector bx(0, ulSize - 1);
|
|
math_Vector by(0, ulSize - 1);
|
|
math_Vector bz(0, ulSize - 1);
|
|
|
|
// Determining the coefficient matrix of the overdetermined LGS
|
|
for (unsigned i = 0; i < ulSize; i++) {
|
|
const gp_Pnt2d& uvValue = (*_pvcUVParam)(i);
|
|
double fU = uvValue.X();
|
|
double fV = uvValue.Y();
|
|
unsigned ulIdx = 0;
|
|
|
|
// Pre-calculation of the values of the base functions
|
|
std::vector<double> basisU(_usUCtrlpoints);
|
|
for (unsigned j = 0; j < _usUCtrlpoints; j++) {
|
|
basisU[j] = _clUSpline.BasisFunction(j, fU);
|
|
}
|
|
std::vector<double> basisV(_usVCtrlpoints);
|
|
for (unsigned k = 0; k < _usVCtrlpoints; k++) {
|
|
basisV[k] = _clVSpline.BasisFunction(k, fV);
|
|
}
|
|
|
|
for (unsigned j = 0; j < _usUCtrlpoints; j++) {
|
|
double valueU = basisU[j];
|
|
if (valueU == 0.0) {
|
|
for (unsigned k = 0; k < _usVCtrlpoints; k++) {
|
|
M(i, ulIdx) = 0.0;
|
|
ulIdx++;
|
|
}
|
|
}
|
|
else {
|
|
for (unsigned k = 0; k < _usVCtrlpoints; k++) {
|
|
M(i, ulIdx) = valueU * basisV[k];
|
|
ulIdx++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Determine the right side
|
|
for (int ii = _pvcPoints->Lower(); ii <= _pvcPoints->Upper(); ii++) {
|
|
const gp_Pnt& pnt = (*_pvcPoints)(ii);
|
|
bx(ii) = pnt.X();
|
|
by(ii) = pnt.Y();
|
|
bz(ii) = pnt.Z();
|
|
}
|
|
|
|
// Solve the over-determined LGS with Householder-Transformation
|
|
math_Householder hhX(M, bx);
|
|
math_Householder hhY(M, by);
|
|
math_Householder hhZ(M, bz);
|
|
|
|
if (!(hhX.IsDone() && hhY.IsDone() && hhZ.IsDone())) {
|
|
// LGS could not be solved
|
|
return false;
|
|
}
|
|
Xx = hhX.AllValues();
|
|
Xy = hhY.AllValues();
|
|
Xz = hhZ.AllValues();
|
|
|
|
unsigned ulIdx = 0;
|
|
for (unsigned j = 0; j < _usUCtrlpoints; j++) {
|
|
for (unsigned k = 0; k < _usVCtrlpoints; k++) {
|
|
_vCtrlPntsOfSurf(j, k) = gp_Pnt(Xx(ulIdx, 0), Xy(ulIdx, 0), Xz(ulIdx, 0));
|
|
ulIdx++;
|
|
}
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
namespace Reen
|
|
{
|
|
class ScalarProduct
|
|
{
|
|
public:
|
|
explicit 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;
|
|
};
|
|
} // namespace Reen
|
|
|
|
bool BSplineParameterCorrection::SolveWithSmoothing(double fWeight)
|
|
{
|
|
unsigned ulSize = _pvcPoints->Length();
|
|
unsigned ulDim = _usUCtrlpoints * _usVCtrlpoints;
|
|
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, ulDim - 1);
|
|
math_Vector Mby(0, ulDim - 1);
|
|
math_Vector Mbz(0, ulDim - 1);
|
|
|
|
// Determining the coefficient matrix of the overdetermined LGS
|
|
for (unsigned i = 0; i < ulSize; i++) {
|
|
const gp_Pnt2d& uvValue = (*_pvcUVParam)(i);
|
|
double fU = uvValue.X();
|
|
double fV = uvValue.Y();
|
|
int ulIdx = 0;
|
|
|
|
// Pre-calculation of the values of the basis functions
|
|
std::vector<double> basisU(_usUCtrlpoints);
|
|
for (unsigned j = 0; j < _usUCtrlpoints; j++) {
|
|
basisU[j] = _clUSpline.BasisFunction(j, fU);
|
|
}
|
|
std::vector<double> basisV(_usVCtrlpoints);
|
|
for (unsigned k = 0; k < _usVCtrlpoints; k++) {
|
|
basisV[k] = _clVSpline.BasisFunction(k, fV);
|
|
}
|
|
|
|
for (unsigned j = 0; j < _usUCtrlpoints; j++) {
|
|
double valueU = basisU[j];
|
|
if (valueU == 0.0) {
|
|
for (unsigned k = 0; k < _usVCtrlpoints; k++) {
|
|
M(i, ulIdx) = 0.0;
|
|
ulIdx++;
|
|
}
|
|
}
|
|
else {
|
|
for (unsigned k = 0; k < _usVCtrlpoints; k++) {
|
|
M(i, ulIdx) = valueU * basisV[k];
|
|
ulIdx++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// The product of its transform and itself results in the quadratic
|
|
// system matrix (slowly)
|
|
#if 0
|
|
math_Matrix MTM = M.TMultiply(M);
|
|
#elif 0
|
|
math_Matrix MTM(0, ulDim - 1, 0, ulDim - 1);
|
|
for (unsigned m = 0; m < ulDim; m++) {
|
|
math_Vector Mm = M.Col(m);
|
|
for (unsigned n = m; n < ulDim; 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);
|
|
// NOLINTBEGIN
|
|
QFuture<std::vector<double>> future =
|
|
QtConcurrent::mapped(columns, std::bind(&ScalarProduct::multiply, &scalar, sp::_1));
|
|
// NOLINTEND
|
|
QFutureWatcher<std::vector<double>> watcher;
|
|
watcher.setFuture(future);
|
|
watcher.waitForFinished();
|
|
|
|
math_Matrix MTM(0, ulDim - 1, 0, ulDim - 1);
|
|
int rowIndex = 0;
|
|
for (const auto& it : future) {
|
|
int colIndex = 0;
|
|
for (std::vector<double>::const_iterator jt = it.begin(); jt != it.end();
|
|
++jt, colIndex++) {
|
|
MTM(rowIndex, colIndex) = *jt;
|
|
}
|
|
rowIndex++;
|
|
}
|
|
#endif
|
|
|
|
// Determine the right side
|
|
for (int ii = _pvcPoints->Lower(); ii <= _pvcPoints->Upper(); ii++) {
|
|
const gp_Pnt& pnt = (*_pvcPoints)(ii);
|
|
bx(ii) = pnt.X();
|
|
by(ii) = pnt.Y();
|
|
bz(ii) = pnt.Z();
|
|
}
|
|
for (unsigned i = 0; i < ulDim; i++) {
|
|
math_Vector Mi = M.Col(i);
|
|
Mbx(i) = Mi * bx;
|
|
Mby(i) = Mi * by;
|
|
Mbz(i) = Mi * bz;
|
|
}
|
|
|
|
// Solve the LGS with the LU decomposition
|
|
math_Gauss mgGaussX(MTM + fWeight * _clSmoothMatrix);
|
|
math_Gauss mgGaussY(MTM + fWeight * _clSmoothMatrix);
|
|
math_Gauss mgGaussZ(MTM + fWeight * _clSmoothMatrix);
|
|
|
|
mgGaussX.Solve(Mbx, Xx);
|
|
if (!mgGaussX.IsDone()) {
|
|
return false;
|
|
}
|
|
|
|
mgGaussY.Solve(Mby, Xy);
|
|
if (!mgGaussY.IsDone()) {
|
|
return false;
|
|
}
|
|
|
|
mgGaussZ.Solve(Mbz, Xz);
|
|
if (!mgGaussZ.IsDone()) {
|
|
return false;
|
|
}
|
|
|
|
unsigned ulIdx = 0;
|
|
for (unsigned j = 0; j < _usUCtrlpoints; j++) {
|
|
for (unsigned k = 0; k < _usVCtrlpoints; k++) {
|
|
_vCtrlPntsOfSurf(j, k) = gp_Pnt(Xx(ulIdx), Xy(ulIdx), Xz(ulIdx));
|
|
ulIdx++;
|
|
}
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
void BSplineParameterCorrection::CalcSmoothingTerms(bool bRecalc,
|
|
double fFirst,
|
|
double fSecond,
|
|
double fThird)
|
|
{
|
|
if (bRecalc) {
|
|
Base::SequencerLauncher seq("Initializing...",
|
|
static_cast<size_t>(3) * static_cast<size_t>(_usUCtrlpoints)
|
|
* static_cast<size_t>(_usUCtrlpoints)
|
|
* static_cast<size_t>(_usVCtrlpoints)
|
|
* static_cast<size_t>(_usVCtrlpoints));
|
|
CalcFirstSmoothMatrix(seq);
|
|
CalcSecondSmoothMatrix(seq);
|
|
CalcThirdSmoothMatrix(seq);
|
|
}
|
|
|
|
_clSmoothMatrix = fFirst * _clFirstMatrix + fSecond * _clSecondMatrix + fThird * _clThirdMatrix;
|
|
}
|
|
|
|
void BSplineParameterCorrection::CalcFirstSmoothMatrix(Base::SequencerLauncher& seq)
|
|
{
|
|
unsigned m = 0;
|
|
for (unsigned k = 0; k < _usUCtrlpoints; k++) {
|
|
for (unsigned l = 0; l < _usVCtrlpoints; l++) {
|
|
unsigned n = 0;
|
|
|
|
for (unsigned i = 0; i < _usUCtrlpoints; i++) {
|
|
for (unsigned j = 0; j < _usVCtrlpoints; j++) {
|
|
_clFirstMatrix(m, n) = _clUSpline.GetIntegralOfProductOfBSplines(i, k, 1, 1)
|
|
* _clVSpline.GetIntegralOfProductOfBSplines(j, l, 0, 0)
|
|
+ _clUSpline.GetIntegralOfProductOfBSplines(i, k, 0, 0)
|
|
* _clVSpline.GetIntegralOfProductOfBSplines(j, l, 1, 1);
|
|
seq.next();
|
|
n++;
|
|
}
|
|
}
|
|
m++;
|
|
}
|
|
}
|
|
}
|
|
|
|
void BSplineParameterCorrection::CalcSecondSmoothMatrix(Base::SequencerLauncher& seq)
|
|
{
|
|
unsigned m = 0;
|
|
for (unsigned k = 0; k < _usUCtrlpoints; k++) {
|
|
for (unsigned l = 0; l < _usVCtrlpoints; l++) {
|
|
unsigned n = 0;
|
|
|
|
for (unsigned i = 0; i < _usUCtrlpoints; i++) {
|
|
for (unsigned 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)
|
|
* _clVSpline.GetIntegralOfProductOfBSplines(j, l, 1, 1)
|
|
+ _clUSpline.GetIntegralOfProductOfBSplines(i, k, 0, 0)
|
|
* _clVSpline.GetIntegralOfProductOfBSplines(j, l, 2, 2);
|
|
seq.next();
|
|
n++;
|
|
}
|
|
}
|
|
m++;
|
|
}
|
|
}
|
|
}
|
|
|
|
void BSplineParameterCorrection::CalcThirdSmoothMatrix(Base::SequencerLauncher& seq)
|
|
{
|
|
unsigned m = 0;
|
|
for (unsigned k = 0; k < _usUCtrlpoints; k++) {
|
|
for (unsigned l = 0; l < _usVCtrlpoints; l++) {
|
|
unsigned n = 0;
|
|
|
|
for (unsigned i = 0; i < _usUCtrlpoints; i++) {
|
|
for (unsigned j = 0; j < _usVCtrlpoints; j++) {
|
|
_clThirdMatrix(m, n) = _clUSpline.GetIntegralOfProductOfBSplines(i, k, 3, 3)
|
|
* _clVSpline.GetIntegralOfProductOfBSplines(j, l, 0, 0)
|
|
+ _clUSpline.GetIntegralOfProductOfBSplines(i, k, 3, 1)
|
|
* _clVSpline.GetIntegralOfProductOfBSplines(j, l, 0, 2)
|
|
+ _clUSpline.GetIntegralOfProductOfBSplines(i, k, 1, 3)
|
|
* _clVSpline.GetIntegralOfProductOfBSplines(j, l, 2, 0)
|
|
+ _clUSpline.GetIntegralOfProductOfBSplines(i, k, 1, 1)
|
|
* _clVSpline.GetIntegralOfProductOfBSplines(j, l, 2, 2)
|
|
+ _clUSpline.GetIntegralOfProductOfBSplines(i, k, 2, 2)
|
|
* _clVSpline.GetIntegralOfProductOfBSplines(j, l, 1, 1)
|
|
+ _clUSpline.GetIntegralOfProductOfBSplines(i, k, 0, 2)
|
|
* _clVSpline.GetIntegralOfProductOfBSplines(j, l, 3, 1)
|
|
+ _clUSpline.GetIntegralOfProductOfBSplines(i, k, 2, 0)
|
|
* _clVSpline.GetIntegralOfProductOfBSplines(j, l, 1, 3)
|
|
+ _clUSpline.GetIntegralOfProductOfBSplines(i, k, 0, 0)
|
|
* _clVSpline.GetIntegralOfProductOfBSplines(j, l, 3, 3);
|
|
seq.next();
|
|
n++;
|
|
}
|
|
}
|
|
m++;
|
|
}
|
|
}
|
|
}
|
|
|
|
void BSplineParameterCorrection::EnableSmoothing(bool bSmooth, double fSmoothInfl)
|
|
{
|
|
EnableSmoothing(bSmooth, fSmoothInfl, 1.0, 0.0, 0.0);
|
|
}
|
|
|
|
void BSplineParameterCorrection::EnableSmoothing(bool bSmooth,
|
|
double fSmoothInfl,
|
|
double fFirst,
|
|
double fSec,
|
|
double fThird)
|
|
{
|
|
if (_bSmoothing && bSmooth) {
|
|
CalcSmoothingTerms(false, fFirst, fSec, fThird);
|
|
}
|
|
else if (bSmooth) {
|
|
CalcSmoothingTerms(true, fFirst, fSec, fThird);
|
|
}
|
|
|
|
ParameterCorrection::EnableSmoothing(bSmooth, fSmoothInfl);
|
|
}
|
|
|
|
const math_Matrix& BSplineParameterCorrection::GetFirstSmoothMatrix() const
|
|
{
|
|
return _clFirstMatrix;
|
|
}
|
|
|
|
const math_Matrix& BSplineParameterCorrection::GetSecondSmoothMatrix() const
|
|
{
|
|
return _clSecondMatrix;
|
|
}
|
|
|
|
const math_Matrix& BSplineParameterCorrection::GetThirdSmoothMatrix() const
|
|
{
|
|
return _clThirdMatrix;
|
|
}
|
|
|
|
void BSplineParameterCorrection::SetFirstSmoothMatrix(const math_Matrix& rclMat)
|
|
{
|
|
_clFirstMatrix = rclMat;
|
|
}
|
|
|
|
void BSplineParameterCorrection::SetSecondSmoothMatrix(const math_Matrix& rclMat)
|
|
{
|
|
_clSecondMatrix = rclMat;
|
|
}
|
|
|
|
void BSplineParameterCorrection::SetThirdSmoothMatrix(const math_Matrix& rclMat)
|
|
{
|
|
_clThirdMatrix = rclMat;
|
|
}
|