31:T42c,advect1 = { let p, q, orderp, orderq; const random = d3.randomNormal(); const { hypot } = Math; return function advect1(values, target, epsilon) { const l = values.length / 4; if (!(p && p.length === l)) { p = new Float32Array(l); q = new Float32Array(l); orderp = new Uint32Array(l); orderq = new Uint32Array(l); for (let i = 0; i < l; i++) orderp[i] = orderq[i] = i; } let x = random(), y = random(), z = random(); const norm = 1 / hypot(x, y, z); x *= norm; y *= norm; z *= norm; p = project_pixels(values, x, y, z, p); q = project_pixels(target, x, y, z, q); orderp.sort((i, j) => p[i] - p[j]); orderq.sort((i, j) => q[i] - q[j]); let energy = 0; for (let k = 0; k < orderp.length; k++) { const i = orderp[k], j = orderq[k]; const a = q[j] - p[i], b = a * epsilon; values[4 * i + 0] += b * x; values[4 * i + 1] += b * y; values[4 * i + 2] += b * z; energy += Math.abs(a); } return energy; }; }32:T542,// Compute a batch of transformations and apply them at once advectBatch = { let orderp, orderq, p, q, batch; const random = d3.randomNormal(); const { hypot } = Math; return function advectBatch(values, target, n, epsilon) { const l = values.length / 4; if (!(p && p.length === l)) { p = new Float32Array(l); q = new Float32Array(l); batch = new Float32Array(values.length); orderp = new Uint32Array(l); orderq = new Uint32Array(l); for (let i = 0; i < l; i++) orderp[i] = orderq[i] = i; } batch.fill(0); for (let s = 1; s < n; s++) { let x = random(), y = random(), z = random(), norm = 1 / hypot(x, y, z); x *= norm; y *= norm; z *= norm; p = project_pixels(values, x, y, z, p); q = project_pixels(target, x, y, z, q); orderp.sort((i, j) => p[i] - p[j]); orderq.sort((i, j) => q[i] - q[j]); for (let k = 0; k < l; k++) { const i = orderp[k], j = orderq[k]; const b = q[j] - p[i]; batch[4 * i + 0] += b * x; batch[4 * i + 1] += b * y; batch[4 * i + 2] += b * z; } } let energy = 0; for (let i = 0; i < batch.length; i++) { values[i] += batch[i] * epsilon; energy += Math.abs(batch[i]); } return energy; }; }33:[] 2f:{"title":"Sliced Optimal Transport","cells":[{"id":0,"value":"md`# Sliced Optimal Transport\n\n> _Color transfer is an image processing application where you want to retarget the color histogram of an input image according to the color histogram of a target one._\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":265,"value":"{\n const context = DOM.context2d(width, h);\n\n const ctx = DOM.context2d(w, h, 1);\n ctx.putImageData(source, 0, 0);\n context.drawImage(ctx.canvas, 0, 0);\n ctx.putImageData(target, 0, 0);\n context.drawImage(ctx.canvas, w * 2, 0);\n\n context.drawImage(transfered, w, 0);\n context.fillStyle = \"white\";\n context.fillText(`source image`, 10, 10);\n context.fillText(\n `${transfered.step} steps, batches of ${\n transfered.batchSize\n }, energy per pixel: ${transfered.energy.toFixed(2)}`,\n w + 10,\n 10\n );\n context.fillText(`target colors`, w * 2 + 10, 10);\n\n return context.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":345,"value":"{\n // histograms\n const context = DOM.context2d(width, h);\n context.fillStyle = \"white\";\n context.fillRect(0, 0, width, h);\n\n const color = [\"red\", \"green\", \"blue\"];\n\n histograms(source.data);\n context.translate(w, 0);\n histograms(transfered.getContext(\"2d\").getImageData(0, 0, w, h).data);\n\n context.translate(w, 0);\n histograms(target.data);\n\n return context.canvas;\n\n function histograms(data) {\n for (let i = 0; i < 3; i++) {\n context.strokeStyle = color[i];\n context.beginPath();\n const n = histogram(data, i);\n for (let i = 0; i < 256; i++) {\n context.lineTo(i, h - n[i] / 5);\n }\n context.stroke();\n }\n }\n function histogram(data, j) {\n const h = new Uint32Array(256);\n for (const c of data.filter((_, i) => i % 4 === j)) h[c]++;\n return h;\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":156,"value":"viewof swapTarget = checkbox([\"swap source and target\"])","pinned":false,"hidden":false,"mode":"ojs"},{"id":434,"value":"## The exercise\n\nI was lucky to (virtually) attend [the transpopt-igrv workshop](https://transpopt-igrv.sciencesconf.org/) on Optimal Transport in Nov. 2020. The introductory lesson was given by [Julie Delon](https://delon.wp.imt.fr/) (Université de Paris). The exercise in this notebook is described in “[Color Transfer via Discrete Optimal Transport using the sliced approach](https://codimd.math.cnrs.fr/s/2eRBqV9zl)”, with many references. I decided to write in JavaScript rather than C++ or python, so I guess I’ll get a bad grade.\n\nPS. Thanks to [Élie Michel](https://exppad.com/) and [Job van der Zwan](/@jobleonard) for suggestions.\n
PPS. (déc. 2023): See [Ricky Reusser’s implementation](https://observablehq.com/@rreusser/sliced-optimal-transport)","pinned":false,"hidden":false,"mode":"md"},{"id":56,"value":"md`## Project pixels on a random direction T`","pinned":false,"hidden":false,"mode":"ojs"},{"id":82,"value":"function project_pixels(pixels, x, y, z, P) {\n for (let i = 0; i < P.length; ++i)\n P[i] = pixels[4 * i] * x + pixels[4 * i + 1] * y + pixels[4 * i + 2] * z;\n return P;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":86,"value":"md`## Project and advect once\n\n(advect is a function that updates the pixels in-place towards the target, along a direction T)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":89,"value":"$31","pinned":true,"hidden":false,"mode":"ojs"},{"id":305,"value":"md`## naïve computation\n\nat each iteration, we take a random direction then advect 30% in that direction.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":142,"value":"naive = {\n await visibility();\n // image context\n const context = DOM.context2d(w, h, 1),\n imagedata = context.getImageData(0, 0, w, h),\n { data } = imagedata;\n\n const [from, to] = [source.data, target.data];\n\n // result context\n const ctx = DOM.context2d(width, h);\n\n yield ctx.canvas;\n // this both converts to float & slices the array\n const values = Float32Array.from(from);\n for (let i = 0; i < 100; i++) {\n const energy = advect1(values, to, .3);\n for (let i = 0; i < data.length; i++) data[i] = values[i];\n context.putImageData(imagedata, 0, 0);\n ctx.drawImage(context.canvas, 0, 0);\n ctx.fillStyle = \"white\";\n ctx.fillText(`${i + 1} steps, ${(energy / w / h).toFixed(2)}`, 10, 10);\n if ((i + 1) % 5 === 0) yield ctx.canvas;\n }\n yield ctx.canvas;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":308,"value":"md`## batch transformation\n\nInstead of following one direction, we take 8 directions and advect in the mean direction.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":193,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":105,"value":"md`## Advect a few times & recreate an image`","pinned":false,"hidden":false,"mode":"ojs"},{"id":24,"value":"viewof ratio = select({\n options: [1, 2, 4, 8],\n value: 4,\n description: \"reduce image by…\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":197,"value":"transfered = {\n const batchSize = 8;\n\n // image context\n const context = DOM.context2d(w, h, 1),\n c = context.canvas,\n imagedata = context.getImageData(0, 0, w, h),\n { data } = imagedata;\n\n const [from, to] = [source.data, target.data];\n\n // this both converts to float & slices the array\n const values = Float32Array.from(from);\n c.step = 0;\n c.batchSize = batchSize;\n do {\n c.energy = advectBatch(values, to, batchSize, .3 / batchSize) / w / h;\n for (let i = 0; i < data.length; i++) data[i] = Math.round(values[i]);\n context.putImageData(imagedata, 0, 0);\n c.step++;\n yield c;\n } while (!(c.energy < batchSize / 2));\n yield c;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":48,"value":"md`## Load image data\n\nreturns *source* and *target*`","pinned":false,"hidden":false,"mode":"ojs"},{"id":15,"value":"W0 = 1280","pinned":true,"hidden":false,"mode":"ojs"},{"id":16,"value":"W1 = 1024","pinned":true,"hidden":false,"mode":"ojs"},{"id":20,"value":"w = W0 / ratio","pinned":true,"hidden":false,"mode":"ojs"},{"id":22,"value":"h = W1 / ratio","pinned":true,"hidden":false,"mode":"ojs"},{"id":4,"value":"// returns the pixels as r,g,b,a,r,g,b,a,r,g,b,a…\npixelData = src =>\n d3.image(src, { crossOrigin: \"anonymous\" }).then(d => {\n const context = DOM.context2d(w, h, 1),\n canvas = context.canvas;\n context.drawImage(d, 0, 0, w, h);\n canvas.value = context.getImageData(0, 0, w, h);\n return canvas;\n })","pinned":true,"hidden":false,"mode":"ojs"},{"id":400,"value":"source = swapTarget ? image2 : image1","pinned":true,"hidden":false,"mode":"ojs"},{"id":401,"value":"target = swapTarget ? image1 : image2","pinned":true,"hidden":false,"mode":"ojs"},{"id":449,"value":"url1 = \"https://rawcdn.githack.com/dcoeurjo/transport-optimal-gdrigrv-11-2020/5507227c7bc569f524069bb2980fcfc72139c1ab/pexelA-0.png\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":33,"value":"viewof image1 = pixelData(url1)","pinned":true,"hidden":false,"mode":"ojs"},{"id":445,"value":"url2 = \"https://rawcdn.githack.com/dcoeurjo/transport-optimal-gdrigrv-11-2020/5507227c7bc569f524069bb2980fcfc72139c1ab/pexelB-0.png\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":35,"value":"viewof image2 = pixelData(url2)","pinned":true,"hidden":false,"mode":"ojs"},{"id":6,"value":"d3 = require(\"d3@6\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":158,"value":"import { checkbox, select } from \"@jashkenas/inputs\"","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}