Cycles: svn merge -r41225:41232 ^/trunk/blender
[blender.git] / intern / cycles / subd / subd_patch.cpp
1 /*
2  * Copyright 2011, Blender Foundation.
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18
19 /* Parts adapted from code in the public domain in NVidia Mesh Tools. */
20
21 #include "mesh.h"
22
23 #include "subd_patch.h"
24
25 #include "util_math.h"
26 #include "util_types.h"
27
28 CCL_NAMESPACE_BEGIN
29
30 /* De Casteljau Evaluation */
31
32 static float3 decasteljau_quadratic(float t, const float3 cp[3])
33 {
34         float3 d0 = cp[0] + t*(cp[1] - cp[0]);
35         float3 d1 = cp[1] + t*(cp[2] - cp[1]);
36
37         return d0 + t*(d1 - d0);
38 }
39
40 static void decasteljau_cubic(float3 *P, float3 *dt, float t, const float3 cp[4])
41 {
42         float3 d0 = cp[0] + t*(cp[1] - cp[0]);
43         float3 d1 = cp[1] + t*(cp[2] - cp[1]);
44         float3 d2 = cp[2] + t*(cp[3] - cp[2]);
45
46         d0 += t*(d1 - d0);
47         d1 += t*(d2 - d1);
48
49         *P = d0 + t*(d1 - d0);
50         if(dt) *dt = d1 - d0;
51 }
52
53 static void decasteljau_bicubic(float3 *P, float3 *du, float3 *dv, const float3 cp[16], float u, float v)
54 {
55         float3 ucp[4], utn[4];
56
57         /* interpolate over u */
58         decasteljau_cubic(ucp+0, utn+0, u, cp);
59         decasteljau_cubic(ucp+1, utn+1, u, cp+4);
60         decasteljau_cubic(ucp+2, utn+2, u, cp+8);
61         decasteljau_cubic(ucp+3, utn+3, u, cp+12);
62
63         /* interpolate over v */
64         decasteljau_cubic(P, dv, v, ucp);
65         if(du) decasteljau_cubic(du, NULL, v, utn);
66 }
67
68 static float3 decasteljau_tangent(const float3 cp[12], float u, float v)
69 {
70         float3 ucp[3];
71
72         decasteljau_cubic(ucp+0, NULL, v, cp);
73         decasteljau_cubic(ucp+1, NULL, v, cp+4);
74         decasteljau_cubic(ucp+2, NULL, v, cp+8);
75
76         return decasteljau_quadratic(u, ucp);
77 }
78
79 /* Linear Quad Patch */
80
81 void LinearQuadPatch::eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v)
82 {
83         float3 d0 = interp(hull[0], hull[1], u);
84         float3 d1 = interp(hull[2], hull[3], u);
85
86         *P = interp(d0, d1, v);
87
88         if(dPdu && dPdv) {
89                 *dPdu = interp(hull[1] - hull[0], hull[3] - hull[2], v);
90                 *dPdv = interp(hull[2] - hull[0], hull[3] - hull[1], u);
91         }
92 }
93
94 BoundBox LinearQuadPatch::bound()
95 {
96         BoundBox bbox;
97
98         for(int i = 0; i < 4; i++)
99                 bbox.grow(hull[i]);
100         
101         return bbox;
102 }
103
104 /* Linear Triangle Patch */
105
106 void LinearTrianglePatch::eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v)
107 {
108         *P = u*hull[0] + v*hull[1] + (1.0f - u - v)*hull[2];
109
110         if(dPdu && dPdv) {
111                 *dPdu = hull[0] - hull[2];
112                 *dPdv = hull[1] - hull[2];
113         }
114 }
115
116 BoundBox LinearTrianglePatch::bound()
117 {
118         BoundBox bbox;
119
120         for(int i = 0; i < 3; i++)
121                 bbox.grow(hull[i]);
122         
123         return bbox;
124 }
125
126 /* Bicubic Patch */
127
128 void BicubicPatch::eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v)
129 {
130         decasteljau_bicubic(P, dPdu, dPdv, hull, u, v);
131 }
132
133 BoundBox BicubicPatch::bound()
134 {
135         BoundBox bbox;
136
137         for(int i = 0; i < 16; i++)
138                 bbox.grow(hull[i]);
139         
140         return bbox;
141 }
142
143 /* Bicubic Patch with Tangent Fields */
144
145 void BicubicTangentPatch::eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v)
146 {
147         decasteljau_bicubic(P, NULL, NULL, hull, u, v);
148
149         if(dPdu) *dPdu = decasteljau_tangent(utan, u, v);
150         if(dPdv) *dPdv = decasteljau_tangent(vtan, v, u);
151 }
152
153 BoundBox BicubicTangentPatch::bound()
154 {
155         BoundBox bbox;
156
157         for(int i = 0; i < 16; i++)
158                 bbox.grow(hull[i]);
159         
160         return bbox;
161 }
162
163 /* Gregory Patch */
164
165 static float no_zero_div(float f)
166 {
167         if(f == 0.0f) return 1.0f;
168         return f;
169 }
170
171 void GregoryQuadPatch::eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v)
172 {
173         float3 bicubic[16];
174
175         float U = 1 - u;
176         float V = 1 - v;
177
178         /*  8     9     10     11
179          * 12   0\1     2/3    13
180          * 14   4/5     6\7    15
181          * 16    17     18     19
182          */
183
184         bicubic[5] = (u*hull[1] + v*hull[0])/no_zero_div(u + v);
185         bicubic[6] = (U*hull[2] + v*hull[3])/no_zero_div(U + v);
186         bicubic[9] = (u*hull[5] + V*hull[4])/no_zero_div(u + V);
187         bicubic[10] = (U*hull[6] + V*hull[7])/no_zero_div(U + V);
188
189         // Map gregory control points to bezier control points.
190         bicubic[0] = hull[8];
191         bicubic[1] = hull[9];
192         bicubic[2] = hull[10];
193         bicubic[3] = hull[11];
194         bicubic[4] = hull[12];
195         bicubic[7] = hull[13];
196         bicubic[8] = hull[14];
197         bicubic[11] = hull[15];
198         bicubic[12] = hull[16];
199         bicubic[13] = hull[17];
200         bicubic[14] = hull[18];
201         bicubic[15] = hull[19];
202
203         decasteljau_bicubic(P, dPdu, dPdv, bicubic, u, v);
204 }
205
206 BoundBox GregoryQuadPatch::bound()
207 {
208         BoundBox bbox;
209
210         for(int i = 0; i < 20; i++)
211                 bbox.grow(hull[i]);
212         
213         return bbox;
214 }
215
216 void GregoryTrianglePatch::eval(float3 *P, float3 *dPdu, float3 *dPdv, float u, float v)
217 {
218         /*                    6
219          *                    
220          *              14   0/1   7
221          *                                                
222          *    13   5/4     3\2   8
223          *                     
224          * 12      11       10      9
225          */
226
227         float w = 1 - u - v;
228         float uu = u * u;
229         float vv = v * v;
230         float ww = w * w;
231         float uuu = uu * u;
232         float vvv = vv * v;
233         float www = ww * w;
234
235         float U = 1 - u;
236         float V = 1 - v;
237         float W = 1 - w;
238
239         float3 C0 = ( v*U * hull[5] + u*V * hull[4] ) / no_zero_div(v*U + u*V);
240         float3 C1 = ( w*V * hull[3] + v*W * hull[2] ) / no_zero_div(w*V + v*W);
241         float3 C2 = ( u*W * hull[1] + w*U * hull[0] ) / no_zero_div(u*W + w*U);
242
243         *P =
244                 (hull[12] * www + 3*hull[11] * ww*u + 3*hull[10] * w*uu + hull[ 9]*uuu) * (w + u) +
245                 (hull[ 9] * uuu + 3*hull[ 8] * uu*v + 3*hull[ 7] * u*vv + hull[ 6]*vvv) * (u + v) +
246                 (hull[ 6] * vvv + 3*hull[14] * vv*w + 3*hull[13] * v*ww + hull[12]*www) * (v + w) -
247                 (hull[12] * www*w + hull[ 9] * uuu*u + hull[ 6] * vvv*v) +
248                 12*(C0 * u*v*ww + C1 * uu*v*w   + C2 * u*vv*w);
249
250         if(dPdu || dPdv) {
251                 float3 E1 = (hull[12]*www + 3*hull[11]*ww*u + 3*hull[10]*w*uu + hull[ 9]*uuu);
252                 float3 E2 = (hull[ 9]*uuu + 3*hull[ 8]*uu*v + 3*hull[ 7]*u*vv + hull[ 6]*vvv);
253                 float3 E3 = (hull[ 6]*vvv + 3*hull[14]*vv*w + 3*hull[13]*v*ww + hull[12]*www);
254
255                 if(dPdu) {
256                         float3 E1u = 3*( - hull[12]*ww + hull[11]*(ww-2*u*w) +   hull[10]*(2*u*w-uu) + hull[ 9]*uu);
257                         float3 E2u = 3*(   hull[ 9]*uu + 2*hull[ 8]*u*v      +   hull[ 7]*vv             );
258                         float3 E3u = 3*(                    - hull[14]*vv                - 2*hull[13]*v*w               - hull[12]*ww);
259                         float3 Su  = 4*( -hull[12]*www + hull[9]*uuu);
260                         float3 Cu  = 12*( C0*(ww*v-2*u*v*w) + C1*(2*u*v*w-uu*v) + C2*vv*(w-u) );
261
262                         *dPdu = E1u*(w+u) + (E2+E2u*(u+v)) + (E3u*(v+w)-E3) - Su + Cu;
263                 }
264
265                 if(dPdv) {
266                         float3 E1v = 3*(-hull[12]*ww  - 2*hull[11]*w*u       -   hull[10]*uu             );
267                         float3 E2v = 3*(                      hull[ 8]*uu                + 2*hull[ 7]*u*v               + hull[ 6]*vv);
268                         float3 E3v = 3*( hull[ 6]*vv  +  hull[14]*(2*w*v-vv) +   hull[13]*(ww-2*w*v) - hull[12]*ww);
269                         float3 Sv  = 4*(-hull[12]*www +  hull[ 6]*vvv);
270                         float3 Cv  = 12*(C0*(u*ww-2*u*v*w) + C1*uu*(w-v) + C2*(2*u*v*w-u*vv));
271
272                         *dPdv = ((E1v*(w+u)-E1) + (E2+E2v*(u+v)) + E3v*(v+w) - Sv + Cv );
273                 }
274         }
275 }
276
277 BoundBox GregoryTrianglePatch::bound()
278 {
279         BoundBox bbox;
280
281         for(int i = 0; i < 20; i++)
282                 bbox.grow(hull[i]);
283         
284         return bbox;
285 }
286
287 CCL_NAMESPACE_END
288