summaryrefslogtreecommitdiff
path: root/src/ocean
diff options
context:
space:
mode:
Diffstat (limited to 'src/ocean')
-rw-r--r--src/ocean/halftone.cpp47
-rw-r--r--src/ocean/halftone.h19
-rw-r--r--src/ocean/ocean_alt.cpp269
-rw-r--r--src/ocean/ocean_alt.h53
4 files changed, 311 insertions, 77 deletions
diff --git a/src/ocean/halftone.cpp b/src/ocean/halftone.cpp
new file mode 100644
index 0000000..6aa59cb
--- /dev/null
+++ b/src/ocean/halftone.cpp
@@ -0,0 +1,47 @@
+#include "halftone.h"
+#include <cmath>
+#include <vector>
+
+halftone::halftone()
+{
+
+}
+
+//void halftone::init(int n, int m, float density){
+// N = n;
+// M = m;
+// // std::vector<std::vector<int>> m_halftone;
+// // 1. intiialize pattern array with density
+// for (int i=0; i<N; i++){
+// for (int j=0; j<M; j++){
+// // determine whether to place a 1 or 0 here based on density
+
+// // generate random number 0-1
+// float r = ((float) rand() / (float) RAND_MAX);
+// if (r < density){
+// m_halftone[i][j] = 1;
+// } else {
+// m_halftone[i][j] = 0;
+// }
+// }
+// }
+
+//}
+
+//void halftone::apply_gaussian(){
+// for (int i=0; i<N; i++){
+// for (int j=0; j<M; j++){
+// m_halftone[i][j] = radial_gaussian(i, j);
+// }
+// }
+
+//}
+
+//float halftone::radial_gaussian(int i, int j){
+// float r = sqrt((i-.5*N)*(i-.5*N) + (j-.5*M)*(j-.5*M));
+// float result = exp(-(r*r) / (2.f*m_sigma*m_sigma));
+
+// return result;
+
+//}
+
diff --git a/src/ocean/halftone.h b/src/ocean/halftone.h
new file mode 100644
index 0000000..e913eac
--- /dev/null
+++ b/src/ocean/halftone.h
@@ -0,0 +1,19 @@
+#ifndef HALFTONE_H
+#define HALFTONE_H
+
+
+#include <vector>
+class halftone
+{
+public:
+ halftone();
+
+private:
+ int N = 0;
+ int M = 0;
+ float m_sigma = 24.f;
+ std::vector<std::vector<int>> m_halftone;
+
+};
+
+#endif // HALFTONE_H
diff --git a/src/ocean/ocean_alt.cpp b/src/ocean/ocean_alt.cpp
index ca080cd..0e44d9c 100644
--- a/src/ocean/ocean_alt.cpp
+++ b/src/ocean/ocean_alt.cpp
@@ -12,6 +12,7 @@ ocean_alt::ocean_alt()
// initializes static constants (aka they are not time dependent)
void ocean_alt::init_wave_index_constants(){
+ float tex_step = 1.f/num_rows;
for (int i=0; i<N; i++){
Eigen::Vector2i m_n = index_1d_to_2d(i);
@@ -22,6 +23,23 @@ void ocean_alt::init_wave_index_constants(){
Eigen::Vector2d k_conj = get_k_vector(-n_prime, m_prime);
+// Eigen::Vector3f v = Eigen::Vector3f(0,0,1);
+// Eigen::Vector3f norm = Eigen::Vector3f(0,1,0);
+// if (abs(norm[1]) < 1.f){
+// v = (Eigen::Vector3f(0,1,0) - norm[1]*norm).normalized();
+// }
+// Eigen::Vector3f u = norm.cross(v).normalized();
+
+// float u_coord = u.dot(Eigen::Vector3f(n_prime, 0, m_prime)) / 64.f;
+// float v_coord = v.dot(Eigen::Vector3f(n_prime, 0, m_prime)) / 64.f;
+
+// //std::cout << u_coord << ", " << v_coord << std::endl;
+
+
+ // texture coord:
+ Eigen::Vector2f texCoord = Eigen::Vector2f(1, 1);
+
+
// store h0'(n,m) and w'(n,m) for every index, to be used for later
Eigen::Vector2d h0_prime = h_0_prime(k);
@@ -36,7 +54,6 @@ void ocean_alt::init_wave_index_constants(){
wave_const.h0_prime = h0_prime;
wave_const.h0_prime_conj = h0_prime_conj;
wave_const.w_prime = w_prime;
- wave_const.w_prime = w_prime;
wave_const.base_horiz_pos = get_horiz_pos(i);
wave_const.k_vector = k;
@@ -45,10 +62,24 @@ void ocean_alt::init_wave_index_constants(){
// initialize m_current_h to be h0 for now
// m_current_h.push_back(h0_prime);
m_current_h.push_back(Eigen::Vector2d(0.0, 0.0));
- m_displacements.push_back(Eigen::Vector2d(0.0, 0.0));
- m_slopes.push_back(Eigen::Vector2d(0.0, 0.0));
+ // m_displacements.push_back(Eigen::Vector2d(0.0, 0.0));
+ // m_slopes.push_back(Eigen::Vector2d(0.0, 0.0));
m_normals.push_back(Eigen::Vector3f(0.0, 1.0, 0.0));
+ // initialize foam constant vectors
+ m_foam_constants.k_vectors.push_back(Eigen::Vector2f(k[0], k[1]));
+ m_foam_constants.positions.push_back(Eigen::Vector3f(0,0,0));
+ m_foam_constants.wavelengths.push_back(0);
+ // m_foam_constants.texCoords.push_back(texCoord);
+
+
+ m_slopes_x.push_back(Eigen::Vector2d(0.0, 0.0));
+ m_slopes_z.push_back(Eigen::Vector2d(0.0, 0.0));
+
+ m_displacements_x.push_back(Eigen::Vector2d(0.0, 0.0));
+ m_displacements_z.push_back(Eigen::Vector2d(0.0, 0.0));
+
+ m_vertices.push_back(Eigen::Vector3f(0.0, 0.0, 0.0));
}
}
@@ -57,8 +88,14 @@ void ocean_alt::fft_prime(double t){
// FFT
std::vector<Eigen::Vector2d> h_tildas = std::vector<Eigen::Vector2d>();
- std::vector<Eigen::Vector2d> ikh = std::vector<Eigen::Vector2d>();
- std::vector<Eigen::Vector2d> neg_ik_hat_h = std::vector<Eigen::Vector2d>();
+ // std::vector<Eigen::Vector2d> ikh = std::vector<Eigen::Vector2d>();
+// std::vector<Eigen::Vector2d> neg_ik_hat_h = std::vector<Eigen::Vector2d>();
+
+ std::vector<Eigen::Vector2d> ikhx = std::vector<Eigen::Vector2d>();
+ std::vector<Eigen::Vector2d> ikhz = std::vector<Eigen::Vector2d>();
+
+ std::vector<Eigen::Vector2d> neg_ik_hat_h_x = std::vector<Eigen::Vector2d>();
+ std::vector<Eigen::Vector2d> neg_ik_hat_h_z = std::vector<Eigen::Vector2d>();
// find each h_tilda at each index, to be used for next for loop
for (int i=0; i<N; i++){
@@ -67,25 +104,47 @@ void ocean_alt::fft_prime(double t){
h_tildas.emplace_back(h_t_prime);
Eigen::Vector2d k_vector = m_waveIndexConstants[i].k_vector;
- ikh.emplace_back(-h_t_prime[1] * k_vector[0], -h_t_prime[1] * k_vector[1]);
+ // ikh.emplace_back(-h_t_prime[1] * k_vector[0], -h_t_prime[1] * k_vector[1]);
+ ikhx.emplace_back(-h_t_prime[1] * k_vector[0], -h_t_prime[1] * k_vector[0]);
+ ikhz.emplace_back(-h_t_prime[1] * k_vector[1], -h_t_prime[1] * k_vector[1]);
+
- Eigen::Vector2d k_normalized = k_vector.normalized();
- Eigen::Vector2d neg_ik_hat_h_val =
- Eigen::Vector2d(k_normalized[1] * h_t_prime[1], k_normalized[0] * h_t_prime[1]);
- neg_ik_hat_h.emplace_back(neg_ik_hat_h_val);
+// Eigen::Vector2d neg_ik_hat_h_val =
+// Eigen::Vector2d(k_normalized[1] * h_t_prime[1], k_normalized[0] * h_t_prime[1]);
+// neg_ik_hat_h.emplace_back(neg_ik_hat_h_val);
+ double len = k_vector.norm();
+ if (len < .000001)
+ {
+ neg_ik_hat_h_x.emplace_back(0.0, 0.0);
+ neg_ik_hat_h_z.emplace_back(0.0, 0.0);
+ }
+ else
+ {
+ Eigen::Vector2d k_normalized = k_vector.normalized();
+ neg_ik_hat_h_x.emplace_back(k_normalized[0] * h_t_prime[1], k_normalized[0] * h_t_prime[1]);
+ neg_ik_hat_h_z.emplace_back(k_normalized[1] * h_t_prime[1], k_normalized[1] * h_t_prime[1]);
+ }
}
bool fast = true;
if (fast)
{
std::vector<Eigen::Vector2d> tmp = fast_fft(h_tildas);
- std::vector<Eigen::Vector2d> tmp2 = fast_fft(ikh);
- std::vector<Eigen::Vector2d> tmp3 = fast_fft(neg_ik_hat_h);
+// std::vector<Eigen::Vector2d> tmp2 = fast_fft(ikh);
+ // std::vector<Eigen::Vector2d> tmp3 = fast_fft(neg_ik_hat_h);
+ std::vector<Eigen::Vector2d> tmp4 = fast_fft(ikhx);
+ std::vector<Eigen::Vector2d> tmp5 = fast_fft(ikhz);
+ std::vector<Eigen::Vector2d> tmp6 = fast_fft(neg_ik_hat_h_x);
+ std::vector<Eigen::Vector2d> tmp7 = fast_fft(neg_ik_hat_h_z);
for (int i = 0; i < N; i++)
{
m_current_h[i] = tmp[i];
- m_slopes[i] = tmp2[i];
- m_displacements[i] = tmp3[i];
+ // m_slopes[i] = tmp2[i];
+ // m_displacements[i] = tmp3[i];
+ m_slopes_x[i] = tmp4[i];
+ m_slopes_z[i] = tmp5[i];
+ m_displacements_x[i] = tmp6[i];
+ m_displacements_z[i] = tmp7[i];
}
return;
@@ -95,8 +154,8 @@ void ocean_alt::fft_prime(double t){
for (int i=0; i<N; i++){
Eigen::Vector2d x_vector = m_waveIndexConstants[i].base_horiz_pos;
m_current_h[i] = Eigen::Vector2d(0.0, 0.0);
- m_displacements[i] = Eigen::Vector2d(0.0, 0.0);
- m_slopes[i] = Eigen::Vector2d(0.0, 0.0);
+ // m_displacements[i] = Eigen::Vector2d(0.0, 0.0);
+ // m_slopes[i] = Eigen::Vector2d(0.0, 0.0);
@@ -116,8 +175,8 @@ void ocean_alt::fft_prime(double t){
Eigen::Vector2d k_normalized = k_vector.normalized();
- m_displacements[i] += k_normalized*imag_comp;
- m_slopes[i] += -k_vector*imag_comp;
+ // m_displacements[i] += k_normalized*imag_comp;
+ // m_slopes[i] += -k_vector*imag_comp;
}
}
@@ -196,7 +255,7 @@ Eigen::Vector2d ocean_alt::get_k_vector(int n_prime, int m_prime){
double M_ = (double)num_cols;
double k_x = (2*M_PI*n_ - M_PI*N_)/Lx;
- double k_z = (2*M_PI*m_ - M_PI*M_)/Lz;
+ double k_z = (2*M_PI*m_ - M_PI*M_)/Lz;
return Eigen::Vector2d(k_x, k_z);
}
@@ -254,78 +313,166 @@ Eigen::Vector2d ocean_alt::complex_exp(double exponent){
return Eigen::Vector2d(real, imag);
}
-std::vector<Eigen::Vector3f> ocean_alt::get_vertices()
+void ocean_alt::update_ocean()
{
std::vector<Eigen::Vector3f> vertices = std::vector<Eigen::Vector3f>();
+ if (iterations < 10){
+ for (int i = 0; i < N; i++){
+ Eigen::Vector2d amplitude = m_current_h[i];
+ float height = amplitude[0];
+
+ if (height < min) min = height;
+ if (height > max) max = height;
+
+ }
+ iterations ++;
+ }
+ // reset normals & vertices arrays for the single tile
+ m_vertices = std::vector<Eigen::Vector3f>(N);
+ m_normals = std::vector<Eigen::Vector3f>(N);
for (int i = 0; i < N; i++){
Eigen::Vector2d horiz_pos = spacing*m_waveIndexConstants[i].base_horiz_pos;
Eigen::Vector2d amplitude = m_current_h[i];
float height = amplitude[0];
- Eigen::Vector2d slope = m_slopes[i] * .3f;
- Eigen::Vector3f s = Eigen::Vector3f(-slope[0], 0.0, -slope[1]);
- Eigen::Vector3f y = Eigen::Vector3f(0.0, 1.0, 0.0);
+ // Eigen::Vector2d slope = m_slopes[i] * .3f;
+ // Eigen::Vector3f s = Eigen::Vector3f(-slope[0], 0.0, -slope[1]);
+ // Eigen::Vector3f y = Eigen::Vector3f(0.0, 1.0, 0.0);
- float xs = 1.f + s[0]*s[0];
- float ys = 1.f + s[1]*s[1];
- float zs = 1.f + s[2]*s[2];
-
- Eigen::Vector3f diff = y - s;
- Eigen::Vector3f norm = Eigen::Vector3f(diff[0]/ sqrt(xs), diff[1]/ sqrt(ys), diff[2]/sqrt(zs));
+// float xs = 1.f + s[0]*s[0];
+// float ys = 1.f + s[1]*s[1];
+// float zs = 1.f + s[2]*s[2];
+//
+// Eigen::Vector3f diff = y - s;
+// Eigen::Vector3f Eigen::Vector3f(diff[0]/ sqrt(xs), diff[1]/ sqrt(ys), diff[2]/sqrt(zs));
// NEW
-// Eigen::Vector3f norm = Eigen::Vector3f(-slope[0], 1.0, -slope[1]);
- norm.normalize();
- //NEW
-
-
-
-
+ Eigen::Vector3f norm = Eigen::Vector3f(-m_slopes_x[i][0], 1.0, -m_slopes_z[i][0]);
+ norm = norm.normalized(); // FIXME: why do I have to be inverted?
//if (i==6) std::cout << amplitude[0] << std::endl;
// calculate displacement
- Eigen::Vector2d disp = lambda*m_displacements[i];
-
- //
+ // Eigen::Vector2d disp = lambda*m_displacements[i];
+ Eigen::Vector2d disp = lambda*Eigen::Vector2d(m_displacements_x[i][0], m_displacements_z[i][0])
+ + Eigen::Vector2d(vertex_displacement, vertex_displacement); // set corner at 0,0 for retiling
// for final vertex position, use the real number component of amplitude vector
- vertices.push_back(Eigen::Vector3f(horiz_pos[0] + disp[0], height, horiz_pos[1] + disp[1]));
+ m_vertices[i] = {horiz_pos[0] + disp[0], height, horiz_pos[1] + disp[1]};
m_normals[i] = norm.normalized();//Eigen::Vector3f(-slope[0], 1.0, -slope[1]).normalized();
//std::cout << "normal: " << m_normals[i] << std::endl
+ Eigen::Vector2i m_n = index_1d_to_2d(i);
+
+
+ // m_foam_constants.wavelengths[i] = 2.f* M_PI * m_slopes[i].dot(m_slopes[i]) / Lx;
+ float h_0 = m_waveIndexConstants[i].h0_prime[0]; // min*.2f;
+ float h_max = max*.01f; // the smaller the constant, the more foam there is
+ m_foam_constants.wavelengths[i] = (height - h_0 ) / (h_max - h_0);
+
+// if (i < 5){
+// std::cout << h_0 << ", " << h_max << std::endl;
+// std::cout << m_foam_constants.wavelengths[i] << std::endl;
+// }
+
+
+
}
- return vertices;
+
+ // populate foam constants
+ m_foam_constants.positions = vertices;
+}
+
+std::vector<Eigen::Vector3f> ocean_alt::get_vertices(){
+ // extend the returned array based on the tilecount
+ std::vector<Eigen::Vector3f> vertices = std::vector<Eigen::Vector3f>();
+ for (int i = 0; i < num_tiles_x; i++)
+ {
+ for (int j = 0; j < num_tiles_z; j++)
+ {
+ for (int k = 0; k < N; k++)
+ {
+ double c = Lx - 2 / (num_rows / Lx);
+ Eigen::Vector3f vertex = m_vertices[k] + Eigen::Vector3f(i*(c), 0.0, (j)*(c));
+ vertices.push_back(vertex);
+ }
+ }
+ }
+
+ return vertices;
}
std::vector<Eigen::Vector3f> ocean_alt::getNormals(){
- return m_normals;
+ // based on the tile count, add more to the normals
+ std::vector<Eigen::Vector3f> normals = std::vector<Eigen::Vector3f>();
+ // do the x 1D direction first
+ for (int i = 0; i < num_tiles_x; i++)
+ {
+ for (int j = 0; j < num_tiles_z; j++)
+ {
+ for (int k = 0; k < N; k++)
+ {
+ normals.push_back(m_normals[k]);
+ }
+ }
+ }
+
+
+ return normals;
}
std::vector<Eigen::Vector3i> ocean_alt::get_faces()
{
// connect the vertices into faces
std::vector<Eigen::Vector3i> faces = std::vector<Eigen::Vector3i>();
- for (int i = 0; i < N; i++)
- {
- int x = i / num_rows;
- int z = i % num_rows;
+ for (int i = 0; i < num_tiles_x; i++)
+ {
+ for (int j = 0; j < num_tiles_z; j++)
+ {
+ for (int k = 0; k < N; k++)
+ {
+ int x = k % num_rows;
+ int z = k / num_rows;
- // connect the vertices into faces
- if (x < num_rows - 1 && z < num_cols - 1)
- {
- int i1 = i;
- int i2 = i + 1;
- int i3 = i + num_rows;
- int i4 = i + num_rows + 1;
+ // connect the vertices into faces
+ if (x < num_rows - 1 && z < num_cols - 1)
+ {
+ int tile_index_offset = (j + num_tiles_z * i) * N;
+ int i1 = k + tile_index_offset;
+ int i2 = k + 1 + tile_index_offset;
+ int i3 = k + num_rows + tile_index_offset;
+ int i4 = k + num_rows + 1 + tile_index_offset;
+
+ faces.emplace_back(i2, i1, i3);
+ faces.emplace_back(i2, i3, i4);
+ }
+ }
+ }
+ }
+ return faces;
+
+
+// for (int i = 0; i < N; i++)
+// {
+// int x = i / num_rows;
+// int z = i % num_rows;
+//
+// // connect the vertices into faces
+// if (x < num_rows - 1 && z < num_cols - 1)
+// {
+// int i1 = i;
+// int i2 = i + 1;
+// int i3 = i + num_rows;
+// int i4 = i + num_rows + 1;
+//
// faces.emplace_back(i2, i1, i3);
// faces.emplace_back(i2, i3, i4);
- faces.emplace_back(i1, i2, i3);
- faces.emplace_back(i3, i2, i4);
- }
- }
- return faces;
+// faces.emplace_back(i1, i2, i3);
+// faces.emplace_back(i3, i2, i4);
+// }
+// }
+// return faces;
}
Eigen::Vector2d muliply_complex(Eigen::Vector2d a, Eigen::Vector2d b)
@@ -442,11 +589,11 @@ std::vector<Eigen::Vector2d> ocean_alt::fast_fft
}
// divide by N*N and add the signs based on the indices
- double sign[] = {1, -1};
+ double sign[] = {1.0, -1.0};
for (int i = 0; i < N; i++)
{
- // h[i] /= N;
- h[i] /= sqrt(N);
+ h[i] /= N;
+ // h[i] /= sqrt(N);
h[i] *= sign[(i / num_rows + i % num_cols) % 2];
}
diff --git a/src/ocean/ocean_alt.h b/src/ocean/ocean_alt.h
index 4bb1d54..dc8888d 100644
--- a/src/ocean/ocean_alt.h
+++ b/src/ocean/ocean_alt.h
@@ -20,16 +20,28 @@ struct WaveIndexConstant{
Eigen::Vector2d k_vector = Eigen::Vector2d(0.f, 0.f); // static horiz pos with no displacement
};
+struct FoamConstants{
+ std::vector<Eigen::Vector3f> positions;
+ std::vector<Eigen::Vector2f> k_vectors;
+ std::vector<float> wavelengths;
+ std::vector<Eigen::Vector2f> texCoords;
+};
+
class ocean_alt
{
public:
ocean_alt();
void updateVertexAmplitudes(double t);
+ void update_ocean();
std::vector<Eigen::Vector3f> get_vertices();
std::vector<Eigen::Vector3i> get_faces();
void fft_prime(double t);
std::vector<Eigen::Vector3f> getNormals();
+ FoamConstants getFoamConstants(){
+ return m_foam_constants;
+ }
+
@@ -48,44 +60,53 @@ private:
Eigen::Vector2d get_horiz_pos(int i);
std::pair<double, double> sample_complex_gaussian();
-
-
-
-
-
-
-
+ // FOAM
+ std::vector<float> m_saturation;
std::map<int, WaveIndexConstant> m_waveIndexConstants; // stores constants that only need to be calculate once for each grid constant
-
-
const double Lx = 512.0;
const double Lz = 512.0;
const int num_rows = 256;
const int num_cols = 256;
+ const int num_tiles_x = 2;
+ const int num_tiles_z = 2;
+
+ const double vertex_displacement = Lx / 2;
+
const int N = num_rows*num_cols; // total number of grid points
- const double lambda = 0; // how much displacement matters
- const double spacing = 35.0; // spacing between grid points
+ const double lambda = .5; // how much displacement matters
+ const double spacing = 1.0; // spacing between grid points
- const double A = 100; // numeric constant for the Phillips spectrum
- const double V = 50; // wind speed
+ const double A = 10000; // numeric constant for the Phillips spectrum
+ const double V = 500; // wind speed
const double gravity = 9.81;
const double L = V*V/gravity;
const Eigen::Vector2d omega_wind = Eigen::Vector2d(1.0, 0.0); // wind direction, used in Phillips equation
std::vector<Eigen::Vector2d> m_current_h; // current height fields for each K
- std::vector<Eigen::Vector2d> m_displacements; // current displacement vector for each K
- std::vector<Eigen::Vector2d> m_slopes; // current displacement vector for each K
+ // std::vector<Eigen::Vector2d> m_displacements; // current displacement vector for each K
+ // std::vector<Eigen::Vector2d> m_slopes; // current displacement vector for each K
+ std::vector<Eigen::Vector2d> m_slopes_x;
+ std::vector<Eigen::Vector2d> m_slopes_z;
+ std::vector<Eigen::Vector2d> m_displacements_x;
+ std::vector<Eigen::Vector2d> m_displacements_z;
//std::vector<Eigen::Vector3f> m_slope_vectors; // current displacement vector for each K
- std::vector<Eigen::Vector3f> m_normals; // current displacement vector for each K
+ std::vector<Eigen::Vector3f> m_normals; // normal calculations
+
+ // FOR FOAM:
+ FoamConstants m_foam_constants;
+ float max = 0;
+ float min = 0;
+ int iterations = 0;
+ std::vector<Eigen::Vector3f> m_vertices; // current displacement vector for each K
const double D = 1.0; // Depth below mean water level (for dispersion relation)
generated by cgit v1.2.3-70-g09d2 (git 2.46.0) at 2025-10-06 07:11:39 +0000