Spherical Harmonics is widely used in Computer Graphics. They are analogue to Fourier basis on a sphere, consists of a set of orthogonal functions to represent functions defined on the surface of a sphere.
However, they are very tricky to implement due to lots of constants and integral functions.
Here is a real-time visualization that forks iq’s raytracing of SH functions, but adding the fourth order of them.
Code:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 |
// Created by inigo quilez - iq/2013 // Added the 4th order of SH - ruofei/2017 // License Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License. // Five bands of Spherical Harmonics functions (or atomic orbitals if you want). // For reference and fun. // antialias level (try 1, 2, 3, ...) #define AA 1 // change TIME to 0 to stop rotating #define TIME iGlobalTime //#define SHOW_SPHERES //--------------------------------------------------------------------------------- // Constants, see here: http://en.wikipedia.org/wiki/Table_of_spherical_harmonics #define k01 0.2820947918 // sqrt( 1/PI)/2 #define k02 0.4886025119 // sqrt( 3/PI)/2 #define k03 1.0925484306 // sqrt( 15/PI)/2 #define k04 0.3153915652 // sqrt( 5/PI)/4 #define k05 0.5462742153 // sqrt( 15/PI)/4 #define k06 0.5900435860 // sqrt( 70/PI)/8 #define k07 2.8906114210 // sqrt(105/PI)/2 #define k08 0.4570214810 // sqrt( 42/PI)/8 #define k09 0.3731763300 // sqrt( 7/PI)/4 #define k10 1.4453057110 // sqrt(105/PI)/4 #define k11 2.5033429418 // Math.sqrt(35/Math.PI) * 3 / 4 #define k12 1.7701307698 // Math.sqrt(70/Math.PI) * 3 / 8 #define k13 0.9461746958 // Math.sqrt(5/Math.PI) * 3 / 4 #define k14 0.6690465436 // Math.sqrt(10/Math.PI) * 3 / 8 #define k15 0.1057855469 // Math.sqrt(1/Math.PI) * 3 / 16 #define k16 0.4730873479 // Math.sqrt(5/Math.PI) * 3 / 8 #define k17 0.6258357354 // Math.sqrt(35/Math.PI) * 3 / 16 // Y_l_m(s), where l is the band and m the range in [-l..l] float SH( in int l, in int m, in vec3 s ) { vec3 n = s.zxy; //---------------------------------------------------------- if( l==0 ) return k01; //---------------------------------------------------------- if( l==1 && m==-1 ) return -k02*n.y; if( l==1 && m== 0 ) return k02*n.z; if( l==1 && m== 1 ) return -k02*n.x; //---------------------------------------------------------- if( l==2 && m==-2 ) return k03*n.x*n.y; if( l==2 && m==-1 ) return -k03*n.y*n.z; if( l==2 && m== 0 ) return k04*(3.0*n.z*n.z-1.0); if( l==2 && m== 1 ) return -k03*n.x*n.z; if( l==2 && m== 2 ) return k05*(n.x*n.x-n.y*n.y); //---------------------------------------------------------- if( l==3 && m==-3 ) return -k06*n.y*(3.0*n.x*n.x-n.y*n.y); if( l==3 && m==-2 ) return k07*n.z*n.y*n.x; if( l==3 && m==-1 ) return -k08*n.y*(5.0*n.z*n.z-1.0); if( l==3 && m== 0 ) return k09*n.z*(5.0*n.z*n.z-3.0); if( l==3 && m== 1 ) return -k08*n.x*(5.0*n.z*n.z-1.0); if( l==3 && m== 2 ) return k10*n.z*(n.x*n.x-n.y*n.y); if( l==3 && m== 3 ) return -k06*n.x*(n.x*n.x-3.0*n.y*n.y); //---------------------------------------------------------- return 0.0; } // unrolled version of the above float SH_0_0( in vec3 s ) { vec3 n = s.zxy; return k01; } float SH_1_0( in vec3 s ) { vec3 n = s.zxy; return -k02*n.y; } float SH_1_1( in vec3 s ) { vec3 n = s.zxy; return k02*n.z; } float SH_1_2( in vec3 s ) { vec3 n = s.zxy; return -k02*n.x; } float SH_2_0( in vec3 s ) { vec3 n = s.zxy; return k03*n.x*n.y; } float SH_2_1( in vec3 s ) { vec3 n = s.zxy; return -k03*n.y*n.z; } float SH_2_2( in vec3 s ) { vec3 n = s.zxy; return k04*(3.0*n.z*n.z-1.0); } float SH_2_3( in vec3 s ) { vec3 n = s.zxy; return -k03*n.x*n.z; } float SH_2_4( in vec3 s ) { vec3 n = s.zxy; return k05*(n.x*n.x-n.y*n.y); } float SH_3_0( in vec3 s ) { vec3 n = s.zxy; return -k06*n.y*(3.0*n.x*n.x-n.y*n.y); } float SH_3_1( in vec3 s ) { vec3 n = s.zxy; return k07*n.z*n.y*n.x; } float SH_3_2( in vec3 s ) { vec3 n = s.zxy; return -k08*n.y*(5.0*n.z*n.z-1.0); } float SH_3_3( in vec3 s ) { vec3 n = s.zxy; return k09*n.z*(5.0*n.z*n.z-3.0); } float SH_3_4( in vec3 s ) { vec3 n = s.zxy; return -k08*n.x*(5.0*n.z*n.z-1.0); } float SH_3_5( in vec3 s ) { vec3 n = s.zxy; return k10*n.z*(n.x*n.x-n.y*n.y); } float SH_3_6( in vec3 s ) { vec3 n = s.zxy; return -k06*n.x*(n.x*n.x-3.0*n.y*n.y); } float SH_4_0( in vec3 s ) { vec3 n = s.zxy; return k11 * (n.x*n.y * (n.x*n.x - n.y*n.y)); } float SH_4_1( in vec3 s ) { vec3 n = s.zxy; return k12 * (3.0*n.x*n.x - n.y*n.y) * n.y * n.z; } float SH_4_2( in vec3 s ) { vec3 n = s.zxy; return k13 * (n.x*n.y * (7.0*n.z*n.z-dot(n,n)) ); } float SH_4_3( in vec3 s ) { vec3 n = s.zxy; return k14 * (n.z*n.y * (7.0*n.z*n.z-3.0*dot(n,n)) ); } float SH_4_4( in vec3 s ) { vec3 n = s.zxy; float z2 = n.z*n.z, r2 = dot(n, n); return k15 * (35.0 * z2*z2 - 30.0 * z2*r2 + 3.0 * r2*r2); } float SH_4_5( in vec3 s ) { vec3 n = s.zxy; return k14 * (n.z*n.x * (7.0*n.z*n.z-3.0*dot(n,n)) ); } float SH_4_6( in vec3 s ) { vec3 n = s.zxy; return k16 * ( (n.x*n.x-n.y*n.y)*(7.0*n.z*n.z-dot(n,n)) ); } float SH_4_7( in vec3 s ) { vec3 n = s.zxy; return k12 * n.x*n.z*(n.x*n.x-3.0*n.y*n.y); } float SH_4_8( in vec3 s ) { vec3 n = s.zxy; float x2 = n.x*n.x, y2 = n.y*n.y; return k17 * (x2 * (x2-3.0*y2) - y2*(3.0*x2 - y2)); } vec3 map( in vec3 p ) { vec3 p00 = p - vec3( 0.00, 3.0, 0.0); vec3 p01 = p - vec3(-1.25, 2.0, 0.0); vec3 p02 = p - vec3( 0.00, 2.0, 0.0); vec3 p03 = p - vec3( 1.25, 2.0, 0.0); vec3 p04 = p - vec3(-2.50, 0.5, 0.0); vec3 p05 = p - vec3(-1.25, 0.5, 0.0); vec3 p06 = p - vec3( 0.00, 0.5, 0.0); vec3 p07 = p - vec3( 1.25, 0.5, 0.0); vec3 p08 = p - vec3( 2.50, 0.5, 0.0); vec3 p09 = p - vec3(-3.75,-1.0, 0.0); vec3 p10 = p - vec3(-2.50,-1.0, 0.0); vec3 p11 = p - vec3(-1.25,-1.0, 0.0); vec3 p12 = p - vec3( 0.00,-1.0, 0.0); vec3 p13 = p - vec3( 1.25,-1.0, 0.0); vec3 p14 = p - vec3( 2.50,-1.0, 0.0); vec3 p15 = p - vec3( 3.75,-1.0, 0.0); vec3 p16 = p - vec3(-5.00,-2.7, 0.0); vec3 p17 = p - vec3(-3.75,-2.7, 0.0); vec3 p18 = p - vec3(-2.50,-2.7, 0.0); vec3 p19 = p - vec3(-1.25,-2.7, 0.0); vec3 p20 = p - vec3( 0.00,-2.7, 0.0); vec3 p21 = p - vec3( 1.25,-2.7, 0.0); vec3 p22 = p - vec3( 2.50,-2.7, 0.0); vec3 p23 = p - vec3( 3.75,-2.7, 0.0); vec3 p24 = p - vec3( 5.00,-2.7, 0.0); float r, d; vec3 n, s, res; #ifdef SHOW_SPHERES #define SHAPE (vec3(d-0.35, -1.0+2.0*clamp(0.5 + 16.0*r,0.0,1.0),d)) #else #define SHAPE (vec3(d-abs(r), sign(r),d)) #endif d=length(p00); n=p00/d; r = SH_0_0( n ); s = SHAPE; res = s; d=length(p01); n=p01/d; r = SH_1_0( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p02); n=p02/d; r = SH_1_1( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p03); n=p03/d; r = SH_1_2( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p04); n=p04/d; r = SH_2_0( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p05); n=p05/d; r = SH_2_1( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p06); n=p06/d; r = SH_2_2( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p07); n=p07/d; r = SH_2_3( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p08); n=p08/d; r = SH_2_4( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p09); n=p09/d; r = SH_3_0( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p10); n=p10/d; r = SH_3_1( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p11); n=p11/d; r = SH_3_2( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p12); n=p12/d; r = SH_3_3( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p13); n=p13/d; r = SH_3_4( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p14); n=p14/d; r = SH_3_5( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p15); n=p15/d; r = SH_3_6( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p16); n=p16/d; r = SH_4_0( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p17); n=p17/d; r = SH_4_1( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p18); n=p18/d; r = SH_4_2( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p19); n=p19/d; r = SH_4_3( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p20); n=p20/d; r = SH_4_4( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p21); n=p21/d; r = SH_4_5( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p22); n=p22/d; r = SH_4_6( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p23); n=p23/d; r = SH_4_7( n ); s = SHAPE; if( s.x<res.x ) res=s; d=length(p24); n=p24/d; r = SH_4_8( n ); s = SHAPE; if( s.x<res.x ) res=s; return vec3( res.x, 0.5+0.5*res.y, res.z ); } vec3 intersect( in vec3 ro, in vec3 rd ) { vec3 res = vec3(1e10,-1.0, 1.0); float maxd = 20.0; float h = 1.0; float t = 0.0; vec2 m = vec2(-1.0); for( int i=0; i<200; i++ ) { if( h<0.001||t>maxd ) break; vec3 res = map( ro+rd*t ); h = res.x; m = res.yz; t += h*0.3; } if( t<maxd && t<res.x ) res=vec3(t,m); return res; } vec3 calcNormal( in vec3 pos ) { vec3 eps = vec3(0.001,0.0,0.0); return normalize( vec3( map(pos+eps.xyy).x - map(pos-eps.xyy).x, map(pos+eps.yxy).x - map(pos-eps.yxy).x, map(pos+eps.yyx).x - map(pos-eps.yyx).x ) ); } void mainImage( out vec4 fragColor, in vec2 fragCoord ) { vec3 tot = vec3(0.0); for( int m=0; m<AA; m++ ) for( int n=0; n<AA; n++ ) { vec2 p = (-iResolution.xy + 2.0*(fragCoord.xy+vec2(float(m),float(n))/float(AA))) / iResolution.y; // camera float an = 0.314*TIME - 10.0*iMouse.x/iResolution.x; float dist = 8.0; vec3 ro = vec3(dist*sin(an),0.0,dist*cos(an)); vec3 ta = vec3(0.0,0.0,0.0); // camera matrix vec3 ww = normalize( ta - ro ); vec3 uu = normalize( cross(ww,vec3(0.0,1.0,0.0) ) ); vec3 vv = normalize( cross(uu,ww)); // create view ray vec3 rd = normalize( p.x*uu + p.y*vv + 2.0*ww ); // background vec3 col = vec3(0.5) * clamp(1.0-length(p)*0.5, 0.0, 1.0); // raymarch vec3 tmat = intersect(ro,rd); if( tmat.y>-0.5 ) { // geometry vec3 pos = ro + tmat.x*rd; vec3 nor = calcNormal(pos); vec3 ref = reflect( rd, nor ); // material vec3 mate = 0.5*mix( vec3(1.0,0.6,0.15), vec3(0.2,0.4,0.5), tmat.y ); float occ = clamp( 2.0*tmat.z, 0.0, 1.0 ); float sss = pow( clamp( 1.0 + dot(nor,rd), 0.0, 1.0 ), 1.0 ); // lights vec3 lin = 2.5*occ*vec3(1.0,1.00,1.00)*(0.6+0.4*nor.y); lin += 1.0*sss*vec3(1.0,0.95,0.70)*occ; // surface-light interacion col = mate.xyz * lin; } // gamma col = pow( clamp(col,0.0,1.0), vec3(0.4545) ); tot += col; } tot /= float(AA*AA); fragColor = vec4( tot, 1.0 ); } |
Applications: Irradiance
One can easily use the equation above and extract high order of spherical harmonics using prefilter.c, then apply them for irradiance
First, convert the coordinate systems from panorama:
x = 2.0 * i / width – 1.0;
y = 2.0 * j / width – 1.0;theta = y * PI;
phi = x * PI / 2.0;x = sin(phi) * sin(theta);
y = cos(theta);
z = sin(phi) * cos(theta);
The delta angle is given by
domega = (2 * PI / width) * (2 * PI / width) * sinc(theta);
Then one can update the integration using the fomulas
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 |
, float domega, float x, float y, float z)"] for (int col = 0; col < 3; col++) { float c; /* A different constant for each coefficient */ /* L_{00}. Note that Y_{00} = 0.282095 */ coeffs[0][col] += rgb[col] * k01 * domega; /* L_{1m}. -1 <= m <= 1. The linear terms */ coeffs[1][col] += rgb[col] * (k02 * y) * domega; /* Y_{1-1} = 0.488603 y */ coeffs[2][col] += rgb[col] * (k02 * z) * domega; /* Y_{10} = 0.488603 z */ coeffs[3][col] += rgb[col] * (k02 * x) * domega; /* Y_{11} = 0.488603 x */ /* The Quadratic terms, L_{2m} -2 <= m <= 2 */ /* First, L_{2-2}, L_{2-1}, L_{21} corresponding to xy,yz,xz */ c = 1.092548; coeffs[4][col] += rgb[col] * (k03 * xy) * domega; /* Y_{2-2} = 1.092548 xy */ coeffs[5][col] += rgb[col] * (k03 * yz) * domega; /* Y_{2-1} = 1.092548 yz */ /* L_{20}. Note that Y_{20} = 0.315392 (3z^2 - 1); c = 0.315392;*/ coeffs[6][col] += rgb[col] * (k04 * (3 * z2 - 1.0)) * domega; coeffs[7][col] += rgb[col] * (k03 * xz) * domega; /* Y_{21} = 1.092548 xz */ /* L_{22}. Note that Y_{22} = 0.546274 (x^2 - y^2); c = 0.546274; */ coeffs[8][col] += rgb[col] * (k05 * (x2 - y2)) * domega; //---------------------------------------------------------- coeffs[9][col] += rgb[col] * (k06 * y * (3.0 * x2 - y2)) * domega; coeffs[10][col] += rgb[col] * (k07 * xy * z) * domega; coeffs[11][col] += rgb[col] * (k08 * y * (5.0 * z2 - 1.0)) * domega; coeffs[12][col] += rgb[col] * (k09 * z * (5.0 * z2 - 3.0)) * domega; coeffs[13][col] += rgb[col] * (k08 * z * (5.0 * z2 - 1.0)) * domega; coeffs[14][col] += rgb[col] * (k10 * z * (x2 - y2)) * domega; coeffs[15][col] += rgb[col] * (k06 * x * (x2 - 3.0 * y2)) * domega; //---------------------------------------------------------- coeffs[16][col] += rgb[col] * (k11 * xy * (x2 - y2)) * domega; coeffs[17][col] += rgb[col] * (k12 * (3.0 * x2 - y2) * yz) * domega; coeffs[18][col] += rgb[col] * (k13 * xy * (7.0 * z2 - r2)) * domega; coeffs[19][col] += rgb[col] * (k14 * yz * (7.0 * z2 - 3.0 * r2)) * domega; coeffs[20][col] += rgb[col] * (k15 * (35.0 * z2 * z2 - 30.0 * z2 * r2 + 3.0 * r2 * r2)) * domega; coeffs[21][col] += rgb[col] * (k14 * xz * (7.0 * z2 - 3.0 * r2)) * domega; coeffs[22][col] += rgb[col] * (k16 * (x2 - y2) * (7.0 * z2 - r2)) * domega; coeffs[23][col] += rgb[col] * (k12 * (x2 - 3.0 * y2) * xz) * domega; coeffs[24][col] += rgb[col] * (k17 * (x2 * (x2 - 3.0 * y2) - y2 * (3.0 * x2 - y2))) * domega; } |
