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