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