44:["$","title",null,{"children":"Discrete Medial Axis | Luis Retondaro | Observable"}] 4c:Tb0f,compute_RDMA = function (SDT) { let Dx = (i) => SDT.map(function(value,index) { return value[i]; }); let D_left = (d) => d.map(function(value,index) {return value==0?0:Math.ceil(index-(Math.sqrt(value)*2+Math.sqrt(value)*0.1))}); let D_right = (d) => d.map(function(value,index) {return value==0?0:Math.floor(index+(Math.sqrt(value)*2+Math.sqrt(value)*0.1))}); let D_obj = (l, r) => (r).map(function(item,i) { return {"left": l[i], "right": r[i], "index": i }; }); let ordering = function (a, b) { return a["left"] - b["left"] || a["right"] - b["right"]; } let Incl = function (ob, x1, x2) { let ret; if ((ob[x1].left == ob[x2].left) && (ob[x1].right == ob[x2].right)) { ret = true; } else { ret = (ob[x2].left >= ob[x1].left && ob[x2].right <= ob[x1].right) && (ob[x1].left != ob[x2].left || ob[x1].right != ob[x2].right); } return ret; } let _Rdma = SDT.map((d)=>d.map((e)=>0)); // A zero matrix let _rows = SDT[0].length; let _cols = SDT.length; let L = []; let _q; let _S; //======== Step 1 -> cols (x-axis) ========= for (var j=0; j < _rows; j++) { let dx = Dx(j); let ddd = dx.map(function(v,i){return -1;}); let left_x = D_left(dx); let right_x = D_right(dx); let dxobj = D_obj(left_x,right_x); dxobj.sort(ordering); _q = 0; _S = [0]; for (var k=1; k < _cols; k++) { if (!Incl(dxobj, _S[_q], k)) { if (Incl(dxobj, k, _S[_q])) { _S[_q] = -1; } _q++; _S[_q] = k; } } for (var l=0; l < _S.length; l++) { if (_S[l] != -1) { _Rdma[dxobj[_S[l]].index][j] = dx[dxobj[_S[l]].index]; ddd[dxobj[_S[l]].index] = dx[dxobj[_S[l]].index]; } } } //======== Step 2 -> rows (y-axis) ========= for (var i=0; i < _cols; i++) { let dy = _Rdma[i].map(function(v,i){return v;}); let ddd = dy.map(function(v,i){return v-1;}); for (var it = 0; it < _rows; it++) { _Rdma[i][it] = 0; } let left_y = D_left(dy); let right_y = D_right(dy); let dyobj = D_obj(left_y,right_y); dyobj.sort(ordering); _q = 0; _S = [0]; for (var k=1; k < _rows; k++) { if (!Incl(dyobj, _S[_q], k)) { if (Incl(dyobj, k, _S[_q])) { _S[_q] = -1; } _q++; _S[_q] = k; } } for (var l=0; l < _S.length; l++) { if (_S[l] != -1) { _Rdma[i][dyobj[_S[l]].index] = dy[dyobj[_S[l]].index]; ddd[dyobj[_S[l]].index] = dy[dyobj[_S[l]].index]; } } } return _Rdma; } 4d:[] 48:{"id":"569b9b65d83fb6f6","slug":"discrete-medial-axis","title":"Discrete Medial Axis","thumbnail":"0c57b4acc3d8b96ea036be14f535dee0fb958a5548a02d6ca2b7249feda4347c","access_level":"public","create_time":"$D2019-06-24T17:18:11.950Z","update_time":"$D2022-03-10T18:04:30.935Z","view_time":null,"first_public_time":"$D2019-06-24T17:30:14.507Z","has_unpublished_edits":null,"fork":false,"stars":2,"starred":false,"owner":{"type":"team","id":"71e5ea5268419c5a","login":"lretondaro","name":"Luis Retondaro","avatar_url":"https://avatars-next.observableusercontent.com/d5e7faa70402437c197456859137585a1f5ae86f61b37fee567d63511c467d09","tier":"starter_2024","bio":".","followed_by_you":false,"follows_you":false},"version":285,"publish_version":285,"role":null,"authored":false,"authors":[{"type":"individual","id":"146ee3f34bbcc61d","login":"lretondaro","name":"Luis Retondaro","avatar_url":"https://avatars-next.observableusercontent.com/d5e7faa70402437c197456859137585a1f5ae86f61b37fee567d63511c467d09","tier":null,"bio":".","followed_by_you":false,"follows_you":false}],"members":[{"type":"individual","id":"146ee3f34bbcc61d","login":"lretondaro","name":"Luis Retondaro","avatar_url":"https://avatars-next.observableusercontent.com/d5e7faa70402437c197456859137585a1f5ae86f61b37fee567d63511c467d09","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"lapsed":null,"subscribed":null,"forks":0,"collection":{"id":"e35da0991f30f391","slug":"computer-graphics","title":"Computer Graphics","owner":{"id":"71e5ea5268419c5a","login":"lretondaro","name":"Luis Retondaro"}},"type":"notebook","body":{"titling":"auto","title":"Discrete Medial Axis","cells":[{"id":0,"value":"md`# Discrete Medial Axis`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":5,"value":"md `This implementation is based on paper: \n\t[Optimal Separable Algorithms to Compute the Reverse Euclidean Distance Transformation\n\tand Discrete Medial Axis in Arbitrary Dimension](https://arxiv.org/pdf/0705.3343.pdf). At:\n\tIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 29, NO. 3, MARCH 2007\n\nAuthors:\n\tDavid Coeurjolly and Annick Montanvert `","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":15,"value":"html `

\n\n
In binary images, the distance transformation (DT) are classic tool for the geometrical skeleton and the medial axis (MA) extraction. In digital space, the MA is a convenient tool to represent shapes since it is reversible: from the MA points, we can exactly reconstruct the original shape.\n
This algorithm is based on the discrete approach proposed for the Reverse Euclidean Distance Transformation\n(REDT).
\n\n\n

`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":58,"value":"md `## 1. The separable techniques to compute the Euclidean distance transformation of *d-dimensional* images.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":65,"value":"md `In the 2D case the authors consider a two-dimensional binary\nimage ${tex `P` } of size ${tex `n \\times n`}; ${tex `\\overline P`} denotes the complementary of ${tex `P`} ,\ni.e., the set of background pixels. The output of the\nalgorithm is a 2D image ${tex `H = (i,j)`} storing the squared\ndistance transformation (SDT). For each point ${tex `(i,j)`} of the image,\nthe squared distance transformation is given by:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":80,"value":"tex.block `h(i,j) = min \\{(i - x)^2 + (j-y)^2; 0 \\leq x, y < and (x,y) \\in \\overline P) \\} `","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":30,"value":"md `### Function SDT_step1\n\nBuild, from the source image ${tex `P`}, a one-dimensional Euclidean Distance Transformation (EDT) according to the first dimension (x-axis) denoted by ${tex `G = \\{g(i,j) \\}`}, where, for a given row ${tex `j`}:\n\n${tex.block `g(i,j) = \\underset {x}{min} \\{ \\left | i-x\\right |; 0 \\leq x < n, and (x,j) \\in \\overline P\\} `} `","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":103,"value":"{ \n let div = html `

`;\n let img = html `

`\n div.append(img);\n return div;\n}\n ","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":9,"value":"SDT_step1 = function(_P) {\n\n\n \tvar _rows = _P[0].length;\n \tvar _cols = _P.length;\n\n var mmap = _P.map(function(v,i){return v;});\n\n\n let _G = mmap.map((d)=>d.map((e)=>0));\n \n for (var j=0; j < _rows; j++) {\n\n\n \tif (mmap[0][j] == 0) { //it is background\n \t\t_G[0][j] = 0;\n \t} else {\n \t\t_G[0][j] = Number.POSITIVE_INFINITY;\n \t}\n \tfor (var i=1; i < _cols; i++) {\n \t\tif (mmap[i][j] == 0) { //it is background\n \t\t\t_G[i][j] = 0;\n \t\t} else {\n _G[i][j] = 1 + _G[i-1][j];\n \t\t}\n \t}\n\n \tfor (var i=_cols-2; i >= 0; i--) {\n \t\tif (_G[i+1][j] < _G[i][j]) {\n \t\t\t_G[i][j] = 1 + _G[i+1][j];\n \t\t}\n \t}\n }\n\n return _G;\n}\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":183,"value":"md `### Function SDT_step2\n\nThen, construct the image ${tex `H = (i,j)`} with a y-axis process:\n\n${tex.block `h(i,j) = \\underset {y}{min} \\{ g(i,y)^2 + (j-y)^2; 0 \\leq y < n \\} `} `","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":12,"value":"SDT_step2 = function(_G) {\n\n\tlet sep = (u, v, i) => ~~(((v**2) - (u**2) + (_G[i][v]**2) - (_G[i][u]**2) ) / (2*(v - u)));\n let F = (i,j,y) => _G[i][y]*_G[i][y]+(j-y)*(j-y);\n\n\tlet _H = _G.map((d)=>d.map((e)=>0)); // A zero matrix\n let _rows = _G[0].length; \n let _cols = _G.length;\n\n\n for (var i=0; i < _cols; i++) {\n \tlet _q = 0, _S = [0], _T = [0];\n\n \tfor (var j=1; j < _rows; j++) {\n \t\twhile (_q >= 0 && F(i,_T[_q],_S[_q],) > F(i, _T[_q], j)) _q--;\n\n \t\tif (_q < 0) {\n \t\t\t_q = 0;\n \t\t\t_S = [j];\n \t\t} else {\n \t\t\tlet _w = 1 + sep(_S[_q], j, i);\n \t\t\t//console.assert (!Number.isNaN(w));\n \t\t\tif (_w < _rows) {\n \t\t\t\t_q++;\n \t\t\t\t_S[_q] = j; \n \t\t\t\t_T[_q] = _w;\n \t\t\t}\n\n \t\t}\n \t}\n \tfor (var j=_rows-1; j >= 0; j--) {\n \t\t_H[i][j] = F(i, j, _S[_q]);\n\n \t\tif (j == _T[_q]) {\n \t\t\t_q--;\n \t\t}\n \t}\n\n }\n\n return _H;\n\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":18,"value":"REDT_step1 = function (_F) {\n\n\tlet sep = (u, v, j) => ~~(((u*u) - (v*v) - _F[u][j] + _F[v][j] ) / (2*(u - v)));\n let F = (j,i,x) => _F[x][j]-((i-x)*(i-x));\n\n\tlet _G = _F.map((d)=>d.map((e)=>0)); // A zero matrix\n let _rows = _F[0].length; \n let _cols = _F.length;\n\n for (var j=0; j < _rows; j++) {\n \tlet _q = 0, _S = [0], _T = [0];\n\n \tfor (var i=1; i < _cols; i++) {\n\n \t\twhile (_q >= 0 && F(j,_T[_q],_S[_q]) < F(j, _T[_q], i)) _q--;\n\n\n \t\tif (_q < 0) {\n \t\t\t_q = 0;\n \t\t\t_S = [i];\n \t\t} else {\n \t\t\tlet _w = 1 + sep(_S[_q], i, j);\n \t\t\t//console.assert (!Number.isNaN(w));\n \t\t\tif (_w < _cols) {\n \t\t\t\t_q++;\n \t\t\t\t_S[_q] = i; \n \t\t\t\t_T[_q] = _w;\n \t\t\t}\n\n \t\t}\n \t}\n \tfor (var i=_cols-1; i >= 0; i--) {\n \t\t_G[i][j] = F(j, i, _S[_q]);\n\n \t\tif (i == _T[_q]) {\n \t\t\t_q--;\n \t\t}\n \t}\n\n }\n\n return _G;\n\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":193,"value":"md `## 2. Optimal SDT`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":198,"value":"md `For the second step, the ${tex `min`} operation corresponds to a lower envelope computation of a set of parabolas. More precisely, let us suppose that we have computed *Step 1* of the algorithm (x-axis SDT) and let ${tex `\\{g(i,y)\\} (0 \\leq y < n)`} be a column of G. If we consider the set of parabolas ${tex `F_y^i(j) = g(i,y)^2 + (j-y)^2`}, the column ${tex `h(i,j)`} after Step 2 is exactly the lower envelope of ${tex `\\{F_y^i(j)\\}`} with ${tex `0 \\leq y < n`}.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":211,"value":"{ \n let div = html `

`;\n let img = html `