NURBS求取basis函数的代码

// NurbsBasis.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include <stdio.h>
#include <iostream>
using namespace std;

int FindSpan(int n,int p,double u,double *U)
{
if (u == U[n+1]) return n;
int low = p; int high = n+1;
int mid = (low+high)/2;
while (u < U[mid] || u >= U[mid+1]) {
if (u < U[mid]) high = mid;
else low = mid;
mid = (low+high)/2;
}
return mid;
}

void BasisFuns(int i, double u, int p, double *U, double *N)
{
double saved = 0.0,temp = 0.0;
double* left = new double[ p + 1 ];
double* right = new double[ p + 1 ];
N[0] = 1.0;
for( int j = 1; j <= p; j++ )
{
left[ j ] = u - U[ i + 1- j ];
right[ j ] = U[i + j ] - u;
saved = 0.0;
for( int r = 0; r < j; r++ )
{
// Lower triangle
temp = right[ r + 1 ] + left[ j - r ];
temp = N[r] / temp;

// Upper triangle
N[ r ] = saved + right[ r + 1 ]*temp;
saved = left[ j - r ] * temp;
}
N[ j ] = saved;
}
}

void DersBasisFuns( int i, double u, int p, int n, double* U, double ** & ders )
{
// Compute nonzero basis functions and their
// derivatives. First section is A2.2 modified
// to store functions and knot differences
// Input: i, u, p, n, U
// Output: ders
double saved,temp;
double ** ndu;
double ** a;

// Set up a [p+1] * [p+1] matrix ndu
ndu = new double* [p+1];
if ( ndu == NULL )
{
cout << " error " << endl;
}

for( int j = 0; j <= p; j++ )
{
ndu[ j ] = new double [ p + 1 ];
if( ndu[ j ] == NULL )
{
cout << " error " << endl;
for( int k = 0; k <= j; k++ )
{
delete [ ] ndu[ k ];
}
delete [ ] ndu;
}
}

// Set up a [2] * [p+1] matrix a

a = new double* [2];
if ( a == NULL )
{
cout << " error " << endl;
}

for( int j = 0; j <= 1; j++ )
{
a[ j ] = new double [ p + 1 ];
if( a[ j ] == NULL )
{
cout << " error " << endl;
for( int k = 0; k <= j; k++ )
{
delete [ ] a[ k ];
}
delete [ ] a;
}
}

// Set up a [n+1] * [p+1] matrix ders

/*ders = new double* [ n + 1 ];
if ( ders == NULL )
{
cout << " error " << endl;
}

for( int j = 0; j <= 1; j++ )
{
ders[ j ] = new double [ p + 1 ];
if( ders[ j ] == NULL )
{
cout << " error " << endl;
for( int k = 0; k <= j; k++ )
{
delete [ ] ders[ k ];
}
delete [ ] ders;
}
}*/

ndu[ 0 ][ 0 ] = 1.0;
double* left = new double[ p + 1 ];
double* right = new double[ p + 1 ];

for( int j = 1; j <= p; j++ )
{
left[ j ] = u - U[ i + 1- j ];
right[ j ] = U[i + j ] - u;
saved = 0.0;
for( int r = 0; r < j; r++ )
{
// Lower triangle
ndu[ j ][ r ] = right[ r + 1 ] + left[ j - r ];
temp = ndu[ r ][ j - 1 ] / ndu[ j ][ r ];

// Upper triangle
ndu[ r ][ j ] = saved + right[ r + 1 ]*temp;
saved = left[ j - r ] * temp;
}
ndu[ j ][ j ] = saved;
}

for( int j = 0; j <= p; j++ )
{
// Load the basis functions
ders[ 0 ][ j ] = ndu[ j ][ p ];
}

int s1, s2, rk, pk, j1, j2;
double d;
// This section computes the derivatives (Eq.[2.9])
for( int r = 0; r <= p; r++ )
{
s1 = 0; s2 = 1; // Alternate rows in array a
a[0][0] = 1.0;
// Loop to compute k'th derivative
for( int k = 1; k <= n; k++ )
{
d = 0.0;
rk = r - k;
pk = p - k;
if( r >= k )
{
a[ s2 ][ 0 ] = a[ s1 ][ 0 ] / ndu[ pk + 1][ rk ];
d = a[ s2 ][ 0 ] * ndu[ rk ][ pk ];
}
if( rk >= -1 )
j1 = 1;
else
j1 = - rk;
if( r - 1 <= pk )
j2 = k - 1;
else
j2 = p - r;
for( int j = j1; j <= j2; j++ )
{
a[ s2 ][ j ] = ( a[ s1 ][ j ] - a[ s1 ][ j - 1 ] ) / ndu[ pk + 1 ][ rk + j ];
d += a[ s2 ][ j ] * ndu[ rk + j ][ pk ];
}
if( r<= pk )
{
a[ s2 ][ k ] = - a[ s1 ][ k - 1 ] / ndu[ pk + 1 ][ r ];
d += a[ s2 ][ k ] * ndu[ r ][ pk ];
}
ders[ k ][ r ] = d;
// Switch rows
int j = s1;
s1 = s2;
s2 = j;
}
}
// Multiply through by the correct factors
// (Eq.2.9)
int r = p;
for( int k =1; k<=n; k++)
{
for( int j =0; j<= p; j++ )
{
ders[k][j] *=r;
}
r *= (p-k);
}
}
int _tmain(int argc, _TCHAR* argv[])
{
int i0 = 0;
int p0 = 2;
double U0[] = {0,0,0,1,2,3,4,4,5,5,5};
double u0 = 5.0/2.0;
int num0 = sizeof(U0)/sizeof(double) - p0 - 1;
i0 = FindSpan(num0,p0,u0,U0);
double *N0 = new double[p0+1];
BasisFuns(i0,u0,p0,U0,N0);
int nders = 2;
double ** ders0;

ders0 = new double* [nders+1];
if ( ders0 == NULL )
{
cout << " error " << endl;
}

for( int j = 0; j < nders+1; j++ )
{
ders0[ j ] = new double [ p0 + 1 ];
if( ders0[ j ] == NULL )
{
cout << " error " << endl;
for( int k = 0; k <= j; k++ )
{
delete [ ] ders0[ k ];
}
delete [ ] ders0;
}
}

DersBasisFuns(i0,u0,p0,nders,U0,ders0);
for( int i = 0; i < nders+1; i++ )
{
for (int j = 0; j < p0+1; j++)
{
double derivative = ders0[i][j];
cout<<ders0[i][j]<<"\t";
}
cout<<endl;
}
return 0;
}

