1 | % Many thanks to Kurt Ditzel
|
---|
2 | % Acuity Research for this bug fix...
|
---|
3 |
|
---|
4 |
|
---|
5 | function p = calc_cubic_coefs(f1, f2)
|
---|
6 | %
|
---|
7 | % use symbolic solver
|
---|
8 | %
|
---|
9 |
|
---|
10 | c11 = f1(1);
|
---|
11 | c12 = f1(2);
|
---|
12 | c13 = f1(3);
|
---|
13 | c21 = f1(4);
|
---|
14 | c22 = f1(5);
|
---|
15 | c23 = f1(6);
|
---|
16 | c31 = f1(7);
|
---|
17 | c32 = f1(8);
|
---|
18 | c33 = f1(9);
|
---|
19 |
|
---|
20 | d11 = f2(1);
|
---|
21 | d12 = f2(2);
|
---|
22 | d13 = f2(3);
|
---|
23 | d21 = f2(4);
|
---|
24 | d22 = f2(5);
|
---|
25 | d23 = f2(6);
|
---|
26 | d31 = f2(7);
|
---|
27 | d32 = f2(8);
|
---|
28 | d33 = f2(9);
|
---|
29 |
|
---|
30 | %syms c11 c12 c13 c21 c22 c23 c31 c32 c33;
|
---|
31 | %syms d11 d12 d13 d21 d22 d23 d31 d32 d33;
|
---|
32 | %syms f1 f2 a X Y D T;
|
---|
33 |
|
---|
34 | %f1 = [c11, c12, c13, c21, c22, c23, c31, c32, c33];
|
---|
35 | %f2 = [d11, d12, d13, d21, d22, d23, d31, d32, d33];
|
---|
36 | %
|
---|
37 | % disp('X')
|
---|
38 | % X = a*f1 + (1-a)*f2;
|
---|
39 | % pretty(X)
|
---|
40 | % disp('Y');
|
---|
41 | % Y = reshape(X,3,3).';
|
---|
42 | % pretty(Y)
|
---|
43 | % D = det(Y);
|
---|
44 | % disp('D');
|
---|
45 | % %pretty(D);
|
---|
46 | %
|
---|
47 | % disp('Collected coeffs of a')
|
---|
48 | % T = collect(simplify(D),a)
|
---|
49 |
|
---|
50 | % solve for given the points
|
---|
51 | %
|
---|
52 | % p =
|
---|
53 | p = [(d11*c22*d33+c21*d13*d32-c21*d13*c32-c21*c13*d32+c21*c13*c32-c21*d12*d33+c21*d12*c33+c21*c12*d33-c21*c12*c33+c11*c22*c33-c11*c22*d33+d11*c23*c32-d11*c23*d32-d11*d23*c32-d11*c22*c33-c11*d22*c33-c11*c23*c32+d11*d22*c33-d21*d13*d32+c11*c23*d32+c11*d22*d33+d21*d12*d33-d31*d13*c22-d31*c13*d22+d31*c13*c22+d31*d12*c23+d31*c12*d23-d31*c12*c23-c31*d13*d22+c31*d13*c22+c31*c13*d22-c31*c13*c22+c31*d12*d23-c31*d12*c23-c31*c12*d23+c31*c12*c23+d21*d13*c32+d21*c13*d32-d21*c13*c32-d21*d12*c33-d21*c12*d33+d21*c12*c33+d31*d13*d22-d31*d12*d23+c11*d23*c32-d11*d22*d33+d11*d23*d32-c11*d23*d32),
|
---|
54 | +(-2*d11*c22*d33-2*c21*d13*d32+c21*d13*c32+c21*c13*d32+2*c21*d12*d33-c21*d12*c33-c21*c12*d33+c11*c22*d33-d11*c23*c32+2*d11*c23*d32+2*d11*d23*c32+d11*c22*c33+c11*d22*c33-2*d11*d22*c33+3*d21*d13*d32-c11*c23*d32-2*c11*d22*d33-3*d21*d12*d33+2*d31*d13*c22+2*d31*c13*d22-d31*c13*c22-2*d31*d12*c23-2*d31*c12*d23+d31*c12*c23+2*c31*d13*d22-c31*d13*c22-c31*c13*d22-2*c31*d12*d23+c31*d12*c23+c31*c12*d23-2*d21*d13*c32-2*d21*c13*d32+d21*c13*c32+2*d21*d12*c33+2*d21*c12*d33-d21*c12*c33-3*d31*d13*d22+3*d31*d12*d23-c11*d23*c32+3*d11*d22*d33-3*d11*d23*d32+2*c11*d23*d32),
|
---|
55 | +(-3*d21*d13*d32-c11*d23*d32+3*d11*d23*d32-c21*d12*d33-c31*d13*d22+d11*c22*d33+d21*d13*c32+c21*d13*d32-d11*d23*c32+3*d21*d12*d33+c11*d22*d33+d11*d22*c33-d21*c12*d33+d21*c13*d32-3*d11*d22*d33-d11*c23*d32-d31*c13*d22+d31*c12*d23+3*d31*d13*d22-d31*d13*c22+c31*d12*d23-3*d31*d12*d23-d21*d12*c33+d31*d12*c23),
|
---|
56 | -d21*d12*d33+d11*d22*d33-d11*d23*d32+d21*d13*d32-d31*d13*d22+d31*d12*d23];
|
---|