三次样条/曲线拟合
Cubic spline / curve fitting
我需要确定光照变化的参数，它由这个连续分段多项式 C(t) 定义，其中 f(t) 是由两个边界点 (t1,c) 定义的三次曲线和 (t2,0)，还有 f'(t1)=0 和 f'(t2)=0。
Original Paper: Texture-Consistent Shadow Removal
强度曲线是从阴影边界的法线采样的，看起来像这样：
每一行都是样本，显示光照 change.So X 是列数，Y 是像素强度。
我有这样的真实数据（一个样本从所有样本中平均得出）：
我总共有N个样本，我需要确定参数(c,t1,t2)
我该怎么做？
我尝试通过在 Matlab 中求解线性方程来实现：
avr_curve为平均曲线，对所有样本取平均值
f(x)= x^3+a2*x^2+a1*x1+a0 是三次函数
%t1,t2 selected by hand
t1= 10;
t2= 15;
offset=10;
avr_curve= [41, 40, 40, 41, 41, 42, 42, 43, 43, 43, 51, 76, 98, 104, 104, 103, 104, 105, 105, 107, 105];
%gradx= convn(avr_curve,[-1 1],'same');
A= zeros(2*offset+1,3);
%b= zeros(2*offset+1,1);
b= avr_curve';
%for i= 1:2*offset+1
for i=t1:t2
i
x= i-offset-1
A(i,1)= x^2; %a2
A(i,2)= x; %a1
A(i,3)= 1; %a0
b(i,1)= b(i,1)-x^3;
end
u= A\b;
figure,plot(avr_curve(t1:t2))
%estimated cubic curve
for i= 1:2*offset+1
x= i-offset-1;
fx(i)=x^3+u(1)*x^2+u(2)*x+u(3);
end
figure,plot(fx(t1:t2))
[t1 t2] avr_curve 的一部分
我得到的三次曲线（不像avr_curve）
所以我做错了什么？
更新：
似乎我的错误是由于我使用 3 个变量对三次多项式进行建模，如下所示：
f(x)= x^3+a2*x^2+a1*x1+a0 - 3 variables
但后来我使用了 4 个变量，一切似乎都正常：
f(x)= a3*x^3+a2*x^2+a1*x1+a0 - 4 variables
这是 Matlab 中的代码：
%defined by hand
t1= 10;
t2= 14;
avr_curve= [41, 40, 40, 41, 41, 42, 42, 43, 43, 43, 51, 76, 98, 104, 104, 103, 104, 105, 105, 107, 105];
x= [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21];
%x= [-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]; %real x axis
%%%model 1
%%f(x)= x^3+a2*x^2+a1*x1+a0 - 3 variables
%A= zeros(4,3);
%b= [43 104]';
%%cubic equation at t1
%A(1,1)= t1^2; %a2
%A(1,2)= t1; %a1
%A(1,3)= 1; %a0
%b(1,1)= b(1,1)-t1^3;
%%cubic equation at t2
%A(2,1)= t2^2; %a2
%A(2,2)= t2; %a1
%A(2,3)= 1; %a0
%b(2,1)= b(2,1)-t1^3;
%%1st derivative at t1
%A(3,1)= 2*t1; %a2
%A(3,2)= 1; %a1
%A(3,3)= 0; %a0
%b(3,1)= -3*t1^2;
%%1st derivative at t2
%A(4,1)= 2*t2; %a2
%A(4,2)= 1; %a1
%A(4,3)= 0; %a0
%b(4,1)= -3*t2^2;
%model 2
%f(x)= a3*x^3+a2*x^2+a1*x1+a0 - 4 variables
A= zeros(4,4);
b= [43 104]';
%cubic equation at t1
A(1,1)= t1^3; %a3
A(1,2)= t1^2; %a2
A(1,3)= t1; %a1
A(1,4)= 1; %a0
b(1,1)= b(1,1);
%cubic equation at t2
A(2,1)= t2^3; %a3
A(2,2)= t2^2; %a2
A(2,3)= t2; %a1
A(2,4)= 1; %a0
b(2,1)= b(2,1);
%1st derivative at t1
A(3,1)= 3*t1^2; %a3
A(3,2)= 2*t1; %a2
A(3,3)= 1; %a1
A(3,4)= 0; %a0
b(3,1)= 0;
%1st derivative at t2
A(4,1)= 3*t2^2; %a3
A(4,2)= 2*t2; %a2
A(4,3)= 1; %a1
A(4,4)= 0; %a0
b(4,1)= 0;
size(A)
size(b)
u= A\b;
u
%estimated cubic curve
%dx=[1:21]; % global view
dx=t1-1:t2+1; % local view in [t1 t2]
for x= dx
%fx(x)=x^3+u(1)*x^2+u(2)*x+u(3); % model 1
fx(x)= u(1)*x^3+u(2)*x^2+u(3)*x+u(4); % model 2
end
err= 0;
for x= dx
err= err+(fx(x)-avr_curve(x))^2;
end
err
figure,plot(dx,avr_curve(dx),dx,fx(dx))
样条区间 [t1-1 t2+1]
并且在整个间隔
免责声明
我无法对下面给出的代码或方法的正确性提供任何保证，在使用任何这些之前请务必使用您的批判性判断力。
0。定义问题
你有这个分段定义的函数
其中f(t)为三次函数，为了唯一标识它，给出如下附加条件
您想恢复参数 t1、t2 和 sigma 的最佳值最小化给定点集的误差。
这本质上是 least squares 意义上的曲线拟合。
1 参数化 f(t) 三次函数
为了计算候选 Cl(t) 函数与点集之间的误差，我们需要计算 f(t)，它的一般形式（立方体）是
看来我们还有四个额外的参数需要考虑。实际上，此参数完全由 free 三个参数 t1、t2 和 定义西格玛.
重要的是不要将 free 参数与相关参数混淆。
给定 f(t) 的附加条件，我们可以建立这个线性系统
给出了一个解决方案（正如预期的那样）
这告诉我们如何在给定三个自由参数的情况下计算立方体的参数。
这样Cl(t)就完全确定了，现在是时候寻找最佳人选了
2 最小化误差
我现在通常会选择最小二乘法。
由于这不是线性函数，因此没有用于计算 sigma、t1 和 t2 的封闭形式.
但是有数值方法，比如 Gauss-Newton 一个。
然而，需要以一种或另一种方式计算关于三个参数的偏导数。
我不知道如何计算像 t1 这样的分离参数的导数。
我搜索了 MathSE，发现 this 个解决相同问题的问题，但没有人回答。
没有偏导数，最小二乘法就结束了。
所以我选择了一条更实用的道路，并在 C 中实现了一个强力函数，尝试了所有可能的参数三元组和 return 最佳匹配。
3 暴力破解函数
考虑到问题的性质，结果是样本数量为O(n^2)。
算法进行如下：将样本集分为三部分，第一个是t1[=140=之前的点的部分，第二个是[=94之间的点=]t1 和 t2 以及 t2 之后的最后一个点。
第一部分仅用于计算sigma，sigma只是第1部分点的算术平均值。
t1和t2是通过一个循环计算的，t1设置为每一个可能的原始点集中的点，从第二个开始往前。
对于t1的每个选择，t2设置为t1之后的每个可能点.
在每次迭代中都会计算一个错误，如果它是有史以来的最小值，则保存使用的参数。
错误是计算机作为残差的绝对值，因为绝对值应该很快（肯定比平方快）并且它符合度量的目的。
4 代码
#include <stdio.h>
#include <math.h>
float point_on_curve(float sigma, float t1, float t2, float t)
{
float a,b,c,d, K;
if (t <= t1)
return sigma;
if (t >= t2)
return 0;
K = (t1-t2)*(t1-t2)*(t1-t2);
a = -2*sigma/K;
b = 3*sigma*(t1+t2)/K;
c = -6*sigma*t1*t2/K;
d = sigma*t2*t2*(3*t1-t2)/K;
return a*t*t*t + b*t*t + c*t + d;
}
float compute_error(float sigma, float t1, float t2, int s, int dx, int* data, int len)
{
float error=0;
unsigned int i;
for (i = 0; i < len; i++)
error += fabs(point_on_curve(sigma, t1, t2, s+i*dx)- data[i]);
return error;
}
/*
* s is the starting time of the samples set
* dx is the separation in time between two sample (a.k.a. sampling period)
* data is the array of samples
* len is the number of samples
* sigma, t1, t2 are pointers to output parameters computed
*
* return 1 if not enough (3) samples, 0 if everything went ok.
*/
int curve_fit(int s, int dx, int* data, unsigned int len, float* sigma, float* t1, float* t2)
{
float l_sigma = 0;
float l_t1, l_t2;
float sum = 0;
float min_error, cur_error;
char error_valid = 0;
unsigned int i, j;
if (len < 3)
return 1;
for (i = 0; i < len; i++)
{
/* Compute sigma as the average of points <= i */
sum += data[i];
l_sigma = sum/(i+1);
/* Set t1 as the point i+1 */
l_t1 = s+(i+1)*dx;
for (j = i+2; j < len; j++)
{
/* Set t2 as the points i+2, i+3, i+4, ... */
l_t2 = s+j*dx;
/* Compute the error */
cur_error = compute_error(l_sigma, l_t1, l_t2, s, dx, data, len);
if (cur_error < min_error || !error_valid)
{
error_valid = 1;
min_error = cur_error;
*sigma = l_sigma;
*t1 = l_t1;
*t2 = l_t2;
}
}
}
return 0;
}
int main()
{
float sigma, t1, t2;
int data[]={41, 40, 40, 41, 41, 42, 42, 43, 43, 43, 51, 76, 98, 104, 104, 103, 104, 105, 105, 107, 105};
unsigned int len = sizeof(data)/sizeof(int);
unsigned int i;
for (i = 0; i < len; i++)
data[i] -= 107; /* Subtract the max */
if (curve_fit(1,1,data, len, &sigma, &t1, &t2))
printf("Not enough data!\n");
else
printf("Parameters: sigma = %.3f, t1 = %.3f, t2 = %.3f\n", sigma, t1, t2);
return 0;
}
请注意 Cl(t) 被定义为右极限为 0，因此代码假设这是案件。
这就是为什么从每个样本中减去最大值（107）的原因，我使用了Cl(t)[=140=的定义] 在开始时给出，直到后来才注意到样本有偏差。
现在我不打算调整代码，您可以轻松地在问题中添加另一个参数并在需要时从 1 开始重做这些步骤，或者简单地使用最大值转换样本。
代码的输出是
Parameters: sigma = -65.556, t1 = 10.000, t2 = 14.000
与给定的点集相匹配，考虑到它垂直平移了-107。
我需要确定光照变化的参数，它由这个连续分段多项式 C(t) 定义，其中 f(t) 是由两个边界点 (t1,c) 定义的三次曲线和 (t2,0)，还有 f'(t1)=0 和 f'(t2)=0。 Original Paper: Texture-Consistent Shadow Removal
强度曲线是从阴影边界的法线采样的，看起来像这样：
每一行都是样本，显示光照 change.So X 是列数，Y 是像素强度。
我有这样的真实数据（一个样本从所有样本中平均得出）：
我总共有N个样本，我需要确定参数(c,t1,t2)
我该怎么做？
我尝试通过在 Matlab 中求解线性方程来实现：
avr_curve为平均曲线，对所有样本取平均值
f(x)= x^3+a2*x^2+a1*x1+a0 是三次函数
%t1,t2 selected by hand
t1= 10;
t2= 15;
offset=10;
avr_curve= [41, 40, 40, 41, 41, 42, 42, 43, 43, 43, 51, 76, 98, 104, 104, 103, 104, 105, 105, 107, 105];
%gradx= convn(avr_curve,[-1 1],'same');
A= zeros(2*offset+1,3);
%b= zeros(2*offset+1,1);
b= avr_curve';
%for i= 1:2*offset+1
for i=t1:t2
i
x= i-offset-1
A(i,1)= x^2; %a2
A(i,2)= x; %a1
A(i,3)= 1; %a0
b(i,1)= b(i,1)-x^3;
end
u= A\b;
figure,plot(avr_curve(t1:t2))
%estimated cubic curve
for i= 1:2*offset+1
x= i-offset-1;
fx(i)=x^3+u(1)*x^2+u(2)*x+u(3);
end
figure,plot(fx(t1:t2))
[t1 t2] avr_curve 的一部分
我得到的三次曲线（不像avr_curve）
所以我做错了什么？
更新： 似乎我的错误是由于我使用 3 个变量对三次多项式进行建模，如下所示：
f(x)= x^3+a2*x^2+a1*x1+a0 - 3 variables
但后来我使用了 4 个变量，一切似乎都正常：
f(x)= a3*x^3+a2*x^2+a1*x1+a0 - 4 variables
这是 Matlab 中的代码：
%defined by hand
t1= 10;
t2= 14;
avr_curve= [41, 40, 40, 41, 41, 42, 42, 43, 43, 43, 51, 76, 98, 104, 104, 103, 104, 105, 105, 107, 105];
x= [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21];
%x= [-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]; %real x axis
%%%model 1
%%f(x)= x^3+a2*x^2+a1*x1+a0 - 3 variables
%A= zeros(4,3);
%b= [43 104]';
%%cubic equation at t1
%A(1,1)= t1^2; %a2
%A(1,2)= t1; %a1
%A(1,3)= 1; %a0
%b(1,1)= b(1,1)-t1^3;
%%cubic equation at t2
%A(2,1)= t2^2; %a2
%A(2,2)= t2; %a1
%A(2,3)= 1; %a0
%b(2,1)= b(2,1)-t1^3;
%%1st derivative at t1
%A(3,1)= 2*t1; %a2
%A(3,2)= 1; %a1
%A(3,3)= 0; %a0
%b(3,1)= -3*t1^2;
%%1st derivative at t2
%A(4,1)= 2*t2; %a2
%A(4,2)= 1; %a1
%A(4,3)= 0; %a0
%b(4,1)= -3*t2^2;
%model 2
%f(x)= a3*x^3+a2*x^2+a1*x1+a0 - 4 variables
A= zeros(4,4);
b= [43 104]';
%cubic equation at t1
A(1,1)= t1^3; %a3
A(1,2)= t1^2; %a2
A(1,3)= t1; %a1
A(1,4)= 1; %a0
b(1,1)= b(1,1);
%cubic equation at t2
A(2,1)= t2^3; %a3
A(2,2)= t2^2; %a2
A(2,3)= t2; %a1
A(2,4)= 1; %a0
b(2,1)= b(2,1);
%1st derivative at t1
A(3,1)= 3*t1^2; %a3
A(3,2)= 2*t1; %a2
A(3,3)= 1; %a1
A(3,4)= 0; %a0
b(3,1)= 0;
%1st derivative at t2
A(4,1)= 3*t2^2; %a3
A(4,2)= 2*t2; %a2
A(4,3)= 1; %a1
A(4,4)= 0; %a0
b(4,1)= 0;
size(A)
size(b)
u= A\b;
u
%estimated cubic curve
%dx=[1:21]; % global view
dx=t1-1:t2+1; % local view in [t1 t2]
for x= dx
%fx(x)=x^3+u(1)*x^2+u(2)*x+u(3); % model 1
fx(x)= u(1)*x^3+u(2)*x^2+u(3)*x+u(4); % model 2
end
err= 0;
for x= dx
err= err+(fx(x)-avr_curve(x))^2;
end
err
figure,plot(dx,avr_curve(dx),dx,fx(dx))
样条区间 [t1-1 t2+1]
并且在整个间隔
免责声明
我无法对下面给出的代码或方法的正确性提供任何保证，在使用任何这些之前请务必使用您的批判性判断力。
0。定义问题
你有这个分段定义的函数
其中f(t)为三次函数，为了唯一标识它，给出如下附加条件
您想恢复参数 t1、t2 和 sigma 的最佳值最小化给定点集的误差。
这本质上是 least squares 意义上的曲线拟合。
1 参数化 f(t) 三次函数
为了计算候选 Cl(t) 函数与点集之间的误差，我们需要计算 f(t)，它的一般形式（立方体）是
看来我们还有四个额外的参数需要考虑。实际上，此参数完全由 free 三个参数 t1、t2 和 定义西格玛.
重要的是不要将 free 参数与相关参数混淆。
给定 f(t) 的附加条件，我们可以建立这个线性系统
这告诉我们如何在给定三个自由参数的情况下计算立方体的参数。
这样Cl(t)就完全确定了，现在是时候寻找最佳人选了
2 最小化误差
我现在通常会选择最小二乘法。
由于这不是线性函数，因此没有用于计算 sigma、t1 和 t2 的封闭形式.
但是有数值方法，比如 Gauss-Newton 一个。
然而，需要以一种或另一种方式计算关于三个参数的偏导数。
我不知道如何计算像 t1 这样的分离参数的导数。
我搜索了 MathSE，发现 this 个解决相同问题的问题，但没有人回答。
没有偏导数，最小二乘法就结束了。
所以我选择了一条更实用的道路，并在 C 中实现了一个强力函数，尝试了所有可能的参数三元组和 return 最佳匹配。
3 暴力破解函数
考虑到问题的性质，结果是样本数量为O(n^2)。
算法进行如下：将样本集分为三部分，第一个是t1[=140=之前的点的部分，第二个是[=94之间的点=]t1 和 t2 以及 t2 之后的最后一个点。
第一部分仅用于计算sigma，sigma只是第1部分点的算术平均值。
t1和t2是通过一个循环计算的，t1设置为每一个可能的原始点集中的点，从第二个开始往前。
对于t1的每个选择，t2设置为t1之后的每个可能点.
在每次迭代中都会计算一个错误，如果它是有史以来的最小值，则保存使用的参数。
错误是计算机作为残差的绝对值，因为绝对值应该很快（肯定比平方快）并且它符合度量的目的。
4 代码
#include <stdio.h>
#include <math.h>
float point_on_curve(float sigma, float t1, float t2, float t)
{
float a,b,c,d, K;
if (t <= t1)
return sigma;
if (t >= t2)
return 0;
K = (t1-t2)*(t1-t2)*(t1-t2);
a = -2*sigma/K;
b = 3*sigma*(t1+t2)/K;
c = -6*sigma*t1*t2/K;
d = sigma*t2*t2*(3*t1-t2)/K;
return a*t*t*t + b*t*t + c*t + d;
}
float compute_error(float sigma, float t1, float t2, int s, int dx, int* data, int len)
{
float error=0;
unsigned int i;
for (i = 0; i < len; i++)
error += fabs(point_on_curve(sigma, t1, t2, s+i*dx)- data[i]);
return error;
}
/*
* s is the starting time of the samples set
* dx is the separation in time between two sample (a.k.a. sampling period)
* data is the array of samples
* len is the number of samples
* sigma, t1, t2 are pointers to output parameters computed
*
* return 1 if not enough (3) samples, 0 if everything went ok.
*/
int curve_fit(int s, int dx, int* data, unsigned int len, float* sigma, float* t1, float* t2)
{
float l_sigma = 0;
float l_t1, l_t2;
float sum = 0;
float min_error, cur_error;
char error_valid = 0;
unsigned int i, j;
if (len < 3)
return 1;
for (i = 0; i < len; i++)
{
/* Compute sigma as the average of points <= i */
sum += data[i];
l_sigma = sum/(i+1);
/* Set t1 as the point i+1 */
l_t1 = s+(i+1)*dx;
for (j = i+2; j < len; j++)
{
/* Set t2 as the points i+2, i+3, i+4, ... */
l_t2 = s+j*dx;
/* Compute the error */
cur_error = compute_error(l_sigma, l_t1, l_t2, s, dx, data, len);
if (cur_error < min_error || !error_valid)
{
error_valid = 1;
min_error = cur_error;
*sigma = l_sigma;
*t1 = l_t1;
*t2 = l_t2;
}
}
}
return 0;
}
int main()
{
float sigma, t1, t2;
int data[]={41, 40, 40, 41, 41, 42, 42, 43, 43, 43, 51, 76, 98, 104, 104, 103, 104, 105, 105, 107, 105};
unsigned int len = sizeof(data)/sizeof(int);
unsigned int i;
for (i = 0; i < len; i++)
data[i] -= 107; /* Subtract the max */
if (curve_fit(1,1,data, len, &sigma, &t1, &t2))
printf("Not enough data!\n");
else
printf("Parameters: sigma = %.3f, t1 = %.3f, t2 = %.3f\n", sigma, t1, t2);
return 0;
}
请注意 Cl(t) 被定义为右极限为 0，因此代码假设这是案件。
这就是为什么从每个样本中减去最大值（107）的原因，我使用了Cl(t)[=140=的定义] 在开始时给出，直到后来才注意到样本有偏差。
现在我不打算调整代码，您可以轻松地在问题中添加另一个参数并在需要时从 1 开始重做这些步骤，或者简单地使用最大值转换样本。
代码的输出是
Parameters: sigma = -65.556, t1 = 10.000, t2 = 14.000
与给定的点集相匹配，考虑到它垂直平移了-107。