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