/src/reduction.c

 1#include "reduction.h"
 2
 3float** bezier_reduction_matrix(int n, int m)
 4{
 5  assert(n > m);
 6
 7  float (*psi)[n+1] = malloc(sizeof(float) * (m + 1) * (n + 1));
 8
 9  float C = factrl(n)*up_power(1-n, m)
10    / (factrl(m)*up_power(m+2, n));
11
12  // 0,0
13  {
14    psi[0][0] = 0;
15    for(int h = 0; h <= m; ++h)
16    {
17      float t = bico(m,h) * up_power(m+2,h) / ((h-n)*factrl(h));
18      if(h % 2 != 0)
19	t *= -1;
20      psi[0][0] += t;
21    }
22    psi[0][0] *= -n*C;
23  }
24
25  // 0,1
26  {
27    psi[0][1] = (m+1)*C;
28    if(m%2 != 0)
29      psi[0][1] *= -1;
30  }
31
32  // first row
33  {
34    for(int j = 1; j <= n - 1; ++j)
35    {
36      float F = ((float)((j-1)*(n-j+1)))/((float) ((n-j)*(j+1)));
37      float E = 1. + F - ((float)(m*(m+2)))/((float) ((n-j)*(j+1)));
38
39      psi[0][j+1] = E*psi[0][j] - F*psi[0][j-1];
40    }
41  }
42
43  // rest of matrix
44  {
45    for(int i = 0; i <= m - 1; ++i)
46    {
47      for(int j = 0; j <= n; ++j)
48      {
49	float Anj = j*(j-n-1);
50	float Ami = i*(i-m-1);
51	float Bnj = (j-n)*(j+1);
52	float Bmi = (i-m)*(i+1);
53	float Cnj = Anj + Bnj;
54	float Cmi = Ami + Bmi;
55
56	float pp0m1 = ((j - 1) < 0) ? 0 : psi[i][j-1];
57	float pp0p0 = psi[i][j];
58	float pm1p0 = ((i - 1) < 0) ? 0 : psi[i-1][j];
59	float pp0p1 = ((j + 1) > n) ? 0 : psi[i][j+1];
60
61	psi[i+1][j] = (Anj*pp0m1 + (Cmi-Cnj)*pp0p0 +
62	  Bnj*pp0p1 - Ami*pm1p0) / Bmi;
63      }
64    }
65  }
66
67  return (float**)psi;
68}
69
70Bezier* bezier_degree_reduction(Bezier* c, int m, float **reduction_matrix)
71{
72  int n = c->n;
73
74  float (*psi)[n+1] = (float(*)[])reduction_matrix;
75  float (*e) = malloc(sizeof(float) * (n + 1));
76
77  for(int i = 0; i <= m; ++i)
78  {
79    e[i] = 0;
80    for(int j = 0; j <= n; ++j)
81      e[i] += c->c[j]*psi[i][j];
82  }
83
84  Bezier* low = bezier_create_with_coeffs(m, e);
85  low->dom = interval_copy(c->dom);
86  
87  return low;
88}
89