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