Merged changes in the trunk up to revision 46557.
[blender.git] / source / blender / blenkernel / intern / ocean.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
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  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
19  * All rights reserved.
20  *
21  * Contributors: Matt Ebb, Hamed Zaghaghi
22  * Based on original code by Drew Whitehouse / Houdini Ocean Toolkit
23  * OpenMP hints by Christian Schnellhammer
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27
28
29 #include <math.h>
30 #include <stdlib.h>
31
32 #include <string.h>
33
34 #include "MEM_guardedalloc.h"
35
36 #include "DNA_scene_types.h"
37
38 #include "BKE_image.h"
39 #include "BKE_ocean.h"
40 #include "BKE_utildefines.h"
41
42 #include "BKE_global.h" // XXX TESTING
43
44 #include "BLI_math_base.h"
45 #include "BLI_math_inline.h"
46 #include "BLI_rand.h"
47 #include "BLI_string.h"
48 #include "BLI_threads.h"
49 #include "BLI_path_util.h"
50 #include "BLI_utildefines.h"
51
52 #include "IMB_imbuf.h"
53 #include "IMB_imbuf_types.h"
54
55 #include "RE_render_ext.h"
56
57 #ifdef WITH_OCEANSIM
58
59 // Ocean code
60 #include "fftw3.h"
61
62 #define GRAVITY  9.81f
63
64 typedef struct Ocean {
65         /* ********* input parameters to the sim ********* */
66         float _V;
67         float _l;
68         float _w;
69         float _A;
70         float _damp_reflections;
71         float _wind_alignment;
72         float _depth;
73
74         float _wx;
75         float _wz;
76
77         float _L;
78
79         /* dimensions of computational grid */
80         int _M;
81         int _N;
82
83         /* spatial size of computational grid */
84         float _Lx;
85         float _Lz;
86
87         float normalize_factor;                                 // init w
88         float time;
89
90         short _do_disp_y;
91         short _do_normals;
92         short _do_chop;
93         short _do_jacobian;
94
95         /* mutex for threaded texture access */
96         ThreadRWMutex oceanmutex;
97
98         /* ********* sim data arrays ********* */
99
100         /* two dimensional arrays of complex */
101         fftw_complex *_fft_in;                  // init w       sim w
102         fftw_complex *_fft_in_x;                        // init w       sim w
103         fftw_complex *_fft_in_z;                        // init w       sim w
104         fftw_complex *_fft_in_jxx;                      // init w       sim w
105         fftw_complex *_fft_in_jzz;                      // init w       sim w
106         fftw_complex *_fft_in_jxz;                      // init w       sim w
107         fftw_complex *_fft_in_nx;                       // init w       sim w
108         fftw_complex *_fft_in_nz;                       // init w       sim w
109         fftw_complex *_htilda;                  // init w       sim w (only once)
110
111         /* fftw "plans" */
112         fftw_plan _disp_y_plan;                 // init w       sim r
113         fftw_plan _disp_x_plan;                 // init w       sim r
114         fftw_plan _disp_z_plan;                 // init w       sim r
115         fftw_plan _N_x_plan;                    // init w       sim r
116         fftw_plan _N_z_plan;                    // init w       sim r
117         fftw_plan _Jxx_plan;                    // init w       sim r
118         fftw_plan _Jxz_plan;                    // init w       sim r
119         fftw_plan _Jzz_plan;                    // init w       sim r
120
121         /* two dimensional arrays of float */
122         double *  _disp_y;                              // init w       sim w via plan?
123         double * _N_x;                                  // init w       sim w via plan?
124         /*float * _N_y; all member of this array has same values, so convert this array to a float to reduce memory usage (MEM01)*/
125         double _N_y;                                    //                      sim w ********* can be rearranged?
126         double * _N_z;                                  // init w       sim w via plan?
127         double * _disp_x;                               // init w       sim w via plan?
128         double * _disp_z;                               // init w       sim w via plan?
129
130         /* two dimensional arrays of float */
131         /* Jacobian and minimum eigenvalue */
132         double * _Jxx;                                  // init w       sim w
133         double * _Jzz;                                  // init w       sim w
134         double * _Jxz;                                  // init w       sim w
135
136         /* one dimensional float array */
137         float * _kx;                                    // init w       sim r
138         float * _kz;                                    // init w       sim r
139
140         /* two dimensional complex array */
141         fftw_complex * _h0;                             // init w       sim r
142         fftw_complex * _h0_minus;               // init w       sim r
143
144         /* two dimensional float array */
145         float * _k;                                             // init w       sim r
146 } Ocean;
147
148
149
150 static float nextfr(float min, float max)
151 {
152         return BLI_frand()*(min-max)+max;
153 }
154
155 static float gaussRand (void)
156 {
157         float x;                // Note: to avoid numerical problems with very small
158         float y;                // numbers, we make these variables singe-precision
159         float length2;  // floats, but later we call the double-precision log()
160         // and sqrt() functions instead of logf() and sqrtf().
161         do
162         {
163                 x = (float) (nextfr (-1, 1));
164                 y = (float)(nextfr (-1, 1));
165                 length2 = x * x + y * y;
166         }
167         while (length2 >= 1 || length2 == 0);
168
169         return x * sqrtf(-2.0f * logf(length2) / length2);
170 }
171
172 /**
173  * Some useful functions
174  * */
175 MINLINE float lerp(float a, float b, float f)
176 {
177         return a + (b-a)*f;
178 }
179
180 MINLINE float catrom(float p0, float p1, float p2, float p3, float f)
181 {
182         return 0.5f *((2.0f * p1) +
183                       (-p0 + p2) * f +
184                       (2.0f*p0 - 5.0f*p1 + 4.0f*p2 - p3) * f*f +
185                       (-p0 + 3.0f*p1- 3.0f*p2 + p3) * f*f*f);
186 }
187
188 MINLINE float omega(float k, float depth)
189 {
190         return sqrt(GRAVITY*k * tanh(k*depth));
191 }
192
193 // modified Phillips spectrum
194 static float Ph(struct Ocean* o, float kx, float kz )
195 {
196         float tmp;
197         float k2 = kx*kx + kz*kz;
198
199         if (k2 == 0.0f) {
200                 return 0.0f; // no DC component
201         }
202
203         // damp out the waves going in the direction opposite the wind
204         tmp = (o->_wx * kx  + o->_wz * kz)/sqrtf(k2);
205         if (tmp < 0) {
206                 tmp *= o->_damp_reflections;
207         }
208
209         return o->_A * expf( -1.0f / (k2*(o->_L*o->_L))) * expf(-k2 * (o->_l*o->_l)) * powf(fabsf(tmp), o->_wind_alignment) / (k2*k2);
210 }
211
212 static void compute_eigenstuff(struct OceanResult *ocr, float jxx, float jzz, float jxz)
213 {
214         float a, b, qplus, qminus;
215         a = jxx + jzz;
216         b = sqrt((jxx - jzz)*(jxx - jzz) + 4 * jxz * jxz);
217
218         ocr->Jminus = 0.5f*(a-b);
219         ocr->Jplus  = 0.5f*(a+b);
220
221         qplus  = (ocr->Jplus  - jxx)/jxz;
222         qminus = (ocr->Jminus - jxx)/jxz;
223
224         a = sqrt(1 + qplus*qplus);
225         b = sqrt(1 + qminus*qminus);
226
227         ocr->Eplus[0] = 1.0f/ a;
228         ocr->Eplus[1] = 0.0f;
229         ocr->Eplus[2] = qplus/a;
230
231         ocr->Eminus[0] = 1.0f/b;
232         ocr->Eminus[1] = 0.0f;
233         ocr->Eminus[2] = qminus/b;
234 }
235
236 /*
237  * instead of Complex.h
238  * in fftw.h "fftw_complex" typedefed as double[2]
239  * below you can see functions are needed to work with such complex numbers.
240  * */
241 static void init_complex(fftw_complex cmpl, float real, float image)
242 {
243         cmpl[0] = real;
244         cmpl[1] = image;
245 }
246
247 #if 0   // unused
248 static void add_complex_f(fftw_complex res, fftw_complex cmpl, float f)
249 {
250         res[0] = cmpl[0] + f;
251         res[1] = cmpl[1];
252 }
253 #endif
254
255 static void add_comlex_c(fftw_complex res, fftw_complex cmpl1, fftw_complex cmpl2)
256 {
257         res[0] = cmpl1[0] + cmpl2[0];
258         res[1] = cmpl1[1] + cmpl2[1];
259 }
260
261 static void mul_complex_f(fftw_complex res, fftw_complex cmpl, float f)
262 {
263         res[0] = cmpl[0]*f;
264         res[1] = cmpl[1]*f;
265 }
266
267 static void mul_complex_c(fftw_complex res, fftw_complex cmpl1, fftw_complex cmpl2)
268 {
269         fftwf_complex temp;
270         temp[0] = cmpl1[0]*cmpl2[0]-cmpl1[1]*cmpl2[1];
271         temp[1] = cmpl1[0]*cmpl2[1]+cmpl1[1]*cmpl2[0];
272         res[0] = temp[0];
273         res[1] = temp[1];
274 }
275
276 static float real_c(fftw_complex cmpl)
277 {
278         return cmpl[0];
279 }
280
281 static float image_c(fftw_complex cmpl)
282 {
283         return cmpl[1];
284 }
285
286 static void conj_complex(fftw_complex res, fftw_complex cmpl1)
287 {
288         res[0] = cmpl1[0];
289         res[1] = -cmpl1[1];
290 }
291
292 static void exp_complex(fftw_complex res, fftw_complex cmpl)
293 {
294         float r = expf(cmpl[0]);
295
296         res[0] = cos(cmpl[1])*r;
297         res[1] = sin(cmpl[1])*r;
298 }
299
300 float BKE_ocean_jminus_to_foam(float jminus, float coverage)
301 {
302         float foam = jminus * -0.005f + coverage;
303         CLAMP(foam, 0.0f, 1.0f);
304         return foam*foam;
305 }
306
307 void BKE_ocean_eval_uv(struct Ocean *oc, struct OceanResult *ocr, float u, float v)
308 {
309         int i0, i1, j0, j1;
310         float frac_x, frac_z;
311         float uu, vv;
312
313         // first wrap the texture so 0 <= (u, v) < 1
314         u = fmodf(u, 1.0f);
315         v = fmodf(v, 1.0f);
316
317         if (u < 0) u += 1.0f;
318         if (v < 0) v += 1.0f;
319
320         BLI_rw_mutex_lock(&oc->oceanmutex, THREAD_LOCK_READ);
321
322         uu = u * oc->_M;
323         vv = v * oc->_N;
324
325         i0 = (int)floor(uu);
326         j0 = (int)floor(vv);
327
328         i1 = (i0 + 1);
329         j1 = (j0 + 1);
330
331         frac_x = uu - i0;
332         frac_z = vv - j0;
333
334         i0 = i0 % oc->_M;
335         j0 = j0 % oc->_N;
336
337         i1 = i1 % oc->_M;
338         j1 = j1 % oc->_N;
339
340
341 #define BILERP(m) (lerp(lerp(m[i0*oc->_N+j0], m[i1*oc->_N+j0], frac_x), lerp(m[i0*oc->_N+j1], m[i1*oc->_N+j1], frac_x), frac_z))
342         {
343                 if (oc->_do_disp_y) {
344                         ocr->disp[1] = BILERP(oc->_disp_y);
345                 }
346
347                 if (oc->_do_normals) {
348                         ocr->normal[0] = BILERP(oc->_N_x);
349                         ocr->normal[1] = oc->_N_y/*BILERP(oc->_N_y) (MEM01)*/;
350                         ocr->normal[2] = BILERP(oc->_N_z);
351                 }
352
353                 if (oc->_do_chop) {
354                         ocr->disp[0] = BILERP(oc->_disp_x);
355                         ocr->disp[2] = BILERP(oc->_disp_z);
356                 }
357                 else {
358                         ocr->disp[0] = 0.0;
359                         ocr->disp[2] = 0.0;
360                 }
361
362                 if (oc->_do_jacobian) {
363                         compute_eigenstuff(ocr, BILERP(oc->_Jxx), BILERP(oc->_Jzz), BILERP(oc->_Jxz));
364                 }
365         }
366 #undef BILERP
367
368         BLI_rw_mutex_unlock(&oc->oceanmutex);
369 }
370
371 // use catmullrom interpolation rather than linear
372 void BKE_ocean_eval_uv_catrom(struct Ocean *oc, struct OceanResult *ocr, float u, float v)
373 {
374         int i0, i1, i2, i3, j0, j1, j2, j3;
375         float frac_x, frac_z;
376         float uu, vv;
377
378         // first wrap the texture so 0 <= (u, v) < 1
379         u = fmod(u, 1.0f);
380         v = fmod(v, 1.0f);
381
382         if (u < 0) u += 1.0f;
383         if (v < 0) v += 1.0f;
384
385         BLI_rw_mutex_lock(&oc->oceanmutex, THREAD_LOCK_READ);
386
387         uu = u * oc->_M;
388         vv = v * oc->_N;
389
390         i1 = (int)floor(uu);
391         j1 = (int)floor(vv);
392
393         i2 = (i1 + 1);
394         j2 = (j1 + 1);
395
396         frac_x = uu - i1;
397         frac_z = vv - j1;
398
399         i1 = i1 % oc->_M;
400         j1 = j1 % oc->_N;
401
402         i2 = i2 % oc->_M;
403         j2 = j2 % oc->_N;
404
405         i0 = (i1-1);
406         i3 = (i2+1);
407         i0 = i0 <   0 ? i0 + oc->_M : i0;
408         i3 = i3 >= oc->_M ? i3 - oc->_M : i3;
409
410         j0 = (j1-1);
411         j3 = (j2+1);
412         j0 = j0 <   0 ? j0 + oc->_N : j0;
413         j3 = j3 >= oc->_N ? j3 - oc->_N : j3;
414
415 #define INTERP(m) catrom(catrom(m[i0*oc->_N+j0], m[i1*oc->_N+j0], m[i2*oc->_N+j0], m[i3*oc->_N+j0], frac_x), \
416         catrom(m[i0*oc->_N+j1], m[i1*oc->_N+j1], m[i2*oc->_N+j1], m[i3*oc->_N+j1], frac_x), \
417         catrom(m[i0*oc->_N+j2], m[i1*oc->_N+j2], m[i2*oc->_N+j2], m[i3*oc->_N+j2], frac_x), \
418         catrom(m[i0*oc->_N+j3], m[i1*oc->_N+j3], m[i2*oc->_N+j3], m[i3*oc->_N+j3], frac_x), \
419         frac_z)
420
421         {
422                 if (oc->_do_disp_y) {
423                         ocr->disp[1] = INTERP(oc->_disp_y);
424                 }
425                 if (oc->_do_normals) {
426                         ocr->normal[0] = INTERP(oc->_N_x);
427                         ocr->normal[1] = oc->_N_y/*INTERP(oc->_N_y) (MEM01)*/;
428                         ocr->normal[2] = INTERP(oc->_N_z);
429                 }
430                 if (oc->_do_chop) {
431                         ocr->disp[0] = INTERP(oc->_disp_x);
432                         ocr->disp[2] = INTERP(oc->_disp_z);
433                 }
434                 else {
435                         ocr->disp[0] = 0.0;
436                         ocr->disp[2] = 0.0;
437                 }
438
439                 if (oc->_do_jacobian) {
440                         compute_eigenstuff(ocr, INTERP(oc->_Jxx), INTERP(oc->_Jzz), INTERP(oc->_Jxz));
441                 }
442         }
443 #undef INTERP
444
445         BLI_rw_mutex_unlock(&oc->oceanmutex);
446
447 }
448
449 void BKE_ocean_eval_xz(struct Ocean *oc, struct OceanResult *ocr, float x, float z)
450 {
451         BKE_ocean_eval_uv(oc, ocr, x/oc->_Lx, z/oc->_Lz);
452 }
453
454 void BKE_ocean_eval_xz_catrom(struct Ocean *oc, struct OceanResult *ocr, float x, float z)
455 {
456         BKE_ocean_eval_uv_catrom(oc, ocr, x/oc->_Lx, z/oc->_Lz);
457 }
458
459 // note that this doesn't wrap properly for i, j < 0, but its
460 // not really meant for that being just a way to get the raw data out
461 // to save in some image format.
462 void BKE_ocean_eval_ij(struct Ocean *oc, struct OceanResult *ocr, int i, int j)
463 {
464         BLI_rw_mutex_lock(&oc->oceanmutex, THREAD_LOCK_READ);
465
466         i = abs(i) % oc->_M;
467         j = abs(j) % oc->_N;
468
469         ocr->disp[1] = oc->_do_disp_y ? oc->_disp_y[i*oc->_N+j] : 0.0f;
470
471         if (oc->_do_chop) {
472                 ocr->disp[0] = oc->_disp_x[i*oc->_N+j];
473                 ocr->disp[2] = oc->_disp_z[i*oc->_N+j];
474         }
475         else {
476                 ocr->disp[0] = 0.0f;
477                 ocr->disp[2] = 0.0f;
478         }
479
480         if (oc->_do_normals) {
481                 ocr->normal[0] = oc->_N_x[i*oc->_N+j];
482                 ocr->normal[1] = oc->_N_y/*oc->_N_y[i*oc->_N+j] (MEM01)*/;
483                 ocr->normal[2] = oc->_N_z[i*oc->_N+j];
484
485                 normalize_v3(ocr->normal);
486         }
487
488         if (oc->_do_jacobian) {
489                 compute_eigenstuff(ocr, oc->_Jxx[i*oc->_N+j], oc->_Jzz[i*oc->_N+j], oc->_Jxz[i*oc->_N+j]);
490         }
491
492         BLI_rw_mutex_unlock(&oc->oceanmutex);
493 }
494
495 void BKE_simulate_ocean(struct Ocean *o, float t, float scale, float chop_amount)
496 {
497         int i, j;
498
499         scale *= o->normalize_factor;
500
501         BLI_rw_mutex_lock(&o->oceanmutex, THREAD_LOCK_WRITE);
502
503         // compute a new htilda
504 #pragma omp parallel for private(i, j)
505         for (i = 0 ; i  < o->_M ; ++i) {
506                 // note the <= _N/2 here, see the fftw doco about
507                 // the mechanics of the complex->real fft storage
508                 for ( j  = 0 ; j  <= o->_N / 2 ; ++j) {
509                         fftw_complex exp_param1;
510                         fftw_complex exp_param2;
511                         fftw_complex conj_param;
512
513
514                         init_complex(exp_param1, 0.0, omega(o->_k[i*(1+o->_N/2)+j], o->_depth)*t);
515                         init_complex(exp_param2, 0.0, -omega(o->_k[i*(1+o->_N/2)+j], o->_depth)*t);
516                         exp_complex(exp_param1, exp_param1);
517                         exp_complex(exp_param2, exp_param2);
518                         conj_complex(conj_param, o->_h0_minus[i*o->_N+j]);
519
520                         mul_complex_c(exp_param1, o->_h0[i*o->_N+j], exp_param1);
521                         mul_complex_c(exp_param2, conj_param, exp_param2);
522
523                         add_comlex_c(o->_htilda[i*(1+o->_N/2)+j], exp_param1, exp_param2);
524                         mul_complex_f(o->_fft_in[i*(1+o->_N/2)+j], o->_htilda[i*(1+o->_N/2)+j], scale);
525                 }
526         }
527
528 #pragma omp parallel sections private(i, j)
529         {
530
531 #pragma omp section
532                 {
533                         if (o->_do_disp_y) {
534                                 // y displacement
535                                 fftw_execute(o->_disp_y_plan);
536                         }
537                 } // section 1
538
539 #pragma omp section
540                 {
541                         if (o->_do_chop) {
542                                 // x displacement
543                                 for ( i = 0 ; i  < o->_M ; ++i) {
544                                         for ( j  = 0 ; j  <= o->_N / 2 ; ++j) {
545                                                 fftw_complex mul_param;
546                                                 fftw_complex minus_i;
547
548                                                 init_complex(minus_i, 0.0, -1.0);
549                                                 init_complex(mul_param, -scale, 0);
550                                                 mul_complex_f(mul_param, mul_param, chop_amount);
551                                                 mul_complex_c(mul_param, mul_param, minus_i);
552                                                 mul_complex_c(mul_param, mul_param, o->_htilda[i*(1+o->_N/2)+j]);
553                                                 mul_complex_f(mul_param, mul_param, (o->_k[i*(1+o->_N/2)+j] == 0.0 ? 0.0 : o->_kx[i] / o->_k[i*(1+o->_N/2)+j]));
554                                                 init_complex(o->_fft_in_x[i*(1+o->_N/2)+j], real_c(mul_param), image_c(mul_param));
555                                         }
556                                 }
557                                 fftw_execute(o->_disp_x_plan);
558                         }
559                 } //section 2
560
561 #pragma omp section
562                 {
563                         if (o->_do_chop) {
564                                 // z displacement
565                                 for ( i = 0 ; i  < o->_M ; ++i) {
566                                         for ( j  = 0 ; j  <= o->_N / 2 ; ++j) {
567                                                 fftw_complex mul_param;
568                                                 fftw_complex minus_i;
569
570                                                 init_complex(minus_i, 0.0, -1.0);
571                                                 init_complex(mul_param, -scale, 0);
572                                                 mul_complex_f(mul_param, mul_param, chop_amount);
573                                                 mul_complex_c(mul_param, mul_param, minus_i);
574                                                 mul_complex_c(mul_param, mul_param, o->_htilda[i*(1+o->_N/2)+j]);
575                                                 mul_complex_f(mul_param, mul_param, (o->_k[i*(1+o->_N/2)+j] == 0.0 ? 0.0 : o->_kz[j] / o->_k[i*(1+o->_N/2)+j]));
576                                                 init_complex(o->_fft_in_z[i*(1+o->_N/2)+j], real_c(mul_param), image_c(mul_param));
577                                         }
578                                 }
579                                 fftw_execute(o->_disp_z_plan);
580                         }
581                 } // section 3
582
583 #pragma omp section
584                 {
585                         if (o->_do_jacobian) {
586                                 // Jxx
587                                 for ( i = 0 ; i  < o->_M ; ++i) {
588                                         for ( j  = 0 ; j  <= o->_N / 2 ; ++j) {
589                                                 fftw_complex mul_param;
590
591                                                 //init_complex(mul_param, -scale, 0);
592                                                 init_complex(mul_param, -1, 0);
593
594                                                 mul_complex_f(mul_param, mul_param, chop_amount);
595                                                 mul_complex_c(mul_param, mul_param, o->_htilda[i*(1+o->_N/2)+j]);
596                                                 mul_complex_f(mul_param, mul_param, (o->_k[i*(1+o->_N/2)+j] == 0.0 ? 0.0 : o->_kx[i]*o->_kx[i] / o->_k[i*(1+o->_N/2)+j]));
597                                                 init_complex(o->_fft_in_jxx[i*(1+o->_N/2)+j], real_c(mul_param), image_c(mul_param));
598                                         }
599                                 }
600                                 fftw_execute(o->_Jxx_plan);
601
602                                 for ( i = 0 ; i  < o->_M ; ++i) {
603                                         for ( j  = 0 ; j  < o->_N ; ++j) {
604                                                 o->_Jxx[i*o->_N+j] += 1.0;
605                                         }
606                                 }
607                         }
608                 } // section 4
609
610 #pragma omp section
611                 {
612                         if (o->_do_jacobian) {
613                                 // Jzz
614                                 for ( i = 0 ; i  < o->_M ; ++i) {
615                                         for ( j  = 0 ; j  <= o->_N / 2 ; ++j) {
616                                                 fftw_complex mul_param;
617
618                                                 //init_complex(mul_param, -scale, 0);
619                                                 init_complex(mul_param, -1, 0);
620
621                                                 mul_complex_f(mul_param, mul_param, chop_amount);
622                                                 mul_complex_c(mul_param, mul_param, o->_htilda[i*(1+o->_N/2)+j]);
623                                                 mul_complex_f(mul_param, mul_param, (o->_k[i*(1+o->_N/2)+j] == 0.0 ? 0.0 : o->_kz[j]*o->_kz[j] / o->_k[i*(1+o->_N/2)+j]));
624                                                 init_complex(o->_fft_in_jzz[i*(1+o->_N/2)+j], real_c(mul_param), image_c(mul_param));
625                                         }
626                                 }
627                                 fftw_execute(o->_Jzz_plan);
628                                 for ( i = 0 ; i  < o->_M ; ++i) {
629                                         for ( j  = 0 ; j  < o->_N ; ++j) {
630                                                 o->_Jzz[i*o->_N+j] += 1.0;
631                                         }
632                                 }
633                         }
634                 } // section 5
635
636 #pragma omp section
637                 {
638                         if (o->_do_jacobian) {
639                                 // Jxz
640                                 for ( i = 0 ; i  < o->_M ; ++i) {
641                                         for ( j  = 0 ; j  <= o->_N / 2 ; ++j) {
642                                                 fftw_complex mul_param;
643
644                                                 //init_complex(mul_param, -scale, 0);
645                                                 init_complex(mul_param, -1, 0);
646
647                                                 mul_complex_f(mul_param, mul_param, chop_amount);
648                                                 mul_complex_c(mul_param, mul_param, o->_htilda[i*(1+o->_N/2)+j]);
649                                                 mul_complex_f(mul_param, mul_param, (o->_k[i*(1+o->_N/2)+j] == 0.0f ? 0.0f : o->_kx[i]*o->_kz[j] / o->_k[i*(1+o->_N/2)+j]));
650                                                 init_complex(o->_fft_in_jxz[i*(1+o->_N/2)+j], real_c(mul_param), image_c(mul_param));
651                                         }
652                                 }
653                                 fftw_execute(o->_Jxz_plan);
654                         }
655                 } // section 6
656
657 #pragma omp section
658                 {
659                         // fft normals
660                         if (o->_do_normals) {
661                                 for ( i = 0 ; i  < o->_M ; ++i) {
662                                         for ( j  = 0 ; j  <= o->_N / 2 ; ++j) {
663                                                 fftw_complex mul_param;
664
665                                                 init_complex(mul_param, 0.0, -1.0);
666                                                 mul_complex_c(mul_param, mul_param, o->_htilda[i*(1+o->_N/2)+j]);
667                                                 mul_complex_f(mul_param, mul_param, o->_kx[i]);
668                                                 init_complex(o->_fft_in_nx[i*(1+o->_N/2)+j], real_c(mul_param), image_c(mul_param));
669                                         }
670                                 }
671                                 fftw_execute(o->_N_x_plan);
672
673                         }
674                 } // section 7
675
676 #pragma omp section
677                 {
678                         if (o->_do_normals) {
679                                 for ( i = 0 ; i  < o->_M ; ++i) {
680                                         for ( j  = 0 ; j  <= o->_N / 2 ; ++j) {
681                                                 fftw_complex mul_param;
682
683                                                 init_complex(mul_param, 0.0, -1.0);
684                                                 mul_complex_c(mul_param, mul_param, o->_htilda[i*(1+o->_N/2)+j]);
685                                                 mul_complex_f(mul_param, mul_param, o->_kz[i]);
686                                                 init_complex(o->_fft_in_nz[i*(1+o->_N/2)+j], real_c(mul_param), image_c(mul_param));
687                                         }
688                                 }
689                                 fftw_execute(o->_N_z_plan);
690
691 #if 0
692                                 for ( i = 0 ; i  < o->_M ; ++i) {
693                                         for ( j  = 0 ; j  < o->_N ; ++j) {
694                                                 o->_N_y[i*o->_N+j] = 1.0f/scale;
695                                         }
696                                 }
697                                 (MEM01)
698 #endif
699                         o->_N_y = 1.0f/scale;
700                         }
701                 } // section 8
702
703         } // omp sections
704
705         BLI_rw_mutex_unlock(&o->oceanmutex);
706 }
707
708 static void set_height_normalize_factor(struct Ocean *oc)
709 {
710         float res = 1.0;
711         float max_h = 0.0;
712
713         int i, j;
714
715         if (!oc->_do_disp_y) return;
716
717         oc->normalize_factor = 1.0;
718
719         BKE_simulate_ocean(oc, 0.0, 1.0, 0);
720
721         BLI_rw_mutex_lock(&oc->oceanmutex, THREAD_LOCK_READ);
722
723         for (i = 0; i < oc->_M; ++i) {
724                 for (j = 0; j < oc->_N; ++j) {
725                         if ( max_h < fabsf(oc->_disp_y[i*oc->_N+j])) {
726                                 max_h = fabsf(oc->_disp_y[i*oc->_N+j]);
727                         }
728                 }
729         }
730
731         BLI_rw_mutex_unlock(&oc->oceanmutex);
732
733         if (max_h == 0.0f) max_h = 0.00001f; // just in case ...
734
735         res = 1.0f / (max_h);
736
737         oc->normalize_factor = res;
738 }
739
740 struct Ocean *BKE_add_ocean(void)
741 {
742         Ocean *oc = MEM_callocN(sizeof(Ocean), "ocean sim data");
743
744         BLI_rw_mutex_init(&oc->oceanmutex);
745
746         return oc;
747 }
748
749 void BKE_init_ocean(struct Ocean* o, int M, int N, float Lx, float Lz, float V, float l, float A, float w, float damp,
750                                         float alignment, float depth, float time, short do_height_field, short do_chop, short do_normals, short do_jacobian, int seed)
751 {
752         int i, j, ii;
753
754         BLI_rw_mutex_lock(&o->oceanmutex, THREAD_LOCK_WRITE);
755
756         o->_M = M;
757         o->_N = N;
758         o->_V = V;
759         o->_l = l;
760         o->_A = A;
761         o->_w = w;
762         o->_damp_reflections = 1.0f - damp;
763         o->_wind_alignment = alignment;
764         o->_depth = depth;
765         o->_Lx = Lx;
766         o->_Lz = Lz;
767         o->_wx = cos(w);
768         o->_wz = -sin(w); // wave direction
769         o->_L = V*V / GRAVITY;  // largest wave for a given velocity V
770         o->time = time;
771
772         o->_do_disp_y = do_height_field;
773         o->_do_normals = do_normals;
774         o->_do_chop = do_chop;
775         o->_do_jacobian = do_jacobian;
776
777         o->_k = (float*) MEM_mallocN(M * (1+N/2) * sizeof(float), "ocean_k");
778         o->_h0 = (fftw_complex*) MEM_mallocN(M * N * sizeof(fftw_complex), "ocean_h0");
779         o->_h0_minus = (fftw_complex*) MEM_mallocN(M * N * sizeof(fftw_complex), "ocean_h0_minus");
780         o->_kx = (float*) MEM_mallocN(o->_M * sizeof(float), "ocean_kx");
781         o->_kz = (float*) MEM_mallocN(o->_N * sizeof(float), "ocean_kz");
782
783         // make this robust in the face of erroneous usage
784         if (o->_Lx == 0.0f)
785                 o->_Lx = 0.001f;
786
787         if (o->_Lz == 0.0f)
788                 o->_Lz = 0.001f;
789
790         // the +ve components and DC
791         for (i = 0 ; i <= o->_M/2 ; ++i)
792                 o->_kx[i] = 2.0f * (float)M_PI * i / o->_Lx;
793
794         // the -ve components
795         for (i = o->_M-1, ii=0 ; i > o->_M/2 ; --i, ++ii)
796                 o->_kx[i] = -2.0f * (float)M_PI * ii / o->_Lx;
797
798         // the +ve components and DC
799         for (i = 0 ; i <= o->_N/2 ; ++i)
800                 o->_kz[i] = 2.0f * (float)M_PI * i / o->_Lz;
801
802         // the -ve components
803         for (i = o->_N-1, ii=0 ; i > o->_N/2 ; --i, ++ii)
804                 o->_kz[i] = -2.0f * (float)M_PI * ii / o->_Lz;
805
806         // pre-calculate the k matrix
807         for (i = 0 ; i  < o->_M ; ++i)
808                 for (j  = 0 ; j  <= o->_N / 2 ; ++j)
809                         o->_k[i*(1+o->_N/2)+j] = sqrt(o->_kx[i]*o->_kx[i] + o->_kz[j]*o->_kz[j] );
810
811         /*srand(seed);*/
812         BLI_srand(seed);
813
814         for (i = 0 ; i  < o->_M ; ++i) {
815                 for (j = 0 ; j  < o->_N ; ++j) {
816                         float r1 = gaussRand();
817                         float r2 = gaussRand();
818
819                         fftw_complex r1r2;
820                         init_complex(r1r2, r1, r2);
821                         mul_complex_f(o->_h0[i*o->_N+j], r1r2, (float)(sqrt(Ph(o,  o->_kx[i], o->_kz[j]) / 2.0f)));
822                         mul_complex_f(o->_h0_minus[i*o->_N+j], r1r2, (float)(sqrt(Ph(o, -o->_kx[i], -o->_kz[j]) / 2.0f)));
823                 }
824         }
825
826         o->_fft_in = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in");
827         o->_htilda = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_htilda");
828
829         if (o->_do_disp_y) {
830                 o->_disp_y = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_disp_y");
831                 o->_disp_y_plan = fftw_plan_dft_c2r_2d(o->_M, o->_N, o->_fft_in, o->_disp_y, FFTW_ESTIMATE);
832         }
833
834         if (o->_do_normals) {
835                 o->_fft_in_nx = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in_nx");
836                 o->_fft_in_nz = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in_nz");
837
838                 o->_N_x = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_N_x");
839                 /*o->_N_y = (float*) fftwf_malloc(o->_M * o->_N * sizeof(float)); (MEM01)*/
840                 o->_N_z = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_N_z");
841
842                 o->_N_x_plan = fftw_plan_dft_c2r_2d(o->_M, o->_N, o->_fft_in_nx, o->_N_x, FFTW_ESTIMATE);
843                 o->_N_z_plan = fftw_plan_dft_c2r_2d(o->_M, o->_N, o->_fft_in_nz, o->_N_z, FFTW_ESTIMATE);
844         }
845
846         if (o->_do_chop) {
847                 o->_fft_in_x = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in_x");
848                 o->_fft_in_z = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in_z");
849
850                 o->_disp_x = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_disp_x");
851                 o->_disp_z = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_disp_z");
852
853                 o->_disp_x_plan = fftw_plan_dft_c2r_2d(o->_M, o->_N, o->_fft_in_x, o->_disp_x, FFTW_ESTIMATE);
854                 o->_disp_z_plan = fftw_plan_dft_c2r_2d(o->_M, o->_N, o->_fft_in_z, o->_disp_z, FFTW_ESTIMATE);
855         }
856         if (o->_do_jacobian) {
857                 o->_fft_in_jxx = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in_jxx");
858                 o->_fft_in_jzz = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in_jzz");
859                 o->_fft_in_jxz = (fftw_complex*) MEM_mallocN(o->_M * (1+o->_N/2) * sizeof(fftw_complex), "ocean_fft_in_jxz");
860
861                 o->_Jxx = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_Jxx");
862                 o->_Jzz = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_Jzz");
863                 o->_Jxz = (double*) MEM_mallocN(o->_M * o->_N * sizeof(double), "ocean_Jxz");
864
865                 o->_Jxx_plan = fftw_plan_dft_c2r_2d(o->_M, o->_N, o->_fft_in_jxx, o->_Jxx, FFTW_ESTIMATE);
866                 o->_Jzz_plan = fftw_plan_dft_c2r_2d(o->_M, o->_N, o->_fft_in_jzz, o->_Jzz, FFTW_ESTIMATE);
867                 o->_Jxz_plan = fftw_plan_dft_c2r_2d(o->_M, o->_N, o->_fft_in_jxz, o->_Jxz, FFTW_ESTIMATE);
868         }
869
870         BLI_rw_mutex_unlock(&o->oceanmutex);
871
872         set_height_normalize_factor(o);
873
874 }
875
876 void BKE_free_ocean_data(struct Ocean *oc)
877 {
878         if (!oc) return;
879
880         BLI_rw_mutex_lock(&oc->oceanmutex, THREAD_LOCK_WRITE);
881
882         if (oc->_do_disp_y) {
883                 fftw_destroy_plan(oc->_disp_y_plan);
884                 MEM_freeN(oc->_disp_y);
885         }
886
887         if (oc->_do_normals) {
888                 MEM_freeN(oc->_fft_in_nx);
889                 MEM_freeN(oc->_fft_in_nz);
890                 fftw_destroy_plan(oc->_N_x_plan);
891                 fftw_destroy_plan(oc->_N_z_plan);
892                 MEM_freeN(oc->_N_x);
893                 /*fftwf_free(oc->_N_y); (MEM01)*/
894                 MEM_freeN(oc->_N_z);
895         }
896
897         if (oc->_do_chop) {
898                 MEM_freeN(oc->_fft_in_x);
899                 MEM_freeN(oc->_fft_in_z);
900                 fftw_destroy_plan(oc->_disp_x_plan);
901                 fftw_destroy_plan(oc->_disp_z_plan);
902                 MEM_freeN(oc->_disp_x);
903                 MEM_freeN(oc->_disp_z);
904         }
905
906         if (oc->_do_jacobian) {
907                 MEM_freeN(oc->_fft_in_jxx);
908                 MEM_freeN(oc->_fft_in_jzz);
909                 MEM_freeN(oc->_fft_in_jxz);
910                 fftw_destroy_plan(oc->_Jxx_plan);
911                 fftw_destroy_plan(oc->_Jzz_plan);
912                 fftw_destroy_plan(oc->_Jxz_plan);
913                 MEM_freeN(oc->_Jxx);
914                 MEM_freeN(oc->_Jzz);
915                 MEM_freeN(oc->_Jxz);
916         }
917
918         if (oc->_fft_in)
919                 MEM_freeN(oc->_fft_in);
920
921         /* check that ocean data has been initialized */
922         if (oc->_htilda) {
923                 MEM_freeN(oc->_htilda);
924                 MEM_freeN(oc->_k);
925                 MEM_freeN(oc->_h0);
926                 MEM_freeN(oc->_h0_minus);
927                 MEM_freeN(oc->_kx);
928                 MEM_freeN(oc->_kz);
929         }
930
931         BLI_rw_mutex_unlock(&oc->oceanmutex);
932 }
933
934 void BKE_free_ocean(struct Ocean *oc)
935 {
936         if (!oc) return;
937
938         BKE_free_ocean_data(oc);
939         BLI_rw_mutex_end(&oc->oceanmutex);
940
941         MEM_freeN(oc);
942 }
943
944 #undef GRAVITY
945
946
947 /* ********* Baking/Caching ********* */
948
949
950 #define CACHE_TYPE_DISPLACE     1
951 #define CACHE_TYPE_FOAM         2
952 #define CACHE_TYPE_NORMAL       3
953
954 static void cache_filename(char *string, const char *path, const char *relbase, int frame, int type)
955 {
956         char cachepath[FILE_MAX];
957         const char *fname;
958
959         switch (type) {
960         case CACHE_TYPE_FOAM:
961                 fname= "foam_";
962                 break;
963         case CACHE_TYPE_NORMAL:
964                 fname= "normal_";
965                 break;
966         case CACHE_TYPE_DISPLACE:
967         default:
968                 fname= "disp_";
969                 break;
970         }
971
972         BLI_join_dirfile(cachepath, sizeof(cachepath), path, fname);
973
974         BKE_makepicstring(string, cachepath, relbase, frame, R_IMF_IMTYPE_OPENEXR, 1, TRUE);
975 }
976
977 /* silly functions but useful to inline when the args do a lot of indirections */
978 MINLINE void rgb_to_rgba_unit_alpha(float r_rgba[4], const float rgb[3])
979 {
980         r_rgba[0]= rgb[0];
981         r_rgba[1]= rgb[1];
982         r_rgba[2]= rgb[2];
983         r_rgba[3]= 1.0f;
984 }
985 MINLINE void value_to_rgba_unit_alpha(float r_rgba[4], const float value)
986 {
987         r_rgba[0]= value;
988         r_rgba[1]= value;
989         r_rgba[2]= value;
990         r_rgba[3]= 1.0f;
991 }
992
993 void BKE_free_ocean_cache(struct OceanCache *och)
994 {
995         int i, f=0;
996
997         if (!och) return;
998
999         if (och->ibufs_disp) {
1000                 for (i=och->start, f=0; i<=och->end; i++, f++) {
1001                         if (och->ibufs_disp[f]) {
1002                                 IMB_freeImBuf(och->ibufs_disp[f]);
1003                         }
1004                 }
1005                 MEM_freeN(och->ibufs_disp);
1006         }
1007
1008         if (och->ibufs_foam) {
1009                 for (i=och->start, f=0; i<=och->end; i++, f++) {
1010                         if (och->ibufs_foam[f]) {
1011                                 IMB_freeImBuf(och->ibufs_foam[f]);
1012                         }
1013                 }
1014                 MEM_freeN(och->ibufs_foam);
1015         }
1016
1017         if (och->ibufs_norm) {
1018                 for (i=och->start, f=0; i<=och->end; i++, f++) {
1019                         if (och->ibufs_norm[f]) {
1020                                 IMB_freeImBuf(och->ibufs_norm[f]);
1021                         }
1022                 }
1023                 MEM_freeN(och->ibufs_norm);
1024         }
1025
1026         if (och->time)
1027                 MEM_freeN(och->time);
1028         MEM_freeN(och);
1029 }
1030
1031 void BKE_ocean_cache_eval_uv(struct OceanCache *och, struct OceanResult *ocr, int f, float u, float v)
1032 {
1033         int res_x = och->resolution_x;
1034         int res_y = och->resolution_y;
1035         float result[4];
1036
1037         u = fmod(u, 1.0);
1038         v = fmod(v, 1.0);
1039
1040         if (u < 0) u += 1.0f;
1041         if (v < 0) v += 1.0f;
1042
1043         if (och->ibufs_disp[f]) {
1044                 ibuf_sample(och->ibufs_disp[f], u, v, (1.0f/(float)res_x), (1.0f/(float)res_y), result);
1045                 copy_v3_v3(ocr->disp, result);
1046         }
1047
1048         if (och->ibufs_foam[f]) {
1049                 ibuf_sample(och->ibufs_foam[f], u, v, (1.0f/(float)res_x), (1.0f/(float)res_y), result);
1050                 ocr->foam = result[0];
1051         }
1052
1053         if (och->ibufs_norm[f]) {
1054                 ibuf_sample(och->ibufs_norm[f], u, v, (1.0f/(float)res_x), (1.0f/(float)res_y), result);
1055                 copy_v3_v3(ocr->normal, result);
1056         }
1057 }
1058
1059 void BKE_ocean_cache_eval_ij(struct OceanCache *och, struct OceanResult *ocr, int f, int i, int j)
1060 {
1061         const int res_x = och->resolution_x;
1062         const int res_y = och->resolution_y;
1063
1064         if (i < 0) i= -i;
1065         if (j < 0) j= -j;
1066
1067         i = i % res_x;
1068         j = j % res_y;
1069
1070         if (och->ibufs_disp[f]) {
1071                 copy_v3_v3(ocr->disp, &och->ibufs_disp[f]->rect_float[4*(res_x*j + i)]);
1072         }
1073
1074         if (och->ibufs_foam[f]) {
1075                 ocr->foam = och->ibufs_foam[f]->rect_float[4*(res_x*j + i)];
1076         }
1077
1078         if (och->ibufs_norm[f]) {
1079                 copy_v3_v3(ocr->normal, &och->ibufs_norm[f]->rect_float[4*(res_x*j + i)]);
1080         }
1081 }
1082
1083 struct OceanCache *BKE_init_ocean_cache(const char *bakepath, const char *relbase,
1084                                         int start, int end, float wave_scale,
1085                                         float chop_amount, float foam_coverage, float foam_fade, int resolution)
1086 {
1087         OceanCache *och = MEM_callocN(sizeof(OceanCache), "ocean cache data");
1088
1089         och->bakepath = bakepath;
1090         och->relbase = relbase;
1091
1092         och->start = start;
1093         och->end = end;
1094         och->duration = (end - start) + 1;
1095         och->wave_scale = wave_scale;
1096         och->chop_amount = chop_amount;
1097         och->foam_coverage = foam_coverage;
1098         och->foam_fade = foam_fade;
1099         och->resolution_x = resolution*resolution;
1100         och->resolution_y = resolution*resolution;
1101
1102         och->ibufs_disp = MEM_callocN(sizeof(ImBuf *)*och->duration, "displacement imbuf pointer array");
1103         och->ibufs_foam = MEM_callocN(sizeof(ImBuf *)*och->duration, "foam imbuf pointer array");
1104         och->ibufs_norm = MEM_callocN(sizeof(ImBuf *)*och->duration, "normal imbuf pointer array");
1105
1106         och->time = NULL;
1107
1108         return och;
1109 }
1110
1111 void BKE_simulate_ocean_cache(struct OceanCache *och, int frame)
1112 {
1113         char string[FILE_MAX];
1114         int f = frame;
1115
1116         /* ibufs array is zero based, but filenames are based on frame numbers */
1117         /* still need to clamp frame numbers to valid range of images on disk though */
1118         CLAMP(frame, och->start, och->end);
1119         f = frame - och->start; // shift to 0 based
1120
1121         /* if image is already loaded in mem, return */
1122         if (och->ibufs_disp[f] != NULL ) return;
1123
1124
1125         cache_filename(string, och->bakepath, och->relbase, frame, CACHE_TYPE_DISPLACE);
1126         och->ibufs_disp[f] = IMB_loadiffname(string, 0);
1127         //if (och->ibufs_disp[f] == NULL) printf("error loading %s\n", string);
1128         //else printf("loaded cache %s\n", string);
1129
1130         cache_filename(string, och->bakepath, och->relbase, frame, CACHE_TYPE_FOAM);
1131         och->ibufs_foam[f] = IMB_loadiffname(string, 0);
1132         //if (och->ibufs_foam[f] == NULL) printf("error loading %s\n", string);
1133         //else printf("loaded cache %s\n", string);
1134
1135         cache_filename(string, och->bakepath, och->relbase, frame, CACHE_TYPE_NORMAL);
1136         och->ibufs_norm[f] = IMB_loadiffname(string, 0);
1137         //if (och->ibufs_norm[f] == NULL) printf("error loading %s\n", string);
1138         //else printf("loaded cache %s\n", string);
1139 }
1140
1141
1142 void BKE_bake_ocean(struct Ocean *o, struct OceanCache *och, void (*update_cb)(void *, float progress, int *cancel), void *update_cb_data)
1143 {
1144         /* note: some of these values remain uninitialized unless certain options
1145          * are enabled, take care that BKE_ocean_eval_ij() initializes a member
1146          * before use - campbell */
1147         OceanResult ocr;
1148
1149         ImageFormatData imf= {0};
1150
1151         int f, i=0, x, y, cancel=0;
1152         float progress;
1153
1154         ImBuf *ibuf_foam, *ibuf_disp, *ibuf_normal;
1155         float *prev_foam;
1156         int res_x = och->resolution_x;
1157         int res_y = och->resolution_y;
1158         char string[FILE_MAX];
1159
1160         if (!o) return;
1161
1162         if (o->_do_jacobian) prev_foam = MEM_callocN(res_x*res_y*sizeof(float), "previous frame foam bake data");
1163         else                 prev_foam = NULL;
1164
1165         BLI_srand(0);
1166
1167         /* setup image format */
1168         imf.imtype= R_IMF_IMTYPE_OPENEXR;
1169         imf.depth=  R_IMF_CHAN_DEPTH_16;
1170         imf.exr_codec= R_IMF_EXR_CODEC_ZIP;
1171
1172         for (f=och->start, i=0; f<=och->end; f++, i++) {
1173
1174                 /* create a new imbuf to store image for this frame */
1175                 ibuf_foam = IMB_allocImBuf(res_x, res_y, 32, IB_rectfloat);
1176                 ibuf_disp = IMB_allocImBuf(res_x, res_y, 32, IB_rectfloat);
1177                 ibuf_normal = IMB_allocImBuf(res_x, res_y, 32, IB_rectfloat);
1178
1179                 ibuf_disp->profile = ibuf_foam->profile = ibuf_normal->profile = IB_PROFILE_LINEAR_RGB;
1180
1181                 BKE_simulate_ocean(o, och->time[i], och->wave_scale, och->chop_amount);
1182
1183                 /* add new foam */
1184                 for (y=0; y < res_y; y++) {
1185                         for (x=0; x < res_x; x++) {
1186
1187                                 BKE_ocean_eval_ij(o, &ocr, x, y);
1188
1189                                 /* add to the image */
1190                                 rgb_to_rgba_unit_alpha(&ibuf_disp->rect_float[4*(res_x*y + x)], ocr.disp);
1191
1192                                 if (o->_do_jacobian) {
1193                                         /* TODO, cleanup unused code - campbell */
1194
1195                                         float /*r, */ /* UNUSED */ pr=0.0f, foam_result;
1196                                         float neg_disp, neg_eplus;
1197
1198                                         ocr.foam = BKE_ocean_jminus_to_foam(ocr.Jminus, och->foam_coverage);
1199
1200                                         /* accumulate previous value for this cell */
1201                                         if (i > 0) {
1202                                                 pr = prev_foam[res_x*y + x];
1203                                         }
1204
1205                                         /* r = BLI_frand(); */ /* UNUSED */ // randomly reduce foam
1206
1207                                         //pr = pr * och->foam_fade;             // overall fade
1208
1209                                         // remember ocean coord sys is Y up!
1210                                         // break up the foam where height (Y) is low (wave valley),
1211                                         // and X and Z displacement is greatest
1212
1213 #if 0
1214                                         vec[0] = ocr.disp[0];
1215                                         vec[1] = ocr.disp[2];
1216                                         hor_stretch = len_v2(vec);
1217                                         CLAMP(hor_stretch, 0.0, 1.0);
1218 #endif
1219
1220                                         neg_disp = ocr.disp[1] < 0.0f ? 1.0f+ocr.disp[1] : 1.0f;
1221                                         neg_disp = neg_disp < 0.0f ? 0.0f : neg_disp;
1222
1223                                         /* foam, 'ocr.Eplus' only initialized with do_jacobian */
1224                                         neg_eplus = ocr.Eplus[2] < 0.0f ? 1.0f + ocr.Eplus[2]:1.0f;
1225                                         neg_eplus = neg_eplus<0.0f ? 0.0f : neg_eplus;
1226
1227                                         //if (ocr.disp[1] < 0.0 || r > och->foam_fade)
1228                                         //      pr *= och->foam_fade;
1229
1230
1231                                         //pr = pr * (1.0 - hor_stretch) * ocr.disp[1];
1232                                         //pr = pr * neg_disp * neg_eplus;
1233
1234                                         if (pr < 1.0f) pr *=pr;
1235
1236                                         pr *= och->foam_fade * (0.75f + neg_eplus * 0.25f);
1237
1238
1239                                         foam_result = pr + ocr.foam;
1240
1241                                         prev_foam[res_x*y + x] = foam_result;
1242
1243                                         value_to_rgba_unit_alpha(&ibuf_foam->rect_float[4*(res_x*y + x)], foam_result);
1244                                 }
1245
1246                                 if (o->_do_normals) {
1247                                         rgb_to_rgba_unit_alpha(&ibuf_normal->rect_float[4*(res_x*y + x)], ocr.normal);
1248                                 }
1249                         }
1250                 }
1251
1252                 /* write the images */
1253                 cache_filename(string, och->bakepath, och->relbase, f, CACHE_TYPE_DISPLACE);
1254                 if (0 == BKE_imbuf_write(ibuf_disp, string, &imf))
1255                         printf("Cannot save Displacement File Output to %s\n", string);
1256
1257                 if (o->_do_jacobian) {
1258                         cache_filename(string, och->bakepath, och->relbase,  f, CACHE_TYPE_FOAM);
1259                         if (0 == BKE_imbuf_write(ibuf_foam, string, &imf))
1260                                 printf("Cannot save Foam File Output to %s\n", string);
1261                 }
1262
1263                 if (o->_do_normals) {
1264                         cache_filename(string, och->bakepath,  och->relbase, f, CACHE_TYPE_NORMAL);
1265                         if (0 == BKE_imbuf_write(ibuf_normal, string, &imf))
1266                                 printf("Cannot save Normal File Output to %s\n", string);
1267                 }
1268
1269                 IMB_freeImBuf(ibuf_disp);
1270                 IMB_freeImBuf(ibuf_foam);
1271                 IMB_freeImBuf(ibuf_normal);
1272
1273                 progress = (f - och->start) / (float)och->duration;
1274
1275                 update_cb(update_cb_data, progress, &cancel);
1276
1277                 if (cancel) {
1278                         if (prev_foam) MEM_freeN(prev_foam);
1279                         return;
1280                 }
1281         }
1282
1283         if (prev_foam) MEM_freeN(prev_foam);
1284         och->baked = 1;
1285 }
1286
1287 #else // WITH_OCEANSIM
1288
1289 /* stub */
1290 typedef struct Ocean {
1291         /* need some data here, C does not allow empty struct */
1292         int stub;
1293 } Ocean;
1294
1295
1296 float BKE_ocean_jminus_to_foam(float UNUSED(jminus), float UNUSED(coverage))
1297 {
1298         return 0.0f;
1299 }
1300
1301 void BKE_ocean_eval_uv(struct Ocean *UNUSED(oc), struct OceanResult *UNUSED(ocr), float UNUSED(u), float UNUSED(v))
1302 {
1303 }
1304
1305 // use catmullrom interpolation rather than linear
1306 void BKE_ocean_eval_uv_catrom(struct Ocean *UNUSED(oc), struct OceanResult *UNUSED(ocr), float UNUSED(u), float UNUSED(v))
1307 {
1308 }
1309
1310 void BKE_ocean_eval_xz(struct Ocean *UNUSED(oc), struct OceanResult *UNUSED(ocr), float UNUSED(x), float UNUSED(z))
1311 {
1312 }
1313
1314 void BKE_ocean_eval_xz_catrom(struct Ocean *UNUSED(oc), struct OceanResult *UNUSED(ocr), float UNUSED(x), float UNUSED(z))
1315 {
1316 }
1317
1318 void BKE_ocean_eval_ij(struct Ocean *UNUSED(oc), struct OceanResult *UNUSED(ocr), int UNUSED(i), int UNUSED(j))
1319 {
1320 }
1321
1322 void BKE_simulate_ocean(struct Ocean *UNUSED(o), float UNUSED(t), float UNUSED(scale), float UNUSED(chop_amount))
1323 {
1324 }
1325
1326 struct Ocean *BKE_add_ocean(void)
1327 {
1328         Ocean *oc = MEM_callocN(sizeof(Ocean), "ocean sim data");
1329
1330         return oc;
1331 }
1332
1333 void BKE_init_ocean(struct Ocean* UNUSED(o), int UNUSED(M), int UNUSED(N), float UNUSED(Lx), float UNUSED(Lz), float UNUSED(V), float UNUSED(l), float UNUSED(A), float UNUSED(w), float UNUSED(damp),
1334                                            float UNUSED(alignment), float UNUSED(depth), float UNUSED(time), short UNUSED(do_height_field), short UNUSED(do_chop), short UNUSED(do_normals), short UNUSED(do_jacobian), int UNUSED(seed))
1335 {
1336 }
1337
1338 void BKE_free_ocean_data(struct Ocean *UNUSED(oc))
1339 {
1340 }
1341
1342 void BKE_free_ocean(struct Ocean *oc)
1343 {
1344         if (!oc) return;
1345         MEM_freeN(oc);
1346 }
1347
1348
1349 /* ********* Baking/Caching ********* */
1350
1351
1352 void BKE_free_ocean_cache(struct OceanCache *och)
1353 {
1354         if (!och) return;
1355
1356         MEM_freeN(och);
1357 }
1358
1359 void BKE_ocean_cache_eval_uv(struct OceanCache *UNUSED(och), struct OceanResult *UNUSED(ocr), int UNUSED(f), float UNUSED(u), float UNUSED(v))
1360 {
1361 }
1362
1363 void BKE_ocean_cache_eval_ij(struct OceanCache *UNUSED(och), struct OceanResult *UNUSED(ocr), int UNUSED(f), int UNUSED(i), int UNUSED(j))
1364 {
1365 }
1366
1367 struct OceanCache *BKE_init_ocean_cache(const char *UNUSED(bakepath), const char *UNUSED(relbase),
1368                                         int UNUSED(start), int UNUSED(end), float UNUSED(wave_scale),
1369                                         float UNUSED(chop_amount), float UNUSED(foam_coverage), float UNUSED(foam_fade), int UNUSED(resolution))
1370 {
1371         OceanCache *och = MEM_callocN(sizeof(OceanCache), "ocean cache data");
1372
1373         return och;
1374 }
1375
1376 void BKE_simulate_ocean_cache(struct OceanCache *UNUSED(och), int UNUSED(frame))
1377 {
1378 }
1379
1380 void BKE_bake_ocean(struct Ocean *UNUSED(o), struct OceanCache *UNUSED(och), void (*update_cb)(void *, float progress, int *cancel), void *UNUSED(update_cb_data))
1381 {
1382         /* unused */
1383         (void)update_cb;
1384 }
1385 #endif // WITH_OCEANSIM