computer-vision - 三次样条/曲线拟合

Question

我需要确定照明变化的参数，它由这个连续分段多项式 C(t) 定义，其中 f(t) 是由两个边界点 (t1,c) 和 (t2,0) 定义的三次曲线)，也 f'(t1)=0 和 f'(t2)=0。 原始论文：纹理一致的阴影去除

强度曲线是从阴影边界上的法线采样的，它看起来像这样：

每一行都是样本，显示光照变化。所以 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​​] 上的样条

并且全间隔

Answer 1

免责声明

我不能对下面给出的代码或方法的正确性做出任何保证，在使用任何代码或方法之前，请始终使用您的批判性意识。

0. 定义问​​题

你有这个分段定义的函数

其中f(t)是三次函数，为了唯一标识它，给出以下附加条件

您希望恢复参数t1、t2和sigma的最佳值，以最小化给定点集的误差。

这本质上是最小二乘意义上的曲线拟合。

1 参数化f(t)三次函数

为了计算候选Cl(t)函数和我们需要计算f(t)的点集之间的误差，它的一般形式（三次）是

因此，我们似乎有四个额外的参数需要考虑。实际上，这个参数完全由三个自由参数t1、t2和sigma定义。
重要的是不要将自由参数与依赖参数混淆。

给定f(t)的附加条件，我们可以建立这个线性系统

其中有一个解决方案（如预期的那样）由

这告诉我们如何在给定三个自由参数的情况下计算三次参数。
这样Cl(t)就完全确定了，现在是时候找到最佳候选者了。

2 最小化错误

我现在通常会选择最小二乘。
由于这不是线性函数，因此没有用于计算sigma、t1和t2的封闭形式。
然而，也有数值方法，例如高斯-牛顿法。

然而，需要以一种或另一种方式计算关于三个参数的偏导数。
我不知道如何计算像t1这样的分离参数的导数。

我搜索了 MathSE，发现这个问题解决了同样的问题，但是没有人回答。

没有偏导数，最小二乘法就结束了。

所以我采取了一条更实际的道路，并在 C 中实现了一个蛮力函数，它尝试所有可能的三元组参数并返回最佳匹配。

3 蛮力功能

鉴于问题的性质，这在样本数量上是O (n^2)。

算法进行如下：将样本集分为三部分，第一部分是t1之前的点，第二部分是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)的定义，只是后来才注意到样本有偏差。

现在我不打算修改代码，您可以轻松地在问题中添加另一个参数，并在需要时从 1 开始重做步骤，或者简单地使用最大值转换样本。

代码的输出是

 Parameters: sigma = -65.556, t1 = 10.000, t2 = 14.000

考虑到它垂直平移了 -107，它与给定的点集相匹配。