*! delaunay v1.2: A Stata package for Delaunay triangulation, Convex hull, and Voronoi tesselations *! Asjad Naqvi (asjadnaqvi@gmail.com, @AsjadNaqvi) * * Ver 1.11 05.03.2022: x and y flipped. faster export to Stata. warnings added. * Ver 1.10 01.03.2022: missing triangles fixed. rays fixed. * Ver 1.02 27.12.2021: rescale added. if made optional * Ver 1.01 20.12.2021: if/in & notes by wbuchanan * Ver 1.00 22.11.2021: first release ************************** ************************** ** ** ** Delaunay ** ** ** ** Asjad Naqvi ** ** ** ** Last updated: ** ** 04 Sep 2022 ** ** ** ** First release: ** ** 22 Nov 2021 ** ** ** ************************** ************************** cap prog drop delaunay prog def delaunay, sortpreserve version 15 syntax varlist(min = 2 max = 2 numeric) [if] [in], /// [ REScale ] /// // normalize x and y to [0,1] calculate stuff and rescale back [ TRIangles addbox ] /// // export triangles with or without a bounding box [ hull replace ] /// // export the hull [ VORonoi(namelist min=1 max=2) ] /// // export the Voronoi polygons [ OFFset(real 0.05) ] ************************ ***** Error checks ***** ************************ marksample touse, strok gettoken y x : varlist local obs1 = _N // check voronoi locals local vorcount = wordcount("`voronoi'") if (`vorcount' == 1) { if ("`voronoi'" != "lines" & "`voronoi'" != "polygons") { di as err "Valid options for voronoi() are {ul:lines} or {ul:polygons} or both." exit } } if (`vorcount' == 2) { tokenize `voronoi' if ("`1'" != "lines" & "`1'" != "polygons") & ("`2'" != "lines" & "`2'" != "polygons") & ("`1'" != "`2'") { di as err "Valid options for voronoi() are {ul:lines} or {ul:polygons} or both." exit } } if "`replace'" == "" { if "`triangles'" != "" { cap confirm file _triangles.dta if _rc == 0 { di as error "File {it:_triangles.dta} already exists. Please use the {ul:replace} option if you want to overwrite it." exit } } if "`voronoi'" != "" { cap confirm file _vorlines.dta if _rc == 0 { di as error "File {it:_vorlines.dta} already exists. Please use the {ul:replace} option if you want to overwrite it." exit } cap confirm file _vorpoly.dta if _rc == 0 { di as error "File {it:_vorpoly.dta} already exists. Please use the {ul:replace} option if you want to overwrite it." exit } } cap confirm variable _ID if _rc == 0 { di as error "Variable {it:_ID} already exists. Please use the {ul:replace} option if you want to overwrite it." exit } } // advice rescale if "`rescale'" == "" { qui su `y' if `touse', meanonly local ymean = r(mean) qui su `x' if `touse', meanonly local xmean = r(mean) if (`ymean' > 2 * `xmean') | (`xmean' > 2 * `ymean') { di in yellow "Warning: X and Y are likely not on the same scale. If X and Y are not geographic coordinates," di in yellow "then the use of the {ul:rescale} option is highly recommended. See help file." di _newline } } ************************** **** End error checks **** ************************** di _newline di in yellow "Initializing >" // generate the _id variable cap confirm variable _ID if !_rc { drop if _ID==. drop _ID } gen _ID = _n if `touse' lab var _ID "Observation ID" mata: points = select(st_data(., ("`x'", "`y'")), st_data(., "`touse'")) if "`addbox'" != "" { mata: trixmin = .; trixmax = .; triymin = .; triymax = . mata: addbounds(points, trixmin, trixmax, triymin, triymax, `offset') // now we create another copy with bounds mata: points = points \ trixmin, triymin \ trixmin, triymax \ trixmax, triymax \ trixmax, triymin } if "`rescale'" != "" { mata: points2 = points mata: points2[.,1] = rescale(points[.,1], 0, 1) mata: points2[.,2] = rescale(points[.,2], 0, 1) } else { mata: points2 = points } mata: eps = 1e-20; edgestack = J(512, 1, .) mata: coords = initialize(points2); num = rows(points2) mata: maxTriangles = max(((2 * num) - 5, 0)) // core arrays mata: triangles = J(maxTriangles * 3, 1, .) mata: halfedges = J(maxTriangles * 3, 1, .) mata: hull = J(num, 1, .) // arrays for tracking mata: hashSize = ceil(sqrt(num)) mata: hullPrev = J(num, 1, .); hullNext = J(num, 1, .); hullTri = J(num, 1, .); hullHash = J(hashSize, 1, -1) mata: ids = J(num, 1, .); dists = J(num, 1, .) // temporary arrays for sorting points mata: ids = J(num, 1, .); dists = J(num, 1, .) di in yellow "Starting core routines >" // run the core routine mata: _delaunay_core(coords, ids, dists, triangles, halfedges, hull, hullNext, /// hullPrev, hullTri, hullHash, hashSize, eps, edgestack) ************************** *** export geometry *** ************************** // delaunay triangles if "`triangles'" != "" { di in yellow "Exporting geometry > Triangles exported as {it:_triangles.dta}." add_triangles } // convex hull if "`hull'" != "" { di in yellow "Exporting geometry > Hull added as variables." add_hull } // voronoi lines and polygons if "`voronoi'" != "" { di in yellow "Exporting geometry > Voronoi exported as {it:_vorlines.dta}/{it:_vorpoly.dta}." if "`offset'" != "" { voronoi `varlist', `voronoi' offset(`offset') } else { voronoi `varlist', `voronoi' } } // clear mata mata: mata drop coords dists edgestack eps halfedges hashSize hull hullHash hullNext hullPrev hullTri ids maxTriangles num points points2 triangles // core cap mata: mata drop mytriangles triindex // triangles cap mata: mata drop myhull // hull cap mata: mata drop vorlines vorpolygons // voronoi end ************************ // _delaunay_core // main routine ************************ cap mata: mata drop _delaunay_core() mata: // _delaunay_core void _delaunay_core(coords,ids,dists,triangles,halfedges,hull,hullNext,hullPrev,hullTri,hullHash,hashSize,eps,edgestack) { num = rows(coords) / 2 minX = 1e16 minY = 1e16 maxX = -1e16 // arbitrary large numbers maxY = -1e16 for (i=1; i<= num; i++) { x = coords[2 * i - 1, 1] y = coords[2 * i , 1] if (x < minX) minX = x if (y < minY) minY = y if (x > maxX) maxX = x if (y > maxY) maxY = y ids[i,1] = i } cx = (minX + maxX) / 2 cy = (minY + maxY) / 2 minDist = . i0 = 0 i1 = 0 i2 = 0 // pick a seed point close to the center (first observation) for (i=1; i<=num; i++) { dval = dist(cx, cy, coords[2 * i - 1], coords[2 * i]) if (dval < minDist) { i0 = i minDist = dval } } i0x = coords[2 * i0 - 1, 1] i0y = coords[2 * i0 , 1] minDist = . // find the point closest to the seed (second observation) for (i=1; i<=num; i++) { if (i == i0) continue dval = dist(i0x, i0y, coords[2 * i - 1], coords[2 * i]) if ((dval < minDist) & (dval > 0)) { i1 = i minDist = dval } } i1x = coords[2 * i1 - 1, 1] i1y = coords[2 * i1 , 1] minRadius = . // find the third point which forms the smallest circumcircle with the first two (third observation) for (i=1; i<=num; i++) { if (i == i0 | i == i1) continue r = circumradius(i0x, i0y, i1x, i1y, coords[2 * i - 1, 1], coords[2 * i, 1]) if (r < minRadius) { i2 = i minRadius = r } } i2x = coords[2 * i2 - 1, 1] i2y = coords[2 * i2 , 1] if (minRadius == .) { // calculate the distances ot the points with the first point for (i=1; i<=num; i++) { dists[i,1] = (coords[2 * i - 1, 1] - coords[1, 1]) | (coords[2 * i, 1] - coords[2, 1]) } // and order them _quicksort(ids, dists, 1, num - 1) hull = J(num, 1, .) j = 1 d0 = -1e20 // minus infinity for (i=1; i <= num; i++) { id = ids[i,1] if (dists[id,1] > d0) { hull[j, 1] = id j = j + 1 d0 = dists[id,1] } } // initiate the hull, triangle, and halfedges matrices hull = hull[1..j,1] triangles = . halfedges = . } // swap the order of the seed points for counter-clockwise orientation if (orient(i0x, i0y, i1x, i1y, i2x, i2y)) { i = i1 x = i1x y = i1y i1 = i2 i1x = i2x i1y = i2y i2 = i i2x = x i2y = y } c0x = . // blanks to be evaluated c0y = . circumcenter(i0x, i0y, i1x, i1y, i2x, i2y, c0x, c0y) for (i=1; i <= num; i++) { // changed from <= to < dists[i,1] = dist(coords[2 * i - 1,1], coords[2 * i,1], c0x, c0y) } // sort the points by distance from the seed triangle circumcenter _quicksort(ids, dists, 1, num - 1) // set up the seed triangle as the starting hull hullStart = i0 hullSize = 3 hullNext[i0] = i1 hullPrev[i2] = i1 hullNext[i1] = i2 hullPrev[i0] = i2 hullNext[i2] = i0 hullPrev[i1] = i0 // initialize the hull triangles hullTri[i0] = 1 hullTri[i1] = 2 hullTri[i2] = 3 // manual index fix to match the py output hullHash[hashKey(i0x,i0y,cx,cy,hashSize) + 1] = i0 hullHash[hashKey(i1x,i1y,cx,cy,hashSize) + 1] = i1 hullHash[hashKey(i2x,i2y,cx,cy,hashSize) + 1] = i2 trianglesLen = 1 // dump is just a dummy to prevent the function from returning a value on the screen dump = addTriangle(triangles,trianglesLen,halfedges, i0, i1, i2, -1, -1, -1) xp=. yp=. /////////////////////////////////// /// main loop below /// /////////////////////////////////// for (k=1; k <= rows(ids); k++) { // evalaute the sorted points i = ids[k,1] x = coords[2 * i - 1, 1] y = coords[2 * i , 1] // skip near-duplicate points if ((k > 1) & (abs(x - xp) <= eps) & (abs(y - yp) <= eps)) continue xp = x yp = y // skip seed triangle points if ((i == i0) | (i == i1) | (i == i2)) continue // find a visible edge on the convex hull using edge hash start = 1 key = hashKey(x, y, c0x, c0y, hashSize) for (j=1; j <= hashSize; j++) { start = hullHash[(mod((key + j - 1), hashSize) + 1),1] if ((start != -1) & (start != hullNext[max((start,1)), 1])) break // added max to prevent the 0 argument } start = hullPrev[start,1] emp = start www = 0 // while while (www == 0) { // here is where points are being dropped q = hullNext[emp,1] if (orient(x, y, coords[2 * emp - 1], coords[2 * emp], coords[2 * q - 1], coords[2 * q])) break emp = q if (emp == start) { emp = -1 break } } if (emp == -1) continue // add the first triangle from the point t = addTriangle(triangles, trianglesLen, halfedges, emp, i, hullNext[emp,1], -1, -1, hullTri[emp,1]) // recursively flip triangles from the point until they satisfy the Delaunay condition hullTri[i , 1] = legalize(t + 2, coords, halfedges, edgestack, triangles, hullStart, hullTri, hullPrev) hullTri[emp, 1] = t hullSize = hullSize + 1 // walk forward through the hull, adding more triangles and flipping recursively num = hullNext[emp, 1] www = 0 // the infinite while while (www == 0) { q = hullNext[num, 1] if (orient(x, y, coords[2 * num - 1,1], coords[2 * num,1], coords[2 * q - 1,1], coords[2 * q,1]) == 0) break t = addTriangle(triangles,trianglesLen, halfedges,num, i, q, hullTri[i,1], -1, hullTri[num,1]) hullTri[i,1] = legalize(t + 2, coords,halfedges,edgestack,triangles,hullStart,hullTri,hullPrev) hullNext[num,1] = num // mark as removed hullSize = hullSize - 1 num = q } // walk backward from the other side, adding more triangles and flipping if (emp == start) { www = 0 while (www == 0) { q = hullPrev[emp,1] if ((orient(x, y, coords[2 * q - 1], coords[2 * q], coords[2 * emp - 1], coords[2 * emp])) == 0) break t = addTriangle(triangles, trianglesLen, halfedges, q, i, emp, -1, hullTri[emp,1], hullTri[q,1]) dump = legalize(t + 2, coords, halfedges, edgestack, triangles, hullStart, hullTri, hullPrev) hullTri[q,1] = t hullNext[emp,1] = emp // mark as removed hullSize = hullSize - 1 emp = q } } // update the hull indices hullStart = emp hullPrev[i,1] = emp hullNext[emp,1] = i hullPrev[num,1] = i hullNext[i,1] = num // save the two new edges in the hash table hullHash[hashKey(x,y,cx,cy,hashSize) + 1,1] = i hullHash[hashKey(coords[2 * emp - 1,1], coords[2 * emp,1],cx,cy,hashSize) + 1,1] = emp } hull = J(hullSize,1,.) emp = hullStart for (i=1; i<=hullSize; i++) { hull[i,1] = emp emp = hullNext[emp,1] } } end **************************************************** **************** SUBROUTINES ******************** **************************************************** ********************* // initialize // ********************* cap mata: mata drop initialize() mata: // initialize real vector initialize(real matrix data) { real scalar num, i real vector coords num = rows(data) coords = J(num*2,1,.) for (i=1; i<=num; i++) { coords[2 * i - 1, 1] = data[i,1] coords[2 * i , 1] = data[i,2] } return (coords) } end ******************** // _hashKey // ******************** cap mata: mata drop hashKey() mata // _hashKey real scalar hashKey(x,y,cx,cy,hashSize) { real scalar tt tt = mod(floor(pseudoAngle(x - cx, y - cy) * hashSize), hashSize) // new code to fix return . (4 Sep 2022) if (tt != .) { return (tt) } else { return (0) } } end ************************ // pseudoAngle // ************************ cap mata: mata drop pseudoAngle() mata: // pseudoAngle function pseudoAngle(real scalar dx, real scalar dy) { real scalar p p = dx / (abs(dx) + abs(dy)) if (dy > 0) { return ((3 - p) / 4) } else { return ((1 + p) / 4) } } end ************************ // _legalize // ************************ cap mata: mata drop legalize() mata // _legalize function legalize(real scalar a, real vector coords, halfedges, edgestack, triangles, hullStart, hullTri, hullPrev) { i = 1 ar = 0 www = 0 while (www==0) { b = halfedges[a,1] /* if the pair of triangles doesn't satisfy the Delaunay condition (p1 is inside the circumcircle of [p0, pl, pr]), flip them, then do the same check/flip recursively for the new pair of triangles pl pl /||\ / \ al/ || \bl al/ \a / || \ / \ / a||b \ flip /___ar___\ p0\ || /p1 => p0\---bl---/p1 \ || / \ / ar\ || /br b\ /br \||/ \ / pr pr */ a0 = a - mod(a - 1, 3) ar = a0 + mod(a + 1, 3) if (b == -1) { // convex hull edge if (i == 1) break i = i - 1 a = edgestack[i,1] // define in the other file. empty array continue } b0 = b - mod(b - 1,3) al = a0 + mod(a, 3) bl = b0 + mod((b + 1),3) p0 = triangles[ar,1] pr = triangles[a ,1] pl = triangles[al,1] p1 = triangles[bl,1] illegal = inCircle(coords[2 * p0 - 1, 1], coords[2 * p0, 1], coords[2 * pr - 1], coords[2 * pr, 1], coords[2 * pl - 1, 1], coords[2 * pl, 1], coords[2 * p1 - 1, 1], coords[2 * p1, 1]) if (illegal) { triangles[a,1] = p1 triangles[b,1] = p0 hbl = halfedges[bl, 1] // edge swapped on the other side of the hull (rare); fix the halfedge reference if (hbl == -1) { eme = hullStart yyy = 0 while (yyy==0) { if (hullTri[eme, 1] == bl) { hullTri[eme, 1] = a break } eme = hullPrev[eme, 1] if (eme == hullStart) break } } _link(halfedges,a, hbl) _link(halfedges,b, halfedges[ar, 1]) _link(halfedges,ar, bl) brme = b0 + mod(b,3) if (i < rows(edgestack)) { edgestack[i,1] = brme i = i + 1 } } else { if (i == 1) break i = i - 1 a = edgestack[i,1] } } return (ar) } end ******************** // _link // ******************** cap mata: mata drop _link() mata: // _link void _link(real vector halfedges, real scalar a, real scalar b) { halfedges[a,1] = b if (b != -1) { halfedges[b,1] = a } } end ************************ // _addTriangle // ************************ cap mata: mata drop addTriangle() mata: // addTriangle function addTriangle(triangles, trianglesLen, halfedges, i0, i1, i2, a, b, c) { real scalar t t = trianglesLen triangles[t ] = i0 triangles[t + 1] = i1 triangles[t + 2] = i2 _link(halfedges, t , a) _link(halfedges, t + 1, b) _link(halfedges, t + 2, c) trianglesLen = trianglesLen + 3 return (t) } end **************** // dist // **************** cap mata: mata drop dist() mata: // dist real vector dist(ax, ay, bx, by) { real scalar dx, dy dx = ax - bx dy = ay - by return (dx*dx + dy*dy) } end ************************ // orientIfSure // ************************ cap mata: mata drop orientIfSure() mata: // orientIfSure function orientIfSure(px, py, rx, ry, qx, qy) { real scalar left, right left = (ry - py) * (qx - px) right = (rx - px) * (qy - py) // J. Shewchuk error bound check if (abs(left - right) >= (3.3306690738754716e-16 * abs(left + right))) { return(left - right) } else { return (0) } } end ******************** // orient // ******************** cap mata: mata drop orient() mata: // orient function orient(px, py, rx, ry, qx, qy) { real scalar x x = (orientIfSure(px, py, rx, ry, qx, qy) < 0 | orientIfSure(rx, ry, qx, qy, px, py) < 0 | orientIfSure(qx, qy, px, py, rx, ry) < 0 ) return (x) // a true false statement } end ******************** // inCircle // ******************** cap mata: mata drop inCircle() mata: // inCircle function inCircle(ax, ay, bx, by, cx, cy, px, py) { real scalar dx, dy, ex, ey, fx, fy, ap, bp, cp, x dx = ax - px dy = ay - py ex = bx - px ey = by - py fx = cx - px fy = cy - py ap = (dx * dx) + (dy * dy) bp = (ex * ex) + (ey * ey) cp = (fx * fx) + (fy * fy) x = (dx * (ey * cp - bp * fy)) - (dy * (ex * cp - bp * fx)) + (ap * (ex * fy - ey * fx)) return (x < 0) } end ************************ // circumradius // ************************ cap mata: mata drop circumradius() mata: // circumradius function circumradius(ax, ay, bx, by, cx, cy) { real scalar dx, dy, ex, ey, bl, cl, d, x, y dx = bx - ax dy = by - ay ex = cx - ax ey = cy - ay bl = (dx * dx) + (dy * dy) cl = (ex * ex) + (ey * ey) d = 0.5 / ((dx * ey) - (dy * ex)) x = ((ey * bl) - (dy * cl)) * d y = ((dx * cl) - (ex * bl)) * d return (x*x + y*y) } end ************************ // circumcenter // ************************ cap mata: mata drop circumcenter() mata: // circumcenter function circumcenter(ax, ay, bx, by, cx, cy, real scalar retx, real scalar rety) { real scalar dx, dy, ex, ey, bl, cl, d dx = bx - ax dy = by - ay ex = cx - ax ey = cy - ay bl = (dx * dx) + (dy * dy) cl = (ex * ex) + (ey * ey) d = 0.5 / ((dx * ey) - (dy * ex)) retx = ax + ((ey * bl) - (dy * cl)) * d // return x rety = ay + ((dx * cl) - (ex * bl)) * d // return y return } end ************************ // quicksort // ************************ cap mata: mata drop _quicksort() // function updated mata: // quicksort void _quicksort(ids, dists, left, right) { real scalar i, temp, tempDist, j, medianme, xxx, yyy, zzz // if ((right - left) <= 20) { // taken out completely. was causing points to be dropped or not sorted correctly for (i = left + 1; i <= right + 1; i++) { temp = ids[i,1] tempDist = dists[temp,1] j = i - 1 while ((j >= left) & (dists[ids[max((j,1)),1],1] > tempDist)) { // if j hits zero, the second condition cannot be checked. so adding the max condition here ids[j + 1,1] = ids[j,1] j = j - 1 } ids[j + 1,1] = temp } } end ************************ // swapme // ************************ cap mata: mata drop swapme() mata: // swapme function swapme(arr,i,j) { real scalar tmp1 tmp1 = arr[i,1] arr[i,1] = arr[j,1] arr[j,1] = tmp1 } end ******************** // rescale // rescale a vector column ******************** cap mata: mata drop rescale() mata: // rescale real vector rescale(points, a, b) { newpoints = (b - a) * (points :- colmin(points)) :/ (colmax(points) - colmin(points)) :+ a return (newpoints) } end ******************************** *** END OF SUBROUTINES *** ******************************** ***************************************** *** routines for adding data back *** ***************************************** ************************ // return index // for exporting and merging with original data ************************ cap mata: mata drop returnindex() mata: // returnindex function returnindex(triangles4) { real vector myindex myindex = J(rows(triangles4),1,.) for (i = 1; i <= rows(triangles4); i++ ) { myindex[i,1] = floor((i - 1) / 4) + 1 } return (myindex) } end ******************** // addbounds // ******************** cap mata: mata drop addbounds() mata // addbounds function addbounds(points, xmin, xmax, ymin, ymax, offs) { displacex = abs((max(points[.,1]) - min(points[.,1])) * offs) displacey = abs((max(points[.,2]) - min(points[.,2])) * offs) xmin = min(points[.,1]) - displacex xmax = max(points[.,1]) + displacex ymin = min(points[.,2]) - displacey ymax = max(points[.,2]) + displacey } end ************************ // triangles // ************************ cap mata: mata drop fixtriangles() mata: // fix triangles function fixtriangles(inputtri,points) { // duplicate the first observation exporttri = inputtri, J(rows(inputtri),2,.) for (i = 1; i <= rows(exporttri); i++ ) { if (exporttri[i,2] != .) { exporttri[i,3..4] = points[exporttri[i,2],.] } } return (exporttri) } end cap program drop add_triangles program define add_triangles, sortpreserve version 15 // these intermediate names will eventually be removed mata: triangles = select(triangles, (triangles[.,1] :< .)) // drop the missing rows mata: triangles2 = colshape(triangles',3) // // convert triplets to quadruplets. need to add a black row for drawing mata: triangles3 = triangles2, J(rows(triangles2),1,.) mata: triangles4 = colshape(triangles3,1) mata: triindex = returnindex(triangles4) mata: triangles5 = triindex,triangles4 mata: mytriangles = fixtriangles(triangles5,points) qui { preserve clear getmata (_group _ID _X _Y) = mytriangles, force double replace gen serial = _n order serial lab var serial "Triangle: Serial" lab var _group "Triangle: group" lab var _ID "Triangle: Observation _ID" lab var _X "Triangle: X" lab var _Y "Triangle: Y" compress save _triangles.dta, replace restore } *mata: st_matrix("triangles", mytriangles) mata: mata drop triangles2 triangles3 triangles4 triangles5 cap drop tri_num // dud end ******************** // hull // ******************** cap mata: mata drop fixhull() mata: // generate hull function fixhull(hull,points) { // duplicate the first observation hull2 = hull \ . hull2[rows(hull2),1] = hull[1,1] hull3 = J(rows(hull2),1,.), hull2, J(rows(hull2),2,.) for (i = 1; i <= rows(hull2); i++ ) { hull3[i,1] = i hull3[i,3..4] = points[hull3[i,2],.] } return(hull3) } end ********************** // export the hull // ********************** cap program drop add_hull program define add_hull version 15 cap drop hull* // make sure the variables are clear mata: myhull = fixhull(hull,points) mata st_matrix("hull", myhull) mat colnames hull = "hull_num" "hull_ID" "hull_X" "hull_Y" svmat hull, n(col) cap drop hull_num lab var hull_ID "Hull: Observation _ID" lab var hull_X "Hull: X" lab var hull_Y "Hull: Y" cap mat drop hull end