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