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