/* Glass ----- This script uses algorithms similar to raytracing to render a glass bead over a checkered background. Author : Andreas Jonsson Created: April 21st, 2005 Updated: July 28th, 2005 Texture Generator version 1.1.0 */ // Snell's law: // n1*sin(a1) = n2*sin(a2) // where n1 is the refraction index for the source material // and n2 the index for the destination material. // Critical angle, where all incoming light is reflected: // a1 = asin(n2*sin(pi/2)/n1) // It can be seen that this critical angle only exists where the light is // travelling from a material with higher refraction index to a lower index. // Fresnel's law: // rs = (sin(a1-a2)/sin(a1+a2))^2 // rp = (tan(a1-a2)/tan(a1+a2))^2 // Where a1 is incoming angle, and a2 is the angle of the refracted ray // rs and rp are the fraction of the light that is reflected, where each // are polarized components, i.e. either parallel to the surface or perpendicular. // For rays perpendicular to the surface the formulae reduce to: // r = ((n2/n1 - 1)/(n2/n1 + 1))^2 // Some material refraction indices // air = 1.00 // water = 1.33 // glass = 1.52 // quartz = 1.54 // diamond = 2.41 const uint texture_start_size = 1024; const uint texture_end_size = 256; const float sphere_radius = 400; const float sphere_cx = 500; const float sphere_cy = 490; const float n_coeff = 1.52f; const float pi = 3.141592f; // Light direction const float light_vx = 0.4f; const float light_vy = 0.8f; const float light_vz = -1; void checkers(image @img) { img.Resize(16, 16); for( int y = 0; y < 16; ++y ) { for( int x = 0; x < 16; ++x ) { float val; if( bits(x + y)&1 == 1 ) val = 1; else val = 0.5; img.Pixel(x,y) = pixel(val, val*0.95, val*0.6); } } img.Resize(texture_start_size, texture_start_size); } void glass(image @img) { image img2 = img; int w = img.width; int h = img.height; for( int y = 0; y < h; ++y ) { for( int x = 0; x < w; ++x ) { float rad = sqrt((x-sphere_cx)*(x-sphere_cx) + (y-sphere_cy)*(y-sphere_cy)); if( rad < sphere_radius ) { // Compute the ingoing angle as the angle between the surface normal and ray (normal points inward) float a1 = asin(rad/sphere_radius); // Compute the outgoing angle with Snell's law float a2 = asin(sin(a1)/n_coeff); // Since incoming ray comes straight in we can determine // the angle between the outgoing ray and incoming ray by // subtracting the outgoing angle from the incoming angle. // Resulting ray float a3 = a1 - a2; // Determine at what height the incoming ray hit the surface float height = sin(pi/2-a1)*sphere_radius; // Determine at what distance from the center the ray hits the bottom float r = rad - sin(a3)*height; // Pick the pixel that is r distance from the center float x2, y2; if( rad > 0 ) { x2 = (x-sphere_cx)/rad * r; y2 = (y-sphere_cy)/rad * r; } else { x2 = x-sphere_cx; y2 = y-sphere_cy; } pixel p = img.Pixel(x2+sphere_cx,y2+sphere_cy); // Compute the Fresnel term to determine how much light goes into the glass float f = 1 - fresnel(a1,a2); // Multiply with the Fresnel term to find how much of the color is really seen p = pixel(p.r*f, p.g*f, p.b*f); img2.Pixel(x,y) = p; // TODO: We should also add light that was reflected from outside the glass sphere } } } img = img2; } float fresnel(float a1, float a2) { float rs; // Compute the Fresnel term // NOTE: We are actually cheating here by only computing one of the terms if( a1+a2 != 0 ) rs = sin(a1-a2)/sin(a1+a2); else // For rays perpendicular to the surface the formula reduce to the following rs = (n_coeff-1)/(n_coeff+1); return rs*rs; } float getAngle(float v1x, float v1y, float v1z, float v2x, float v2y, float v2z) { // Normalize the two vectors float m = 1/sqrt(v1x*v1x + v1y*v1y + v1z*v1z); v1x *= m; v1y *= m; v1z *= m; m = 1/sqrt(v2x*v2x + v2y*v2y + v2z*v2z); v2x *= m; v2y *= m; v2z *= m; // Compute the dot product float d = v1x*v2x + v1y*v2y + v1z*v2z; // Return the angle between the vectors return acos(d); } void refractVector(float ix, float iy, float iz, float nx, float ny, float nz, float &out ox, float &out oy, float &out oz, float &out a1, float &out a2) { // Normalize the two vectors float m = 1/sqrt(ix*ix + iy*iy + iz*iz); ix *= m; iy *= m; iz *= m; m = 1/sqrt(nx*nx + ny*ny + nz*nz); nx *= m; ny *= m; nz *= m; // Determine the ingoing angle a1 = acos(ix*nx + iy*ny + iz*nz); if( a1 > 0 ) { // Compute the outgoing angle with Snell's law a2 = asin(sin(a1)/n_coeff); // Compute the direction of the refracted ray. It lies in // the same plane as the ingoing ray and surface normal. float X = sin(a2)/cos(a2); ox = ix - nx; oy = iy - ny; oz = iz - nz; m = X/sqrt(ox*ox + oy*oy + oz*oz); ox = ox*m + nx; oy = oy*m + ny; oz = oz*m + nz; } else { ox = ix; oy = iy; oz = iz; } } void light(image @img) { randomizer rnd; const float num_rays = 4; const float intensity = 1/(num_rays*1.5); int w = img.width; int h = img.height; image img2(w,h); float vx0 = light_vx; float vy0 = light_vy; float vz0 = light_vz; // Initial position outside the sphere float px0 = -vx0*sphere_radius; float py0 = -vy0*sphere_radius; float pz0 = -vz0*sphere_radius; // Normalize the vector float m = 1/sqrt(vx0*vx0 + vy0*vy0 + vz0*vz0); vx0 *= m; vy0 *= m; vz0 *= m; for( int y = 0; y < h; ++y ) { for( int x = 0; x < w; ++x ) { for( int i = 0; i < num_rays; ++i ) { // The parallel light comes in sligthly from above // Formulae: R(t) = P + t*V float vx = vx0, vy = vy0, vz = vz0; float px = x + px0 + rnd.getNumber(); float py = y + py0 + rnd.getNumber(); float pz = pz0; // The sphere is located at (cx, cy, 0) // Formulae: (x-cx)^2 + (y-cy)^2 + (z-cz)^2 = radius^2 // Solve the linear equation for t to determine if the ray // intersects the sphere // x = x0 + xd * t P = (x0, y0, z0) // y = y0 + yd * t V = (xd, yd, zd) // z = z0 + zd * t // (x0 + xd*t - cx)^2 + (y0 + yd*t - cy)^2 + (z0 + zd*t - cz)^2 = radius^2 // A*t^2 + B*t + C = 0 // A = xd^2 + yd^2 + zd^2 // B = 2*(xd * (x0 - cx) + yd * (y0 - cy) + zd * (z0 - cz)) // C = (x0 - cx)^2 + (y0 - cy)^2 + (z0 - cz)^2 - radius^2 // If the ray vector is normalized then A = 1 // Solution: t0, t1 = 0.5*(-B (+|-) sqrt(B^2 - 4*C)) float rx = px-sphere_cx; float ry = py-sphere_cy; float rz = pz; float B = 2*(vx*rx + vy*ry + vz*rz); float C = rx*rx + ry*ry + rz*rz - sphere_radius*sphere_radius; float D = B*B - 4*C; if( D >= 0 ) { // Intersection occurs somewhere along the path // Verify if it intersects the sphere before the background plane // NOTE: We only need to verify one of the possible t values, // since we know which one is relevant float t = 0.5*(-B - sqrt(D)); pz = pz + t*vz; if( pz > 0 ) { // Intersection is above the plane px = px + t*vx; py = py + t*vy; // Surface normal (pointing into the sphere) float nx = sphere_cx - px; float ny = sphere_cy - py; float nz = -pz; // Refract the light vector to see where it hits the plane float a1, a2; refractVector(vx,vy,vz, // ingoing vector nx,ny,nz, // normal vx,vy,vz, // outgoing vector a1,a2); // angles // Determine where the refracted ray hits the surface float t = -pz/vz; px = px + vx*t; py = py + vy*t; // Use Fresnel's law to determine the fraction of // the light intensity that hits the surface. img2.Alpha(px,py) += intensity * (1 - fresnel(a1, a2)); } else { // Intersection is below the plane img2.Alpha(x,y) += intensity; } } else { img2.Alpha(x,y) += intensity; } } } } boxBlurAlpha(@img2, 3, 3); for( int y = 0; y < h; ++y ) { for( int x = 0; x < w; ++x ) { pixel p = img.Pixel(x,y); float a = img2.Alpha(x,y); img.Pixel(x,y) = pixel(p.r*a, p.g*a, p.b*a); } } } void boxBlurAlpha(image @img, int bw, int bh) { int x, y; int w = img.width; int h = img.height; //-------------- // smooth x float[] line(w); y = 0; while( y < h ) { // Copy original pixels to the line buffer x = 0; while( x < w ) { line[x] = img.Alpha(x, y); ++x; } // Compute sum of first bw pixels x = 0; float sum = 0.0; while( x < bw ) { sum += line[x]; ++x; } // Set blurred pixels x = 0; while( x < w ) { img.Alpha((x+bw/2)%w, y) = sum/float(bw); sum -= line[x]; sum += line[(x+bw)%w]; ++x; } ++y; } //--------------- // smooth y line.resize(h); x = 0; while( x < w ) { // Copy original pixels to the line buffer y = 0; while( y < h ) { line[y] = img.Alpha(x,y); ++y; } // Compute sum of first bw pixels y = 0; float sum = 0.0; while( y < bh ) { sum += line[y]; ++y; } // Set blurred pixels y = 0; while( y < h ) { img.Alpha(x, (y+bh/2)%h) = sum/float(bh); sum -= line[y]; sum += line[(y+bh)%h]; ++y; } ++x; } } void reflection(image @img) { image img2 = img; int w = img.width; int h = img.height; for( int y = 0; y < h; ++y ) { for( int x = 0; x < w; ++x ) { float rad = sqrt((x-sphere_cx)*(x-sphere_cx) + (y-sphere_cy)*(y-sphere_cy)); if( rad < sphere_radius ) { // Compute the ingoing angle as the angle between the surface normal and ray (normal points inward) float a1 = asin(rad/sphere_radius); // Determine at what height the incoming ray hit the surface float height = sin(pi/2-a1)*sphere_radius; // Compute the reflected vector by adding two normals to it float nx = x-sphere_cx; float ny = y-sphere_cy; float nz = height; float m = 1/sqrt(nx*nx + ny*ny + nz*nz); nx *= m; ny *= m; nz *= m; float d = -1*nz; nx = 0+2*nx; ny = 0+2*ny; nz = -1+2*nz; // Compute the angle between the reflected ray and the incoming light vector float a = getAngle(-nx, -ny, -nz, light_vx, light_vy, light_vz); if( a < 0.5f ) { a = 0.5f*cos(a/0.5f*pi)+0.5f; pixel p = img2.Pixel(x,y); img2.Pixel(x,y) = pixel(p.r+a,p.g+a,p.b+a); } } } } img = img2; } void main() { image img; // Generate a simple background image checkers(@img); // A light source would cast shadows on the background, and the // glass sphere would concentrate the light in a focal point light(@img); // Compute how the glass sphere refracts the light from the background glass(@img); // Add the direct reflected light reflection(@img); img.Resize(texture_end_size, texture_end_size); img.Show(); }