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