NURBS求取basis函数的代码相关推荐

BP神经网络+遗传算法：求取非线性函数极值（一）
1.问题描述对于未知的非线性函数,仅通过函数的输入输出数据难以准确寻找其函数极值,这类问题可以通过利用BP神经网络的非线性拟合能力和遗传算法的非线性寻优能力来解决. 设非线性函数表达式如下: 其函数 ...
Python中ArcPy实现对大量长时间序列栅格遥感影像批量逐像元求取像素平均值
本文介绍基于Python中ArcPy模块,对大量长时间序列栅格遥感影像文件的每一个像元进行多时序平均值的求取. 在遥感应用中,我们经常需要对某一景遥感影像中的全部像元的像素值进行平均值求取-- ...
海伦公式用计算机语言怎么写,python中海伦公式求取三角形面积的示例
python中海伦公式求取三角形面积的示例发布时间:2020-12-07 10:01:44 来源:亿速云阅读:143 作者:小新这篇文章将为大家详细讲解有关python中海伦公式求取三角形面积的 ...
python海伦公式求三角形面积_python编程实战：海伦公式求取三角形的面积
之前小编向大家介绍了在python中求取三角形面积的方法:三角形面积代码.大家对三角形面积的求取有了一定的了解,我们也知道计算机可以进行高精度的计算,那如果说在测量土地的面积的时候,不测三角形的高,只 ...
874. 筛法求欧拉函数
874. 筛法求欧拉函数题目代码题目给定一个正整数 n,求 1∼n 中每个数的欧拉函数之和. 输入格式共一行,包含一个整数 n. 输出格式共一行,包含一个整数,表示 1∼n 中每个数的欧拉 ...
11.9 至 11.17 四道典型题记录： Counter 弹出 | map函数 | 子集求取 | 有序字符桶分装
11.9 至 11.17 四道典型题记录: Counter 弹出 | map函数 | 子集求取 | 有序字符桶分装昨天休息的时候一直在想应该学习哪种语言,我想这也是好多人发愁无法下手的原因之一 ...
计算机数据的平均函数是,excel软件中数据的平均值怎么求取
excel软件中数据的平均值怎么求取腾讯视频/爱奇艺/优酷/外卖充值4折起 Excel是我们经常使用的数据处理工具,我们在编写表格的时候经常会遇到需要求平均值的情况,接下来小编就教大家怎么在exc ...
c# 椭圆拟合库_利用C#版OpenCV实现圆心求取实例代码
前言 OpenCVSharp是OpenCV的.NET wrapper,是一名日本工程师开发的,项目地址为:https://github.com/shimat/opencvsharp. 该源码是 BSD ...
matlab求阈值的函数,小波分析中matlab阈值获取函数及其应用附程序代码.doc
小波分析中matlab阈值获取函数及其应用附程序代码.doc 1.小波分析中MATLAB阈值获取函数MATLAB中实现阈值获取的函数有DDENCMP.THSELECT.WBMPEN和WWDCBM,下面 ...

NURBS求取basis函数的代码

NURBS求取basis函数的代码相关推荐

最新文章

热门文章