SH Shapes

Spherical Harmonic Shapes

Log in to post a comment.

// You can find the Turtle API reference here: https://turtletoy.net/syntax
Canvas.setpenopacity(0.05);

// Global code will be evaluated once.
const turtle = new Turtle();

const PI = Math.PI;
const PI2 = Math.PI*2;
const sqrt2 = Math.sqrt(2);

const RANDOM = 1;

const WIREFRAME = 1

// Change COEF to combine the linear combination of
let COEF = [
            0.1,
          0.1,0.1,0.1,
        1,0.1,0.1,0.1,1,
      0.01,0.1,0.1,0.01,0.01,0.01,0.1,
    0.1,0.001,0.001,0.001,0.001,0.001,0.001,0.001,1,
  0.01,0.01,0.01,0.01,0.001,0.01,0.01,0.01,0.001,0.001,0.001,
1,0.001,0.001,0.001,0.0001,0.001,0.0001,0.001,0.0001,0.001,0.0001,0.001,1,
//0.0001,0.0001,0.0001,0.0001,1,0.0001,0.0001,0.0001,0.0001,0.0001,0.0001,0.0001,0.0001,0.0001,0.0001
];

const LEVEL = Math.sqrt(COEF.length); 

let scale = [1.0,1.0,1.0]
let rotation = quat(0.0,3.1,3.1);
let position = [0,0,1.5];

// Begin SH Math
function factorial(x)
{
    if(x==0){return 1;}
    else{return x*factorial(x-1);}
}

function P(l,m,x)
{
    let pmm = 1.0;
    if(m>0)
    {
        let somx2 = Math.sqrt((1-x)*(1+x));
        let fact = 1.0;
        for(let i=1;i<=m;i++)
        {
            pmm *= (-fact)*somx2;
            fact += 2.0;
        }
        
    }
    if(l==m) return pmm;
    let pmmp1 = x*(2*m+1.0)*pmm;
    if(l==m+1){return pmmp1;}
    let p11 = 0.0;
    
    for(let i = m+2; i<=l;i++)
    {
        pll = ((2.0*11-1.0)*x*pmmp1 - (i+m-1)*pmm) / (i-m);
        pmm = pmmp1;
        pmmp1 = pll;
    }
    return pll;
}

function K(l,m)
{
    let temp = ((2.0*l+1.0)*factorial(l-m))/(4.0*PI*factorial(l+m));
    return Math.sqrt(temp);
}

function SH(l,m,theta,phi)
{
    // theta in range [0...PI]
    // phi in range [0...2PI]
    if(m==0){return K(1,0)*P(l,m,Math.cos(theta));}
    else if(m>0){return sqrt2*K(l,m)*Math.cos(m*phi)*P(l,m,Math.cos(theta));}
    else
    {
        return sqrt2*K(l,-m)*Math.sin(-m*phi)*P(l,-m,Math.cos(theta));
    }
}

function SHBasisIndex(l,m)
{
    return l*(l+1) + m;
}




if(RANDOM)
{
    for(let i=0;i<COEF.length;i++)
    {
        COEF[i] = COEF[i]*Math.random()*0.5;
    }
}


function GetSHScale(theta,phi)
{
    let sum = 0;
    for(var l=0;l<LEVEL;l++)
    {
        for(var m=-l; m<=l;m++)
        {
            let index = SHBasisIndex(l,m);
            sum += COEF[index]*SH(l,m,theta,phi);
        }
    }
    
    return Math.abs(sum);
}

// End SH Math

// Quaternion from Eular angles
// https://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles
function quat(yaw,pitch,roll)
{
    let cy = Math.cos(yaw*0.5);
    let sy = Math.sin(yaw*0.5);
    let cp = Math.cos(pitch*0.5);
    let sp = Math.sin(pitch*0.5);
    let cr = Math.cos(roll*0.5);
    let sr = Math.sin(roll*0.5);
    
    let w = cy * cp * cr + sy * sp * sr;
    let x = cy * cp * sr - sy * sp * cr;
    let y = sy * cp * sr + cy * sp * cr;
    let z = sy * cp * cr - cy * sp * sr;
    
    return [x,y,z,w];
}

// ref https://github.com/toji/gl-matrix/blob/master/src/mat4.js
// fromRotationTranslationScale
// matrix is stored in the row pattern
function make_transform(o,pos,quat,scale)
{
    // scale temp vars
    let sx = scale[0];
    let sy = scale[1];
    let sz = scale[2];
    
    // translation temp vars
    let tx = pos[0];
    let ty = pos[1];
    let tz = pos[2];
    
    //Quat calc temp variables
    let x = quat[0];
    let y = quat[1];
    let z = quat[2];
    let w = quat[3];
    
    let xx = x*x*2.0;
    let xy = x*y*2.0;
    let xz = x*z*2.0;
    let wx = x*w*2.0;
    
    let yy = y*y*2.0;
    let yz = y*z*2.0;
    let wy = y*w*2.0;
    
    let zz = z*z*2.0;
    let wz = z*w*2.0;

    o[0] = (1-yy-zz)*sx;
    o[1] = (xy-wz)*sy;
    o[2] = (xz+wy)*sz;
    o[3] = tx;
    o[4] = (xy+wz)*sx;
    o[5] = (1-xx-zz)*sy;
    o[6] = (yz-wx)*sz;
    o[7] = ty;
    o[8] = (xz-wy)*sx;
    o[9] = (yz+wx)*sy;
    o[10] = (1-xx-yy)*sz;
    o[11] = tz;
    o[12] = 0;
    o[13] = 0;
    o[14] = 0;
    o[15] = 1;
    
    return o;
}

