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