30:T67a,viz={ const start=new Date(); const colors = d3[colorScheme]; const {points,intersects}=geometry; const twopi=2*Math.PI; const wrap = angle=>(angle+twopi)%twopi; const [xlo,xhi] = [d3.min(circles, C=>C.x-C.r),d3.max(circles, C=>C.x+C.r)]; const [ylo,yhi] = [d3.min(circles, C=>C.y-C.r),d3.max(circles, C=>C.y+C.r)]; const aspect = (yhi - ylo)/(xhi - xlo); const w=width; const h=Math.round(w*aspect); const X=d3.scaleLinear().domain([xlo,xhi]).range([0,w]); const Y=d3.scaleLinear().domain([ylo,yhi]).range([h,0]); const scale=X(1)-X(0); const c=DOM.context2d(w,h); c.textAlign = "center"; c.textBaseline = "middle"; c.strokeStyle = "skyblue"; c.fillStyle = "skyblue"; regions.regions.forEach((R,j) => { const n=R.inside.size; c.beginPath(); R.edges.forEach(({from,to,dir},k) => { const [xc,yc] = [X(from.self.x),Y(from.self.y)]; let [ang1,ang2] = [wrap(-from.angle), wrap(-to.angle)]; if(ang1==ang2) ang2=ang1-twopi; // isolated circle c.arc(xc, yc, scale*from.self.r, ang1, ang2, dir>0); }); c.fillStyle = n==0? "lightgray" : (colors[(n-1) % colors.length]); c.fill(); }); c.strokeStyle = c.fillStyle = "#000"; circles.forEach((C,i) => { c.beginPath(); c.arc(X(C.x), Y(C.y), scale*C.r, 0, 2*Math.PI); c.stroke(); //c.fillText(`${i}`, X(C.x), Y(C.y)); }); //points.forEach(P => c.fillText(`${P.idx}`, X(P.pt[0]), Y(P.pt[1]))); const elapsed=new Date() - start; c.fillStyle="#000"; c.textAlign="left"; c.textBaseline="bottom"; c.fillText(`${geometry.elapsed} + ${regions.elapsed} + ${elapsed} ms`, 2, h-2); return c.canvas; }31:T487,drawCircles=function(circles, { w=width, inset=4, colorScheme="Category10", pointFill=null, pointRadius=4 }={}) { const {PI}=Math; // Compute the geometry. const {points,intersects}=getGeometry(circles) // Fit a canvas to the circle bounds. const [xlo,xhi] = [d3.min(circles, C=>C.x-C.r),d3.max(circles, C=>C.x+C.r)]; const [ylo,yhi] = [d3.min(circles, C=>C.y-C.r),d3.max(circles, C=>C.y+C.r)]; const aspect = (yhi - ylo)/(xhi - xlo); const h=Math.round(w*aspect); const X=d3.scaleLinear().domain([xlo,xhi]).range([inset,w-inset]); const Y=d3.scaleLinear().domain([ylo,yhi]).range([h-inset,inset]); const scale=X(1)-X(0); console.log(X(1)-X(0), Y(0)-Y(1)); const c=DOM.context2d(w,h); // Draw circles. c.strokeStyle = c.fillStyle = "#000"; circles.forEach((C,i) => { c.beginPath(); c.arc(X(C.x), Y(C.y), scale*C.r, 0, 2*Math.PI); c.stroke(); //c.fillText(`${i}`, X(C.x), Y(C.y)); }); if(pointFill) { // Draw points. c.fillStyle = pointFill; points.forEach(P => { c.beginPath(); c.arc(X(P.pt[0]), Y(P.pt[1]), pointRadius, 0, 2*PI); c.fill(); }); } return c.canvas; }32:T1b5e,getGeometry=function(circles, {rMinAbs=1e-10,rMinRel=1e-10,discEps=1e-8}={}) { const start=new Date(); const {min,max,abs,sqrt,cos,sin,atan2,PI} = Math; const [rmin,rmax]=d3.extent(circles, C=>C.r); if((rmin(angle+twopi)%twopi; const n=circles.length; const intersects=Array.from({length:n}, ()=>[]), points=[], nests=new Map(); const membership=Array.from({length:n}, (_,i)=>i); for(let i=1; i +1 : the circles are non overlapping // Below, we merge [1] into [0] and [3] into [4]. const disc = (d2 - ri2 - rj2) / (2*ri*rj); //console.log("disc",i,j,disc); if(disc <= -1+discEps) { // circles are nested (and also tangent if disc==-1) const [inner,outer] = (ri > rj) ? [j,i] : [i,j]; console.log(`circle ${inner} is inside ${outer}`); if(nests.has(inner)) nests.get(inner).push(outer); else nests.set(inner, [outer]); continue; } if(disc >= +1-discEps) { // circles are non overlapping (and also tangent if disc==+1) continue; } // If we get here, the circles have some overlap with 2 distinct intersection points. membership[i] = membership[j]; const d=sqrt(d2); const x = (d2 - rj2 + ri2) / (2*d); const a = sqrt(max(0,ri2 - x**2)); // Calculate angle of ray from center of i to center of j. const alpha = atan2(dy, dx); // Calculate half angle between intersection points relative to center of i. const beta = atan2(a, x); // Calculate half angle between intersection points relative to center of j. const gamma = atan2(a, d-x); // Check if this intersection is almost tangent. if(min(beta, gamma) < 1e-7) { // This is probably no longer needed...does it ever happen? console.log(`circles(${i},${j}) are (almost) tangent`); // Average the rays from each circle center. const [C,S]=[cos(alpha),sin(alpha)]; const xij=[xi+ri*C,xj-rj*C]; const yij=[yi+ri*S,yj-rj*S]; const pavg={idx:[i,j], pt:[(xij[0]+xij[1])/2,(yij[0]+yij[1])/2]}; const Ii = {angle:wrap(alpha), self:circles[i], p12:pavg}; const Ij = {angle:wrap(alpha+PI), self:circles[j], p12:pavg}; // Add (circular) self references that we will remove after sorting, to avoid memory leaks. Ii.other= {circle:i, idx:Ii}; Ij.other= {circle:j, idx:Ij}; // Collapse the two points down to one. points.push(pavg); intersects[i].push(Ii); intersects[j].push(Ij); } else { // Calculate and save intersection coordinates. const phi=[alpha-beta,alpha+beta]; const xij = phi.map(p=>xi+ri*cos(p)); const yij = phi.map(p=>yi+ri*sin(p)); const pij={idx:[i,j], pt:[xij[0],yij[0]]}; const pji={idx:[j,i], pt:[xij[1],yij[1]]} points.push(pij, pji); // Save intersection angles relative to each circle center. const Ii1 = {angle:wrap(alpha-beta), self:circles[i], p12:pij}; const Ii2 = {angle:wrap(alpha+beta), self:circles[i], p12:pji}; const Ij1 = {angle:wrap(alpha+gamma+PI), self:circles[j], p12:pij}; const Ij2 = {angle:wrap(alpha-gamma+PI), self:circles[j], p12:pji}; // Add (circular) cross references that we will remove after sorting, to avoid memory leaks. Ii1.other = {circle:j, idx:Ij1}; Ij1.other = {circle:i, idx:Ii1}; Ii2.other = {circle:j, idx:Ij2}; Ij2.other = {circle:i, idx:Ii2}; intersects[i].push(Ii1, Ii2); intersects[j].push(Ij1, Ij2); } } } // Sort intersections along each circle. intersects.forEach(I => { I.sort((A,B)=>A.angle-B.angle); // sort in place for this circle I.forEach((A,i)=>A.idx=i); // record our sorted index along this circle }); // Remove circular references with a second pass. intersects.forEach(I => { I.forEach(A => A.other.idx=A.other.idx.idx); }); // Assign each circle to a group of non-intersecting circles. const edges=intersects.map(I => new Set(I.map(({other}) => other.circle))); const belongs=new Map(); const grps=Array.from({length:circles.length},()=>null); let grpIdx=null; const applyLabel = (grpIdx,circleIdx) => { if(grps[circleIdx]==null) { grps[circleIdx] = grpIdx; belongs.get(grpIdx).push(circleIdx); edges[circleIdx].forEach(neighborIdx => applyLabel(grpIdx, neighborIdx)); } }; grps.forEach((gid,circleIdx) => { if(gid==null) { grpIdx=belongs.size; belongs.set(grpIdx, []); applyLabel(grpIdx,circleIdx); } }); edges.forEach((neighbors,idx) => { if(grps[idx]==null) { // Start a new group. grpIdx = belongs.length; belongs.push([idx]); } }); // Handle isolated circles. intersects.forEach((I,i) => { if(I.length==0) { const {x,y,r}=circles[i]; const pt={idx:[i,i], pt:[x+r,y]}; points.push(pt); I.push({angle:0, self:circles[i], p12:pt, other:{circle:i, idx:0}}); } }); // Build an inverse of the nesting map. const contains=new Map(); nests.forEach((outers,inner)=>{ outers.forEach(outer=>{ const inners=contains?.get(outer) ?? []; inners.push(inner); contains.set(outer, inners); }); }); // Calculate arc areas (zero for non-intersected circles). intersects.forEach((I,circleIdx) => { const rsqby2=circles[circleIdx].r**2/2; const phi=I.map(edge=>edge.angle); const last=phi.length-1; I.forEach((edge,idx) => { edge.dphi=(idx==last) ? phi[0]-phi[last]+twopi : phi[idx+1]-phi[idx]; edge.arcArea = rsqby2*(edge.dphi - sin(edge.dphi)); }); }); // Invert membership to identify the circles within each non-intersecting group. const groups=new Map(); membership.forEach((groupIdx,circleIdx) => { if(groups.has(groupIdx)) groups.get(groupIdx).push(circleIdx); else groups.set(groupIdx, [circleIdx]); }); const elapsed=new Date() - start; return {points:points, intersects:intersects, nests:nests, contains:contains, //groups:groups, membership:membership, edges:edges, grps:grps, belongs:belongs}; membership:grps, groups:belongs, elapsed:elapsed}; }33:T1569,getRegions = function(geometry) { const start=new Date(); const {intersects,nests,contains,membership} = geometry; const regions=[]; // Track which arc edges have been used to create a region so far. Each // edge should be used exactly once in each (ccw/cw) direction. const used=intersects.map(I=>Array.from({length:I.length},()=>[false,false])); // Loop over all intersections to generate all possible regions. [+1,-1].forEach(dir0 => { intersects.forEach((I,circleIdx0) => { I.forEach((A,intersectIdx0)=>{ if(!used[circleIdx0][intersectIdx0][(1-dir0)/2]) { //console.log("start", circleIdx0, intersectIdx0, dir0); const region={ label:"", edges:[], inside:new Set(), outside:new Set(), maybe:new Set(), area:0, membership:membership[circleIdx0], }; let circleIdx=circleIdx0, intersectIdx=intersectIdx0, dir=dir0; // +1=ccw, -1=cw const vertices=[]; while(region.edges.length < 10000) { // protect against a bug leading to an inf loop const edge={from:intersects[circleIdx][intersectIdx], dir:dir}; region.label += "("+[circleIdx,intersectIdx,dir]+")"; vertices.push(intersects[circleIdx][intersectIdx].p12.pt); const nint=intersects[circleIdx].length; const inside=contains.get(circleIdx) ?? []; if(dir>0) { // We are inside this circle since we are following this edge ccw. region.inside.add(circleIdx); // We are also inside any circles this circle is nested within. if(nests.has(circleIdx)) nests.get(circleIdx).forEach(idx=>region.inside.add(idx)); } else { // We are outside this circle since we are following this edge cw. region.outside.add(circleIdx); } // We *might* be inside any circles this circle intersects. const others=intersects[circleIdx].map(({other})=>other.circle); others.forEach(idx => region.maybe.add(idx)); // Record that we have used this edge in this direction. used[circleIdx][intersectIdx][(1-dir)/2]=true; // Update the net arc area contribution to this region's total area. // When dir=-1, the area is associated with the previous intersection. const arcIdx=(intersectIdx-(1-dir)/2+nint) % nint; //console.log("arc", circleIdx, arcIdx, "dir", dir, "dphi", intersects[circleIdx][arcIdx].dphi); region.area += dir*intersects[circleIdx][arcIdx].arcArea; // Travel along circleIdx from intersectIdx. intersectIdx=(intersectIdx+dir+nint) % nint; edge.to=intersects[circleIdx][intersectIdx]; // move to the intersecting circle. const other=intersects[circleIdx][intersectIdx].other; [circleIdx, intersectIdx] = [other.circle,other.idx]; region.edges.push(edge); // back where we started? //console.log("next", circleIdx, intersectIdx); if((circleIdx==circleIdx0) && (intersectIdx==intersectIdx0)) break; // what direction do we take next to keep turning "left"? // Get current location. const [vx,vy]=intersects[circleIdx][intersectIdx].p12.pt; // Get vector to center of circle we got here on. const [cx1,cy1]=[edge.from.self.x-vx,edge.from.self.y-vy]; // Get center of circle we will leave on. const [cx2,cy2]=[circles[circleIdx].x-vx,circles[circleIdx].y-vy]; // Calculate cross product. const xprod = cx1*cy2 - cx2*cy1; dir *= Math.sign(xprod); } // Resolve the "maybe" cases now. const [vx,vy]=region.edges[0].from.p12.pt; [...region.maybe].filter(idx => !region.inside.has(idx) && !region.outside.has(idx)) .forEach(idx => { const {x,y,r}=intersects[idx][0].self; if((vx-x)**2 + (vy-y)**2 < r*r) region.inside.add(idx); }); delete region.outside; delete region.maybe; // Add the polygon area to this region's total area. region.area += polygonArea(...vertices); regions.push(region); } }); //circleIdx0 }); //intersectIdx0 }); //dir0 // sort regions in order of increasing number of overlaps. regions.sort((R1,R2)=>R1.inside.size-R2.inside.size); // Calculate the total area for each number of overlaps. //const maxOverlaps = regions[regions.length-1].inside.size; const maxOverlaps=d3.max(regions, R=>R.inside.size); const areas=d3.range(1, maxOverlaps+1).map( n => d3.sum(regions.filter(R=>R.inside.size==n), R=>R.area)); const groups=new Set(regions.map(R=>R.membership)); groups.forEach(gid => { // Select interior regions belonging to this group. const groupRegions=regions.filter(R=>(R.membership==gid)&&(R.inside.size>=1)); const embedLevel=d3.min(groupRegions, R=>R.inside.size); if(embedLevel > 1) { // This non-intersecting group is embedded within another group, so its area // has been double counted. Correct for that in the areas array here. const groupArea = d3.sum(groupRegions, R=>R.area); areas[embedLevel - 2] -= groupArea; } }); const elapsed=new Date() - start; return {regions:regions, areas:areas, elapsed:elapsed}; }34:[] 2e:{"title":"Circle Overlaps","cells":[{"id":0,"value":"# Circle Overlaps","pinned":false,"hidden":false,"mode":"md"},{"id":1301,"value":"md`How do overlapping circles divide up a flat space? First, define some circles:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1372,"value":"C=[\n // Each circle is specified by its center (x,y) and radius r.\n {x:0,y:0,r:10}, {x:-10,y:10,r:15}, {x:5,y:5,r:12}, {x:25,y:10,r:5},\n {x:45,y:10,r:12}, {x:40,y:5,r:3}, {x:-23,y:10,r:2}, {x:25,y:0,r:5}\n]","pinned":true,"hidden":false,"mode":"ojs"},{"id":1303,"value":"drawCircles(C)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1352,"value":"md`Each distinct pair of circles must be either:\n - non overlapping, or\n - fully nested, with no points of intersection, or\n - overlapping at:\n - two distinct points, or\n - a single point, where the circles are tangent, with:\n - one circle fully nested within the other, or\n - no overlapping area between the circles.\n\nEach of these cases occurs in the example above.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1380,"value":"drawCircles(C, {pointFill:\"red\"})","pinned":true,"hidden":false,"mode":"ojs"},{"id":446,"value":"md`**TODO:**\n- allow circle accessor fns\n- save circle x,y,r in intersects\n- remove refs to circles in getRegions\n- does \\`nests.has(circleIdx)\\` logic apply to both dirs?\n- remove any unused logic from getGeometry\n- generalize to circles on sphere\n- publish as ES6 module on npm?\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1189,"value":"circles3=[\n {x:0,y:0,r:1},{x:0,y:0,r:1.5}\n]","pinned":false,"hidden":false,"mode":"ojs"},{"id":92,"value":"circles2={\n const {cos,sin,PI}=Math;\n const r=1;\n const n=4;\n const s=1.3;\n const phi=d3.range(n).map(k=>2*PI*k/n);\n const circles=phi.map(phi=>{\n return {x:s*r*cos(phi), y:s*r*sin(phi), r:r};\n });\n circles.push({x:s*r,y:s*r,r:r/5});\n circles.push({x:0,y:0,r:r*2.4});\n return circles;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1075,"value":"circles4={\n const rng=d3.randomUniform; //.source(d3.randomLcg(123));\n const [xgen,ygen,rgen] = [[-5,5], [-3,3], [0.1,0.3]].map(([lo,hi])=>rng(lo,hi));\n return Array.from({length:1000}, ()=>({x:xgen(), y:ygen(), r:rgen()}));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1404,"value":"circles=C","pinned":true,"hidden":false,"mode":"ojs"},{"id":861,"value":"md`Area cross check: ${tex.block`\n\\sum_i \\pi r_i^2 = ${Number(circleArea).toFixed(5)}\n\\quad, \\quad\n\\sum_k n_k A_k = ${Number(regionArea).toFixed(5)}\n`}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":11,"value":"$30","pinned":true,"hidden":false,"mode":"ojs"},{"id":1309,"value":"$31","pinned":true,"hidden":false,"mode":"ojs"},{"id":572,"value":"{\n const maxOverlaps=d3.max(regions.regions, R=>R.inside.size);\n return colorbar(d3.range(1, maxOverlaps+1).map(n=>`n=${n}`), {w:width});\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":688,"value":"circleArea=d3.sum(circles,C=>Math.PI*C.r**2)","pinned":false,"hidden":false,"mode":"ojs"},{"id":767,"value":"regionArea = d3.sum(regions.areas, (A,n)=>A*(n+1))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1157,"value":"totalArea = d3.sum(regions.regions, R=>Math.max(0,R.area));","pinned":false,"hidden":false,"mode":"ojs"},{"id":467,"value":"regions = getRegions(geometry)","pinned":false,"hidden":false,"mode":"ojs"},{"id":21,"value":"geometry = getGeometry(circles)","pinned":false,"hidden":false,"mode":"ojs"},{"id":19,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":462,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":638,"value":"polygonArea = (...V) => {\n // Using https://en.wikipedia.org/wiki/Shoelace_formula\n // Note that the calculated area is signed, with values >0 for ccw non-intersecting edges.\n // When edges intersect, the result is the net (ccw) - (cw) area.\n const n=V.length;\n const pair=(i,j)=>(V[i][0]+V[j][0])*(V[i][1]-V[j][1]);\n let area=pair(0,n-1);\n for(let i=0; iMath.round(i*w/n));\n c.textBaseline=\"middle\"; \n for(let i=0; i