function transform_vector(v,m)
{
    let x = v[0], y = v[1], z = v[2], w = v[3];
    let o = [0, 0, 0, 0];
    o[0] = m[0] * x + m[1] * y + m[2] * z + m[3] * w;
    o[1] = m[4] * x + m[5] * y + m[6] * z + m[7] * w;
    o[2] = m[8] * x + m[9] * y + m[10] * z + m[11] * w;
    o[3] = m[12] * x + m[13] * y + m[14] * z + m[15] * w;
    return o;
}

let focus = 200;
function project(v){
    let x = v[0];
    let y = v[1];
    let z = v[2];
    let s = focus/z;
    return [x*s,y*s];
}

//transfrom matrix
let transform = new Float32Array(16);

let Mesh = function(verts,edges,tsf)
{
    let out = {}
    out.verts = verts;
    out.edges = edges;
    out.transform = tsf
    return out;
}



// fill in the SH shape vertice and edge list
function make_sh_shape(nu,nv)
{
    let o = {};
    o.edges = [];
    o.verts = [];
    let du = 1.0/nu;
    let dv = 1.0/nv;
    nu +=1;
    nv +=1;
    let u = 0,v = 0;
    let i = 0,j = 0;
    
    if(!WIREFRAME)
    {
        o.verts.push([0,0,0,1.0]);
    }
    
    // push vertice
    for(let i=0;i<nu;i++)
    {
        u = i*du;
        for(let j=0;j<nv;j++)
        {
            v = j*dv;
            
            let cu = Math.cos(u*PI2);
            let su = Math.sin(u*PI2);
            let cv = Math.cos(v*PI);
            let sv = Math.sin(v*PI);
            
            let x = sv*cu;
            let y = sv*su;
            let z = cv;
            
            let sh_scale = GetSHScale(v*PI,u*PI2);
            
            
            o.verts.push([x*sh_scale,y*sh_scale,z*sh_scale,1.0]);
            
            //o.verts.push([u,v,1.0,1.0]);
        }
    }
    
    //link edges
    for(let i=0;i<nu;i++)
    {
        u = i*du;
        for(let j=0;j<nv;j++)
        {
            v = j*dv;
            let id = i*nv+j;
            let inext = (i+1)%nu;
            let jnext = (j+1)%nv;
            let id_down = i*nv+jnext;
            let id_right = inext*nv+j;
            if(WIREFRAME)
            {
                o.edges.push(id,id_down);
                o.edges.push(id,id_right);
            }
            else
            {
                o.edges.push(0,id);
            }
            
        }
        
    }
    return o;
}

function draw_mesh(tctx,mesh)
{
    let edges = mesh.edges;
    let verts = mesh.verts;
    let tsf = mesh.transform;
    
    for (let i = 0, len = edges.length/2; i < len; i++) {
        
        let i0 = edges[2*i];
        let i1 = edges[2*i+1];
        
        // perform transform
        let v0 = transform_vector(verts[i0],tsf);
        let v1 = transform_vector(verts[i1],tsf);
        
        // projection
        v0 = project(v0);
        v1 = project(v1);
        
        // draw edge
        tctx.penup();
		tctx.goto([v0[0],v0[1]]);
		tctx.pendown();
		tctx.goto([v1[0],v1[1]]);
		
    }
}


function draw_mesh_substep(tctx,mesh,step)
{
    let edges = mesh.edges;
    let len = edges.length/2;
    if(step>=len-2)
    {
        return false;
    }
    
    
    let verts = mesh.verts;
    let tsf = mesh.transform;
    
    let i0 = edges[2*step];
    let i1 = edges[2*step+1];
        
        // perform transform
    let v0 = transform_vector(verts[i0],tsf);
    let v1 = transform_vector(verts[i1],tsf);
        
    // projection
    v0 = project(v0);
    v1 = project(v1);
        
    // draw edge
    tctx.penup();
	tctx.goto([v0[0],v0[1]]);
	tctx.pendown();
	tctx.goto([v1[0],v1[1]]);
	return true;
}


make_transform(transform,position,rotation,scale);

let btl = make_sh_shape(200,100);
let bottle = Mesh(btl.verts,btl.edges,transform);

// The walk function will be called until it returns false.
function walk(i) {
   
    //draw_mesh(turtle,bottle);
    let stop = draw_mesh_substep(turtle,bottle,i);
    return stop;
}