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