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