did you solve it already? if not here is another way, a bit more verbose then the wikipedia version:
what happens is:
- intersect: sphere1/sphere2 -> returns intersection plane1
- intersect: sphere1/sphere3 -> returns intersection plane2
- intersect: plane1/plane2 -> returns intersection ray
- intersect: ray/sphere1 (or any other) -> returns intersection points
// diewald, 12.08.2013 float[] sphere_1 = { 100, 150, 100, 100 }; // yellow, sphere[px,py,pz,r] float[] sphere_2 = { 50, 40, 50, 170 }; // green, sphere[px,py,pz,r] float[] sphere_3 = { 20, 200, 150, 80 }; // blue, sphere[px,py,pz,r] float cam_rx=0, cam_rz=0; public void setup() { size(800, 800, P3D); camera( 0.0f ,0.0f, 700.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f ); textFont(createFont("Monospaced.plain", 12)); textMode(SCREEN); // processing 1.5.1 } void updateCam(){ if( mousePressed ){ cam_rz += (pmouseX-mouseX)*0.01f; cam_rx += (pmouseY-mouseY)*0.01f; } rotateX(cam_rx); rotateZ(cam_rz); // default settings float fov = PI/3.0f; float cameraZ = (height/2.0f) / tan(fov/2.0f); perspective(fov, float(width)/float(height), cameraZ/10.0f, cameraZ*10.0f); } public void draw() { background(255); updateCam(); lights(); sphere_1[0] = map(mouseX, 0, width, -200, 200); strokeWeight(2); gizmo(255); // ... draw spheres ... noStroke(); fill(255,255,0,255); drawSphere(sphere_1); noStroke(); fill(0,255,0,255); drawSphere(sphere_2); noStroke(); fill(0,0,255,255); drawSphere(sphere_3); // ... compute intersections ... boolean isec_sphere_sphere_1 = false; boolean isec_sphere_sphere_2 = false; boolean isec_plane_plane = false; boolean isec_ray_sphere = false; // 1a) intersect: sphere1 - sphere2 -> plane_s1s2 float[] plane_s1s2 = new float[4]; // plane[nx,ny,nz,d] isec_sphere_sphere_1 = intersectSphereSphere(sphere_1, sphere_2, plane_s1s2); // 1b)intersect: sphere1 - sphere3 -> plane_s1s3 float[] plane_s1s3 = new float[4]; // plane[nx,ny,nz,d] isec_sphere_sphere_2 = intersectSphereSphere(sphere_1, sphere_3, plane_s1s3); // 2) intersect: plane_s1s2 - plane_s1s3 -> ray float[] plane_isec_pnt = new float[3]; // point[px,py,pz] float[] plane_isec_dir = new float[3]; // vector[x,y,z] isec_plane_plane = intersectPlanePlane(plane_s1s2, plane_s1s3, plane_isec_pnt, plane_isec_dir); // 3) intersect: ray - sphere1 -> 2 points float[][] points = new float[2][3]; isec_ray_sphere = intersectRaySphere(plane_isec_pnt, plane_isec_dir, sphere_1, points); // ... draw intersection results ... if( isec_sphere_sphere_1){ fill(150, 100); noStroke(); drawPlane(plane_s1s2); } if( isec_sphere_sphere_2){ fill(150, 100); noStroke(); drawPlane(plane_s1s3); } if(isec_plane_plane){ stroke(255,0,0); strokeWeight(3); drawRay(plane_isec_pnt, plane_isec_dir); } if( isec_ray_sphere ){ fill(0,255,255); noStroke(); drawPoint(points[0]); drawPoint(points[1]); } pushMatrix(); { resetMatrix(); translate(-width/2,-height/2); ortho(0, width, 0, height, 0,1); // difference in 1.5.1 and 2.xx ??? String info = ""; info += String.format("isec_sphere_sphere_1: %b\n", isec_sphere_sphere_1); info += String.format("isec_sphere_sphere_2: %b\n", isec_sphere_sphere_2); info += String.format("isec_plane_plane: %b\n", isec_plane_plane ); info += String.format("isec_ray_sphere: %b\n", isec_ray_sphere ); fill(0); text(info, 20, 20); } popMatrix(); } void gizmo(float s){ stroke(255,0,0); line(0,0,0,s,0,0); stroke(0,255,0); line(0,0,0,0,s,0); stroke(0,0,255); line(0,0,0,0,0,s); } void drawPoint(float[] p){ pushMatrix(); { translate(p[0], p[1], p[2]); sphere(4); } popMatrix(); } void drawRay(float[] plane_isec_pnt, float[] plane_isec_dir){ pushMatrix(); { float[] d = plane_isec_dir; float[] p = plane_isec_pnt; translate(p[0], p[1], p[2]); scale(400); line(0,0,0, +d[0], +d[1], +d[2]); line(0,0,0, -d[0], -d[1], -d[2]); } popMatrix(); } void drawPlane(float[] plane){ pushMatrix(); { applyMatrix( orthoNormalBasis(plane) ); rectMode(CENTER); rect(0,0, 700, 700); } popMatrix(); } void drawSphere(float[] sphere){ pushMatrix(); { translate(sphere[0], sphere[1], sphere[2]); sphere(sphere[3]); } popMatrix(); } // ... intersections ... boolean intersectSphereSphere(float[] s1, float[] s2, float[] isec_plane_out){ float[] dis = sub(s1,s2); float dis_len_sq = lenSq(dis); if( dis_len_sq == 0 ){ // spheres have same center. // so if they have a different radius, there's no intersection return false; } float dis_len = (float)Math.sqrt(dis_len_sq); float sum_radi = s1[3] + s2[3]; if( dis_len > sum_radi){ // no intersection, spheres are too far apart. return false; } else if (dis_len == sum_radi){ //spheres intersect at exactly one point } // plane normal (normalized) float[] normal = mult(dis, 1.0f/dis_len); // intersection distances from sphere1-center to plane float isec_d = (s1[3]*s1[3] - s2[3]*s2[3] + dis_len_sq)/(2f*dis_len); // intersection point: on "dis" and on plane float[] isec_p = add(s1, mult(normal, -isec_d)); // intersection plane float[] plane = {normal[0], normal[1], normal[2], 0}; plane[3] = dot(plane, isec_p); copy(normal, isec_plane_out); isec_plane_out[3] = dot(plane, isec_p); return true; } boolean intersectPlanePlane(float[] p1, float[] p2, float[] isec_pnt_out, float[] isec_dir_out){ float[] isec_dir = cross(p1,p2); float denom = dot(isec_dir,isec_dir); if( denom < 0.00001 ) return false; float[] tmp1 = mult(p2, p1[3]); float[] tmp2 = mult(p1, p2[3]); float[] tmp3 = sub(tmp1, tmp2); float[] isec_pnt = cross(tmp3,isec_dir); isec_pnt = mult(isec_pnt, 1f/denom); isec_dir = normalize(isec_dir); // return normalized direction copy(isec_dir, isec_dir_out); copy(isec_pnt, isec_pnt_out); return true; } boolean intersectRaySphere(float[] ray_p, float[] ray_d, float[] sphere, float[][] isec_points_out){ float[] m = sub(ray_p, sphere); float b = dot(m,ray_d); // ray_d must be normalized here float c = dot(m,m) - sphere[3]*sphere[3]; // if( c > 0.0f && b > 0.0f) return false; float discr = b*b - c; if( discr < 0) return false; discr = (float)Math.sqrt(discr); float t1 = -b + discr; float t2 = -b - discr; isec_points_out[0] = add(ray_p, mult(ray_d, t1)); isec_points_out[1] = add(ray_p, mult(ray_d, t2)); return true; } // ... vector functions ... void copy(float[] src, float[] dst){ dst[0] = src[0]; dst[1] = src[1]; dst[2] = src[2]; } float[] negate(float[] a){ return new float[]{-a[0], -a[1], -a[2] }; } float[] add(float[] a, float[] b){ return new float[]{a[0]+b[0], a[1]+b[1], a[2]+b[2] }; } float[] sub(float[] a, float[] b){ return new float[]{a[0]-b[0], a[1]-b[1], a[2]-b[2] }; } float[] mult(float[] a, float val){ return new float[]{a[0]*val, a[1]*val, a[2]*val }; } float lenSq(float[] a){ return dot(a,a); } float dot(float[] a, float[] b){ return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; } float[] cross(float[] a, float[] b) { return new float[]{a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]}; } float[] normalize(float[] a) { return mult(a, (float)(1f/Math.sqrt(lenSq(a)))); } // create a new matrix, having the plane normal as new z-axis PMatrix3D orthoNormalBasis(float[] plane) { float[] up = {0, 1, 0}; float[] z = plane; float[] x = cross(up, z); x = normalize(x); float[] y = cross( z, x); y = normalize(y); PMatrix3D mat = new PMatrix3D(); mat.m00 = x[0]; mat.m10 = x[1]; mat.m20 = x[2]; mat.m01 = y[0]; mat.m11 = y[1]; mat.m21 = y[2]; mat.m02 = z[0]; mat.m12 = z[1]; mat.m22 = z[2]; mat.translate(0,0,plane[3]); return mat; };
