import * as THREE from 'three';
import { OrbitControls } from 'OrbitControls';
var gui = new dat.GUI();
const canvas = document.querySelector('#c');
var scene = new THREE.Scene();
var raycaster = new THREE.Raycaster();
//create some camera
var camera = new THREE.PerspectiveCamera(55, window.innerWidth / window.innerHeight, 0.1, 1000);
camera.position.x = 0;
camera.position.y = 0;
camera.position.z = 1;
camera.lookAt(0, 0, 0);
var renderer = new THREE.WebGLRenderer({
canvas,
antialias: true
});
renderer.setSize(window.innerWidth, window.innerHeight);
renderer.setPixelRatio(2);
renderer.setClearColor(new THREE.Color(0x595959));
document.body.appendChild(renderer.domElement);
var controls = new OrbitControls(camera, renderer.domElement);
// white spotlight shining from the side, casting a shadow
var spotLight = new THREE.SpotLight(0xffffff, 2.5, 25, Math.PI / 6);
spotLight.position.set(4, 10, 7);
scene.add(spotLight);
//
function lerp(v0, v1, u) {
return v0 + (v1-v0) * u
}
function mod(x, y){
return ((x % y) + y) % y;
}
function clamp(a, min, max){
return Math.max(Math.min(a, max), min);
}
function fitrange(value, min1, max1, min2, max2){
return min2 + (value - min1) * (max2 - min2) / (max1 - min1);
}
function arg_min_fancy(a){ //slow
return a.reduce((iMin, x, i, arr) => x < arr[iMin] ? i : iMin, 0);
}
function arg_min(a){
let min = 1e10;
let index = -1;
for(let i = 0; i < a.length; i++){
if(a[i] < min){
min = a[i];
index = i;
}
}
return index;
}
function index_into_list(Q, A){
let out = [];
for(let i = 0; i < Q.length; i++){
out.push(A[Q[i]]);
}
return out;
}
//does argmin and index into list all at once for speed
function get_nearest_Q(Q, dists){
let min = 1e10;
let index = -1;
for(let i = 0; i < Q.length; i++){
if(dists[Q[i]] < min){
min = dists[Q[i]];
index = i;
}
}
return index;
}
function perlin2_grad(pos){
let h = .01;
let dx = noise.perlin2(pos.x - h, pos.z) - noise.perlin2(pos.x + h, pos.z);
let dy = noise.perlin2(pos.x, pos.z - h) - noise.perlin2(pos.x, pos.z + h);
return new THREE.Vector3(dx / h, 0, dy / h);
}
function perlin2_grad2(pos){
let h = .01;
let dx = noise.perlin2(pos.x - h, pos.y) - noise.perlin2(pos.x + h, pos.y);
let dy = noise.perlin2(pos.x, pos.y - h) - noise.perlin2(pos.x, pos.y + h);
return new THREE.Vector2(dx / h, dy / h);
}
function outerproduct2(u){
return [u.x * u.x, u.x * u.y,
u.y * u.x, u.y * u.y];
}
function add2(u, v){
return [u[0] + v[0], u[1] + v[1], u[2] + v[2], u[3] + v[3]]
}
function determinant2(q){
return q[0] * q[3] - q[1] * q[2];
}
function eigvals2(q){
//https://www.johndcook.com/blog/2021/05/07/trick-for-2x2-eigenvalues/
let m = (q[0] + q[3]) / 2.0;
let d = determinant2(q);
return [m + Math.sqrt(m * m - d), m - Math.sqrt(m * m - d)];
}
function eigvectors2(q, eigvals){
//https://people.math.harvard.edu/~knill/teaching/math21b2004/exhibits/2dmatrices/index.html
let [g0, g1] = eigvals;
if(q[3] != 0){
return [[g0 - q[3], q[2]], [g1 - q[3], q[2]]];
}
if(q[2] != 0){
return [ [q[1], g0 - q[0]], [q[1], g1 - q[0]] ];
}
return [[1,0],[0,1]];
}
function normalize2(u){
let mag = Math.sqrt(u[0] * u[0] + u[1] * u[1]);
if(mag == 0.0){
return [0,0];
}
return [u[0] / mag, u[1] / mag]
}
function make_metric(vec2_in, anisotropy){
let q = outerproduct2(vec2_in);
let [g0,g1] = eigvals2(q);
let [e0, e1] = eigvectors2(q, [g0,g1]);
e0 = normalize2(e0); e1 = normalize2(e1);
function gamma(x, n){
return Math.pow(n + Math.abs(x), -1);
}
let n = 1 - anisotropy;
if(n == 0) n += 1e-7;
let gamma_g0 = gamma(g0, n);
let gamma_g1 = gamma(g1, n);
let h0 = outerproduct2(new THREE.Vector2(e0[0] * gamma_g0, e0[1] * gamma_g0));
let h1 = outerproduct2(new THREE.Vector2(e1[0] * gamma_g1, e1[1] * gamma_g1));
return add2(h0, h1);
}
function vec2m2x2(A, u){
return new THREE.Vector2(A[0] * u.x + A[2] * u.y, A[1] * u.x + A[3] * u.y);
}
function generate_pointset(point_count, size){
var points2d = [];
for (let i = 0; i < point_count; i++) {
let x = THREE.MathUtils.randFloatSpread(size.x);
let y = THREE.MathUtils.randFloatSpread(size.y);
points2d.push(new THREE.Vector2(x, y));
}
return points2d;
}
function curl_noise(pos, scale = 4){
let noise = perlin2_grad2(pos.clone().multiplyScalar(scale));
return new THREE.Vector2(noise.y, -noise.x);
//return new THREE.Vector2(1,0);
}
/*
left to implement:
validate the metric stuff works by comparing output values with stuff from houdini!!
*/
class Partio{
constructor(Geometry){
this.geometry = Geometry;
this.vert_count = this.geometry.getAttribute("position").count;
this.directions = Array(this.vert_count);
this.metrics = Array(this.vert_count);
this.npts = [];
this.npts_dists = [];
}
//might be faster to do this inline with the other calls instead of as one big ol thingy
nearpoints_all(radius, max_nb_count){
let n = this.vert_count;
let dists = Array.from({length: n}, e => Array(max_nb_count).fill(1e14)); //make nested filled arrays
let nearpoints = Array.from({length: n}, e => Array(max_nb_count).fill(-1));
for( let i = 0; i < n; i++){
for(let j = 0; j < n; j++){
if(i == j) continue;
let j_dist = this.P(i).distanceToSquared(this.P(j));
if(j_dist > radius) continue;
if(j_dist > dists[i][max_nb_count -1] && dists[i][max_nb_count -1] != -1) continue; //early skip
for(let k = 0; k < max_nb_count; k++){
if(j_dist < dists[i][k] || dists[i][k] == -1){
dists[i].splice(k, 0, j_dist);
nearpoints[i].splice(k,0,j);
dists[i].pop();
nearpoints[j].pop();
break;
}
}
}
}
this.npts = nearpoints;
this.npts_dists = dists;
}
nearpoints(position, radius, max_nb_count){
let n = this.vert_count;
let dists = Array(max_nb_count).fill(1e14); //make nested filled arrays
let nearpoints = Array(max_nb_count).fill(-1);
for(let j = 0; j < n; j++){
let j_dist = position.distanceToSquared(this.P(j));
if(j_dist > radius) continue;
if(j_dist > dists[max_nb_count -1] && dists[max_nb_count -1] != -1) continue; //early skip
for(let k = 0; k < max_nb_count; k++){
if(j_dist < dists[k] || dists[k] == -1){
dists.splice(k, 0, j_dist);
nearpoints.splice(k,0,j);
dists.pop();
nearpoints.pop();
break;
}
}
}
return nearpoints;
}
kd_dist(a, b){
//return a.distanceToSquared(b);
return Math.pow(a.x - b.x, 2) + Math.pow(a.y - b.y, 2);
}
init_kd(){
let positions = [];
for(let i = 0; i < this.vert_count; i++){
let pos = this.P(i);
positions.push({index: i, x:pos.x, y:pos.y});
}
return new kdTree(positions, this.kd_dist, ["x", "y"])
}
build_direction_attrib(funct){
for(let i = 0; i < this.vert_count; i++){
this.directions[i] = funct(this.P(i)); //needs to be an array of vector twos
}
}
generate_metrics(anisotropy){
for(let i = 0; i < this.vert_count; i++){
this.metrics[i] = make_metric(this.directions[i], anisotropy);
}
}
update_step(radius, max_nb_count, anisotropy, time_inc){
let kd = this.init_kd();
this.generate_metrics(anisotropy);
let outP = Array(this.vert_count);
let outDir = Array(this.vert_count);
for(let i = 0; i < this.vert_count; i++){
let out_vel = new THREE.Vector2(0,0);
let out_aniso_dir = this.directions[i].clone();//new THREE.Vector2(0,0);
let weights = 1.0;
let p0 = this.P(i);
//let npts = this.nearpoints(p0, radius, max_nb_count);
let npts = kd.nearest({x:p0.x, y:p0.y}, max_nb_count, [radius]);
for(let j = 0; j < npts.length; j++){
let npt = npts[j][0].index; //near point
//let npt = npts[j];
if(npt == i) continue;
/*
if(npt == -1){
continue;
}*/
let npt_metric = this.metrics[npt];
let npt_aniso_dir = this.directions[npt];
let p1 = this.P(npt);
//console.log(npt, npt_metric);
let dir = p1.clone().sub(p0);
let weight = dir.dot(vec2m2x2(npt_metric, dir));
out_aniso_dir.add(npt_aniso_dir.clone().multiplyScalar(weight));
let normalized_aniso_dir = npt_aniso_dir.clone().normalize();
let projection = normalized_aniso_dir.dot(dir.clone().negate());
let scaled_projection = normalized_aniso_dir.clone().multiplyScalar(projection);
out_vel.add((dir.clone().add(scaled_projection).multiplyScalar(weight))); //inplace updates suck!
weights += weight;
}
if(weights != 0){
out_vel.divideScalar(weights);
out_aniso_dir.divideScalar(weights);
outDir[i] = this.directions[i].clone().lerp(out_aniso_dir, time_inc);
}else{
//console.log("weight is zero, NPTS: " + npts);
outDir[i] = this.directions[i];
}
outP[i] = p0.clone().add(out_vel.multiplyScalar(time_inc));
}
this.geometry.setFromPoints(outP);
this.geometry.attributes.position.needsUpdate = true;
for(let i = 0; i < this.vert_count; i++){
this.directions[i] = outDir[i];
}
}
neighbours(index){
return this.neighbours_list[index];
}
P(index){
return new THREE.Vector2(this.geometry.getAttribute("position").array[index * 3 + 0],
this.geometry.getAttribute("position").array[index * 3 + 1]);
//this.geometry.getAttribute("position").array[index * 3 + 2]);
}
setFromPoints(points2d){
this.geometry.setFromPoints(points2d);
this.geometry.attributes.position.needsUpdate = true;
this.vert_count = this.geometry.getAttribute("position").count;
this.directions = Array(this.vert_count);
this.metrics = Array(this.vert_count);
}
}
var size = { x: 1, y: 1 };
var pointsCount = 2000;
let points2d = generate_pointset(pointsCount, size);
var geometry = new THREE.BufferGeometry().setFromPoints(points2d);
geometry.attributes.position.setUsage( THREE.DynamicDrawUsage );
let partio = new Partio(geometry);
partio.build_direction_attrib(x => {return curl_noise(x, 2)})
var cloud = new THREE.Points(
geometry,
new THREE.PointsMaterial({ color: 0x99ccff, size: .01 })
);
scene.add( cloud );
function make_viz_geo_from_partio(partio, viz_scale){
var viz_geo = new THREE.BufferGeometry()
let indices = [];
let verts = [];
let index = 0;
for(let i = 0; i < partio.vert_count; i++){
indices.push(index++, index++);
let p0 = partio.P(i);
let dir = partio.directions[i].clone();
let p1 = p0.clone().add(dir.multiplyScalar(viz_scale));
verts.push(p0, p1);
}
viz_geo.setFromPoints(verts);
viz_geo.setIndex(indices);
let viz_material = new THREE.LineBasicMaterial( {
transparent: true,
color: "yellow",
linewidth: 1,
linecap: 'round', //ignored by WebGLRenderer
linejoin: 'round' //ignored by WebGLRenderer
} );
return new THREE.LineSegments(viz_geo, viz_material);
}
let viz_scale = .01;
let vector_visualizer = make_viz_geo_from_partio(partio, viz_scale)
vector_visualizer.name = "vec_viz";
scene.add(vector_visualizer)
function update_viz_from_partio(viz_geometry, partio, viz_scale){
let verts = [];
for(let i = 0; i < partio.vert_count; i++){
let p0 = partio.P(i);
let dir = partio.directions[i].clone();
let p1 = p0.clone().add(dir.multiplyScalar(viz_scale));
verts.push(p0, p1);
}
viz_geometry.setFromPoints(verts);
viz_geometry.attributes.position.needsUpdate = true;
}
function reset_full_simulation(scene, partio, num_partios, viz_geometry, viz_scale, dir_funct){
var selectedObject = scene.getObjectByName(viz_geometry.name);
scene.remove( selectedObject );
let points2d = generate_pointset(num_partios, size);
partio.setFromPoints(points2d);
partio.build_direction_attrib(dir_funct);
partio.geometry.attributes.position.needsUpdate = true;
viz_geometry = make_viz_geo_from_partio(partio, viz_scale);
viz_geometry.name = "vec_viz";
scene.add(viz_geometry);
}
let direction_function_selector = {
0: x => {return curl_noise(x, gui_parms_static["Noise Scale"])},
1: x => {return new THREE.Vector2(1,0)},
2: x => {return new THREE.Vector2(0,1)},
}
let gui_parms_static = {
"Number of Particles": 2000,
"Toggle Vector Visualizers": true,
"Nearpoints Radius": .1,
"Max Nearpoints": 10,
"Visualizer Scale": .01,
"Noise Function": 0,
"Noise Scale": 2,
}
let gui_parms_dynamic = {
partio_updates: 1,
"Play/Pause Simulation": function() { this.partio_updates = (this.partio_updates + 1) % 2 },
"Reset Simulation": function() { reset_full_simulation(scene,
partio,
gui_parms_static["Number of Particles"],
vector_visualizer,
.01,
direction_function_selector[gui_parms_static["Noise Function"]]) },
}
gui.add(gui_parms_dynamic, "Play/Pause Simulation");
gui.add(gui_parms_dynamic, "Reset Simulation");
const noise_folder = gui.addFolder("Noise Settings")
noise_folder.add(gui_parms_static, "Noise Function", {"Curl Noise":0, "X+":1, "Y+":2}).onChange(x => { reset_full_simulation(scene, partio, gui_parms_static["Number of Particles"], vector_visualizer, .01, direction_function_selector[gui_parms_static["Noise Function"]])} );
noise_folder.add(gui_parms_static, "Noise Scale", 1, 5).onChange(x => { reset_full_simulation(scene, partio, gui_parms_static["Number of Particles"], vector_visualizer, .01, direction_function_selector[gui_parms_static["Noise Function"]])} );
gui.add(gui_parms_static, "Number of Particles", 100, 10000, 1).onChange(num_partios => { reset_full_simulation(scene, partio, num_partios, vector_visualizer, .01, direction_function_selector[gui_parms_static["Noise Function"]])} );
gui.add(gui_parms_static, "Max Nearpoints", 3, 100, 1);
gui.add(gui_parms_static, "Nearpoints Radius", .01, 1);
const viz_folder = gui.addFolder("Visualizer Settings")
viz_folder.add(gui_parms_static, "Visualizer Scale", .001, .1);
viz_folder.add(gui_parms_static, "Toggle Vector Visualizers");
function animate() {
requestAnimationFrame(animate);
var viz_geo = scene.getObjectByName(vector_visualizer.name);
if(gui_parms_dynamic.partio_updates){
partio.update_step(gui_parms_static["Nearpoints Radius"], gui_parms_static["Max Nearpoints"], 1, 1.0/24.0);
partio.geometry.attributes.position.needsUpdate = true;
update_viz_from_partio(viz_geo.geometry, partio, gui_parms_static["Visualizer Scale"]);
}
if(gui_parms_static["Toggle Vector Visualizers"]){
viz_geo.material.opacity = 1;
}else{
viz_geo.material.opacity = 0;
}
controls.update();
renderer.render(scene, camera);
};
animate();
//*/
!
