diff options
Diffstat (limited to 'source/blender/blenkernel')
-rw-r--r-- | source/blender/blenkernel/BKE_ocean.h | 8 | ||||
-rw-r--r-- | source/blender/blenkernel/CMakeLists.txt | 1 | ||||
-rw-r--r-- | source/blender/blenkernel/intern/ocean.c | 159 | ||||
-rw-r--r-- | source/blender/blenkernel/intern/ocean_intern.h | 137 | ||||
-rw-r--r-- | source/blender/blenkernel/intern/ocean_spectrum.c | 224 |
5 files changed, 428 insertions, 101 deletions
diff --git a/source/blender/blenkernel/BKE_ocean.h b/source/blender/blenkernel/BKE_ocean.h index 3ba6486e6c5..d3ac825039d 100644 --- a/source/blender/blenkernel/BKE_ocean.h +++ b/source/blender/blenkernel/BKE_ocean.h @@ -84,6 +84,9 @@ void BKE_ocean_init(struct Ocean *o, float alignment, float depth, float time, + int spectrum, + float fetch_jonswap, + float sharpen_peak_jonswap, short do_height_field, short do_chop, short do_normals, @@ -122,6 +125,11 @@ void BKE_ocean_cache_eval_ij(struct OceanCache *och, struct OceanResult *ocr, in void BKE_ocean_free_cache(struct OceanCache *och); void BKE_ocean_free_modifier_cache(struct OceanModifierData *omd); +/* ocean_spectrum.c */ +float BLI_ocean_spectrum_piersonmoskowitz(const struct Ocean *oc, const float kx, const float kz); +float BLI_ocean_spectrum_texelmarsenarsloe(const struct Ocean *oc, const float kx, const float kz); +float BLI_ocean_spectrum_jonswap(const struct Ocean *oc, const float kx, const float kz); + #ifdef __cplusplus } #endif diff --git a/source/blender/blenkernel/CMakeLists.txt b/source/blender/blenkernel/CMakeLists.txt index 58f05f65bef..112933e40be 100644 --- a/source/blender/blenkernel/CMakeLists.txt +++ b/source/blender/blenkernel/CMakeLists.txt @@ -178,6 +178,7 @@ set(SRC intern/object_facemap.c intern/object_update.c intern/ocean.c + intern/ocean_spectrum.c intern/outliner_treehash.c intern/packedFile.c intern/paint.c diff --git a/source/blender/blenkernel/intern/ocean.c b/source/blender/blenkernel/intern/ocean.c index c27fb59835f..26485d10fbd 100644 --- a/source/blender/blenkernel/intern/ocean.c +++ b/source/blender/blenkernel/intern/ocean.c @@ -38,11 +38,11 @@ #include "BLI_path_util.h" #include "BLI_rand.h" #include "BLI_task.h" -#include "BLI_threads.h" #include "BLI_utildefines.h" #include "BKE_image.h" #include "BKE_ocean.h" +#include "ocean_intern.h" #include "IMB_imbuf.h" #include "IMB_imbuf_types.h" @@ -54,95 +54,6 @@ #ifdef WITH_OCEANSIM /* Ocean code */ -# include "fftw3.h" - -# define GRAVITY 9.81f - -typedef struct Ocean { - /* ********* input parameters to the sim ********* */ - float _V; - float _l; - float _w; - float _A; - float _damp_reflections; - float _wind_alignment; - float _depth; - - float _wx; - float _wz; - - float _L; - - /* dimensions of computational grid */ - int _M; - int _N; - - /* spatial size of computational grid */ - float _Lx; - float _Lz; - - float normalize_factor; /* init w */ - float time; - - short _do_disp_y; - short _do_normals; - short _do_chop; - short _do_jacobian; - - /* mutex for threaded texture access */ - ThreadRWMutex oceanmutex; - - /* ********* sim data arrays ********* */ - - /* two dimensional arrays of complex */ - fftw_complex *_fft_in; /* init w sim w */ - fftw_complex *_fft_in_x; /* init w sim w */ - fftw_complex *_fft_in_z; /* init w sim w */ - fftw_complex *_fft_in_jxx; /* init w sim w */ - fftw_complex *_fft_in_jzz; /* init w sim w */ - fftw_complex *_fft_in_jxz; /* init w sim w */ - fftw_complex *_fft_in_nx; /* init w sim w */ - fftw_complex *_fft_in_nz; /* init w sim w */ - fftw_complex *_htilda; /* init w sim w (only once) */ - - /* fftw "plans" */ - fftw_plan _disp_y_plan; /* init w sim r */ - fftw_plan _disp_x_plan; /* init w sim r */ - fftw_plan _disp_z_plan; /* init w sim r */ - fftw_plan _N_x_plan; /* init w sim r */ - fftw_plan _N_z_plan; /* init w sim r */ - fftw_plan _Jxx_plan; /* init w sim r */ - fftw_plan _Jxz_plan; /* init w sim r */ - fftw_plan _Jzz_plan; /* init w sim r */ - - /* two dimensional arrays of float */ - double *_disp_y; /* init w sim w via plan? */ - double *_N_x; /* init w sim w via plan? */ - /* all member of this array has same values, - * so convert this array to a float to reduce memory usage (MEM01). */ - /*float * _N_y; */ - double _N_y; /* sim w ********* can be rearranged? */ - double *_N_z; /* init w sim w via plan? */ - double *_disp_x; /* init w sim w via plan? */ - double *_disp_z; /* init w sim w via plan? */ - - /* two dimensional arrays of float */ - /* Jacobian and minimum eigenvalue */ - double *_Jxx; /* init w sim w */ - double *_Jzz; /* init w sim w */ - double *_Jxz; /* init w sim w */ - - /* one dimensional float array */ - float *_kx; /* init w sim r */ - float *_kz; /* init w sim r */ - - /* two dimensional complex array */ - fftw_complex *_h0; /* init w sim r */ - fftw_complex *_h0_minus; /* init w sim r */ - - /* two dimensional float array */ - float *_k; /* init w sim r */ -} Ocean; static float nextfr(RNG *rng, float min, float max) { @@ -285,7 +196,7 @@ float BKE_ocean_jminus_to_foam(float jminus, float coverage) { float foam = jminus * -0.005f + coverage; CLAMP(foam, 0.0f, 1.0f); - return foam * foam; + return foam; } void BKE_ocean_eval_uv(struct Ocean *oc, struct OceanResult *ocr, float u, float v) @@ -893,6 +804,9 @@ void BKE_ocean_init_from_modifier(struct Ocean *ocean, struct OceanModifierData omd->wave_alignment, omd->depth, omd->time, + omd->spectrum, + omd->fetch_jonswap, + omd->sharpen_peak_jonswap, do_heightfield, do_chop, do_normals, @@ -913,6 +827,9 @@ void BKE_ocean_init(struct Ocean *o, float alignment, float depth, float time, + int spectrum, + float fetch_jonswap, + float sharpen_peak_jonswap, short do_height_field, short do_chop, short do_normals, @@ -940,6 +857,13 @@ void BKE_ocean_init(struct Ocean *o, o->_L = V * V / GRAVITY; /* largest wave for a given velocity V */ o->time = time; + /* Spectrum to use. */ + o->_spectrum = spectrum; + + /* Common JONSWAP parameters. */ + o->_fetch_jonswap = fetch_jonswap; + o->_sharpen_peak_jonswap = sharpen_peak_jonswap; + o->_do_disp_y = do_height_field; o->_do_normals = do_normals; o->_do_chop = do_chop; @@ -1001,10 +925,46 @@ void BKE_ocean_init(struct Ocean *o, fftw_complex r1r2; init_complex(r1r2, r1, r2); - mul_complex_f( - o->_h0[i * o->_N + j], r1r2, (float)(sqrt(Ph(o, o->_kx[i], o->_kz[j]) / 2.0f))); - mul_complex_f( - o->_h0_minus[i * o->_N + j], r1r2, (float)(sqrt(Ph(o, -o->_kx[i], -o->_kz[j]) / 2.0f))); + switch (o->_spectrum) { + case MOD_OCEAN_SPECTRUM_JONSWAP: + mul_complex_f(o->_h0[i * o->_N + j], + r1r2, + (float)(sqrt(BLI_ocean_spectrum_jonswap(o, o->_kx[i], o->_kz[j]) / 2.0f))); + mul_complex_f( + o->_h0_minus[i * o->_N + j], + r1r2, + (float)(sqrt(BLI_ocean_spectrum_jonswap(o, -o->_kx[i], -o->_kz[j]) / 2.0f))); + break; + case MOD_OCEAN_SPECTRUM_TEXEL_MARSEN_ARSLOE: + mul_complex_f( + o->_h0[i * o->_N + j], + r1r2, + (float)(sqrt(BLI_ocean_spectrum_texelmarsenarsloe(o, o->_kx[i], o->_kz[j]) / 2.0f))); + mul_complex_f( + o->_h0_minus[i * o->_N + j], + r1r2, + (float)(sqrt(BLI_ocean_spectrum_texelmarsenarsloe(o, -o->_kx[i], -o->_kz[j]) / + 2.0f))); + break; + case MOD_OCEAN_SPECTRUM_PIERSON_MOSKOWITZ: + mul_complex_f( + o->_h0[i * o->_N + j], + r1r2, + (float)(sqrt(BLI_ocean_spectrum_piersonmoskowitz(o, o->_kx[i], o->_kz[j]) / 2.0f))); + mul_complex_f( + o->_h0_minus[i * o->_N + j], + r1r2, + (float)(sqrt(BLI_ocean_spectrum_piersonmoskowitz(o, -o->_kx[i], -o->_kz[j]) / + 2.0f))); + break; + default: + mul_complex_f( + o->_h0[i * o->_N + j], r1r2, (float)(sqrt(Ph(o, o->_kx[i], o->_kz[j]) / 2.0f))); + mul_complex_f(o->_h0_minus[i * o->_N + j], + r1r2, + (float)(sqrt(Ph(o, -o->_kx[i], -o->_kz[j]) / 2.0f))); + break; + } } } @@ -1517,12 +1477,6 @@ void BKE_ocean_bake(struct Ocean *o, #else /* WITH_OCEANSIM */ -/* stub */ -typedef struct Ocean { - /* need some data here, C does not allow empty struct */ - int stub; -} Ocean; - float BKE_ocean_jminus_to_foam(float UNUSED(jminus), float UNUSED(coverage)) { return 0.0f; @@ -1591,6 +1545,9 @@ void BKE_ocean_init(struct Ocean *UNUSED(o), float UNUSED(alignment), float UNUSED(depth), float UNUSED(time), + int UNUSED(spectrum), + float UNUSED(fetch_jonswap), + float UNUSED(sharpen_peak_jonswap), short UNUSED(do_height_field), short UNUSED(do_chop), short UNUSED(do_normals), diff --git a/source/blender/blenkernel/intern/ocean_intern.h b/source/blender/blenkernel/intern/ocean_intern.h new file mode 100644 index 00000000000..7da88419219 --- /dev/null +++ b/source/blender/blenkernel/intern/ocean_intern.h @@ -0,0 +1,137 @@ +/* + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#ifndef __BKE_OCEAN_INTERN_H__ +#define __BKE_OCEAN_INTERN_H__ + +/** \file + * \ingroup bli + */ + +#ifdef __cplusplus +extern "C" { +#endif + +#ifdef WITH_OCEANSIM +# include "BLI_threads.h" +# include "fftw3.h" +# define GRAVITY 9.81f + +typedef struct Ocean { + /* ********* input parameters to the sim ********* */ + float _V; + float _l; + float _w; + float _A; + float _damp_reflections; + float _wind_alignment; + float _depth; + + float _wx; + float _wz; + + float _L; + + /* dimensions of computational grid */ + int _M; + int _N; + + /* spatial size of computational grid */ + float _Lx; + float _Lz; + + float normalize_factor; /* init w */ + float time; + + short _do_disp_y; + short _do_normals; + short _do_chop; + short _do_jacobian; + + /* Which spectral model we are using. */ + int _spectrum; + + /* JONSWAP common parameters. */ + float _fetch_jonswap; + float _sharpen_peak_jonswap; + + /* mutex for threaded texture access */ + ThreadRWMutex oceanmutex; + + /* ********* sim data arrays ********* */ + + /* two dimensional arrays of complex */ + fftw_complex *_fft_in; /* init w sim w */ + fftw_complex *_fft_in_x; /* init w sim w */ + fftw_complex *_fft_in_z; /* init w sim w */ + fftw_complex *_fft_in_jxx; /* init w sim w */ + fftw_complex *_fft_in_jzz; /* init w sim w */ + fftw_complex *_fft_in_jxz; /* init w sim w */ + fftw_complex *_fft_in_nx; /* init w sim w */ + fftw_complex *_fft_in_nz; /* init w sim w */ + fftw_complex *_htilda; /* init w sim w (only once) */ + + /* fftw "plans" */ + fftw_plan _disp_y_plan; /* init w sim r */ + fftw_plan _disp_x_plan; /* init w sim r */ + fftw_plan _disp_z_plan; /* init w sim r */ + fftw_plan _N_x_plan; /* init w sim r */ + fftw_plan _N_z_plan; /* init w sim r */ + fftw_plan _Jxx_plan; /* init w sim r */ + fftw_plan _Jxz_plan; /* init w sim r */ + fftw_plan _Jzz_plan; /* init w sim r */ + + /* two dimensional arrays of float */ + double *_disp_y; /* init w sim w via plan? */ + double *_N_x; /* init w sim w via plan? */ + /* all member of this array has same values, + * so convert this array to a float to reduce memory usage (MEM01). */ + /*float * _N_y; */ + double _N_y; /* sim w ********* can be rearranged? */ + double *_N_z; /* init w sim w via plan? */ + double *_disp_x; /* init w sim w via plan? */ + double *_disp_z; /* init w sim w via plan? */ + + /* two dimensional arrays of float */ + /* Jacobian and minimum eigenvalue */ + double *_Jxx; /* init w sim w */ + double *_Jzz; /* init w sim w */ + double *_Jxz; /* init w sim w */ + + /* one dimensional float array */ + float *_kx; /* init w sim r */ + float *_kz; /* init w sim r */ + + /* two dimensional complex array */ + fftw_complex *_h0; /* init w sim r */ + fftw_complex *_h0_minus; /* init w sim r */ + + /* two dimensional float array */ + float *_k; /* init w sim r */ +} Ocean; +#else +/* stub */ +typedef struct Ocean { + /* need some data here, C does not allow empty struct */ + int stub; +} Ocean; +#endif + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/source/blender/blenkernel/intern/ocean_spectrum.c b/source/blender/blenkernel/intern/ocean_spectrum.c new file mode 100644 index 00000000000..93018f4a5b0 --- /dev/null +++ b/source/blender/blenkernel/intern/ocean_spectrum.c @@ -0,0 +1,224 @@ +/* + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/** \file + * \ingroup bke + */ + +#include "BLI_math.h" +#include "BKE_ocean.h" +#include "ocean_intern.h" + +#ifdef WITH_OCEANSIM + +/* -------------------------------------------------------------------- */ +/** \name Ocean Spectrum from EncinoWaves + * \{ */ + +/* + * Original code from EncinoWaves project Copyright (c) 2015 Christopher Jon Horvath + * Modifications made to work within blender. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + */ + +/** + * alpha_beta_spectrum is a common algorithm for the Pierson-Moskowitz, JONSWAP and TMA models. + * This is a modified implementation from the EncinoWaves project. + */ +static float alpha_beta_spectrum(const float alpha, + const float beta, + const float gamma, + const float omega, + const float peakomega) +{ + return (alpha * sqrt(gamma) / pow(omega, 5.0)) * exp(-beta * pow(peakomega / omega, 4.0)); +} + +static float peak_sharpen(const float omega, const float m_peakomega, const float m_gamma) +{ + const float peak_sharpening_sigma = (omega < m_peakomega) ? 0.07 : 0.09; + const float peak_sharpening = pow( + m_gamma, exp(-sqrt((omega - m_peakomega) / (peak_sharpening_sigma * m_peakomega)) / 2.0)); + + return peak_sharpening; +} + +/** + * Spectrum-type independent modifications. + */ +static float ocean_spectrum_wind_and_damp(const Ocean *oc, + const float kx, + const float kz, + const float val) +{ + const float k2 = kx * kx + kz * kz; + const float k_mag_inv = 1.0f / k2; + const float k_dot_w = (kx * k_mag_inv * oc->_wx) + (kz * k_mag_inv * oc->_wz); + + /* Bias towards wind dir. */ + float newval = val * pow(fabs(k_dot_w), oc->_wind_alignment); + + /* Eliminate wavelengths smaller than cutoff. */ + /* val *= exp(-k2 * m_cutoff); */ + + /* Reduce reflected waves. */ + if (k_dot_w < 0.0f) { + if (oc->_wind_alignment > 0.0) { + newval *= oc->_damp_reflections; + } + } + + return newval; +} + +static float jonswap(const Ocean *oc, const float k2) +{ + /* Get our basic JONSWAP value from #alpha_beta_spectrum. */ + const float k_mag = sqrt(k2); + + const float m_omega = GRAVITY * k_mag * tanh(k_mag * oc->_depth); + const float omega = sqrt(m_omega); + + const float m_fetch = oc->_fetch_jonswap; + + /* Strictly, this should be a random value from a Gaussian (mean 3.3, variance 0.67), + * clamped 1.0 to 6.0. */ + float m_gamma = oc->_sharpen_peak_jonswap; + if (m_gamma < 1.0) { + m_gamma = 1.00; + } + if (m_gamma > 6.0) { + m_gamma = 6.0; + } + + const float m_windspeed = oc->_V; + + const float m_dimensionlessFetch = fabs(GRAVITY * m_fetch / sqrt(m_windspeed)); + const float m_alpha = 0.076 * pow(m_dimensionlessFetch, -0.22); + + const float m_tau = M_PI * 2; + const float m_peakomega = m_tau * 3.5 * fabs(GRAVITY / oc->_V) * + pow(m_dimensionlessFetch, -0.33); + + const float beta = 1.25f; + + float val = alpha_beta_spectrum(m_alpha, beta, GRAVITY, omega, m_peakomega); + + /* Peak sharpening */ + val *= peak_sharpen(m_omega, m_peakomega, m_gamma); + + return val; +} + +/** + * Pierson-Moskowitz model, 1964, assumes waves reach equilibrium with wind. + * Model is intended for large area 'fully developed' sea, where winds have been steadily blowing + * for days over an area that includes hundreds of wavelengths on a side. + */ +float BLI_ocean_spectrum_piersonmoskowitz(const Ocean *oc, const float kx, const float kz) +{ + const float k2 = kx * kx + kz * kz; + + if (k2 == 0.0f) { + /* No DC component. */ + return 0.0f; + } + + /* Get Pierson-Moskowitz value from #alpha_beta_spectrum. */ + const float peak_omega_PM = 0.87f * GRAVITY / oc->_V; + + const float k_mag = sqrt(k2); + const float m_omega = GRAVITY * k_mag * tanh(k_mag * oc->_depth); + + const float omega = sqrt(m_omega); + const float alpha = 0.0081f; + const float beta = 1.291f; + + float val = alpha_beta_spectrum(alpha, beta, GRAVITY, omega, peak_omega_PM); + + val = ocean_spectrum_wind_and_damp(oc, kx, kz, val); + + return val; +} + +/** + * TMA extends the JONSWAP spectrum. + * This spectral model is best suited to shallow water. + */ +float BLI_ocean_spectrum_texelmarsenarsloe(const Ocean *oc, const float kx, const float kz) +{ + const float k2 = kx * kx + kz * kz; + + if (k2 == 0.0f) { + /* No DC component. */ + return 0.0f; + } + + float val = jonswap(oc, k2); + + val = ocean_spectrum_wind_and_damp(oc, kx, kz, val); + + /* TMA modifications to JONSWAP. */ + const float m_depth = oc->_depth; + const float gain = sqrt(m_depth / GRAVITY); + + const float k_mag = sqrt(k2); + + const float m_omega = GRAVITY * k_mag * tanh(k_mag * oc->_depth); + const float omega = sqrt(m_omega); + + const float kitaigorodskiiDepth_wh = omega * gain; + const float kitaigorodskiiDepth = 0.5 + (0.5 * tanh(1.8 * (kitaigorodskiiDepth_wh - 1.125))); + + val *= kitaigorodskiiDepth; + + val = ocean_spectrum_wind_and_damp(oc, kx, kz, val); + + return val; +} + +/** + * Hasselmann et al, 1973. This model extends the Pierson-Moskowitz model with a peak sharpening + * function This enhancement is an artificial construct to address the problem that the wave + * spectrum is never fully developed. + * + * The fetch parameter represents the distance from a lee shore, + * called the fetch, or the distance over which the wind blows with constant velocity. + */ +float BLI_ocean_spectrum_jonswap(const Ocean *oc, const float kx, const float kz) +{ + const float k2 = kx * kx + kz * kz; + + if (k2 == 0.0f) { + /* No DC component. */ + return 0.0f; + } + + float val = jonswap(oc, k2); + + val = ocean_spectrum_wind_and_damp(oc, kx, kz, val); + + return val; +} + +/** \} */ + +#endif /* WITH_OCEANSIM */ |