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