carved = { const image = await FileAttachment("The_Scream_by_Edvard_Munch,_1893_-_Nasjonalgalleriet small@1.jpg").image(); const w = Math.min(width, image.naturalWidth); const h = Math.round(w / image.naturalWidth * image.naturalHeight); const context = DOM.context2d(w, h, 1); context.canvas.style = "max-width: 100%;"; context.drawImage(image, 0, 0, w, h); context.lineCap = "round"; context.lineJoin = "round"; context.lineWidth = 2; const originalImageData = context.getImageData(0, 0, w, h); const input = context.getImageData(0, 0, w, h); const inputData = input.data; const inputData32 = new Uint32Array(inputData.buffer); // In order to put the removed seams in the right locations // we have to store the original x-position of each pixel, // so that we can re-insert them in the right order const originalColumns = new Uint16Array(w * h); const newColumns = new Uint16Array(w * h); const removedSeamLoc = new Uint16Array(h); const insertedSeamLoc = new Uint16Array(h); function initBuffers() { inputData.set(originalImageData.data); for (let y = 0; y < h; y++) { for (let x = 0; x < w; x++) { originalColumns[x + y * w] = x; } } newColumns.fill(0); removedSeamLoc.fill(0); insertedSeamLoc.fill(0); } const { sqrt, min, max, SQRT1_2 } = Math; // =============== // = Pixel delta = // =============== const pixelOrder = (function checkEndian() { const arrayBuffer = new ArrayBuffer(2); const uint8Array = new Uint8Array(arrayBuffer); const uint16array = new Uint16Array(arrayBuffer); uint8Array[0] = 0xAA; // set first byte uint8Array[1] = 0xBB; // set second byte if (uint16array[0] === 0xBBAA) return "ABGR"; // little endian if (uint16array[0] === 0xAABB) return "RGBA"; // big endian else throw new Error("Something crazy just happened"); })(); // This is probably slow compared to inlining it, but let's get it working first const dPixel = pixelOrder == "RGBA" ? function dPixel(p0, p1) { const rgba0 = inputData32[p0]; const rgba1 = inputData32[p1]; const r0 = rgba0 >>> 24; const g0 = (rgba0 >> 16) & 0xFF; const b0 = (rgba0 >> 8) & 0xFF; const r1 = rgba1 >>> 24; const g1 = (rgba1 >> 16) & 0xFF; const b1 = (rgba1 >> 8) & 0xFF; // // Taking the difference between the squares // // is a better approximation of linear rgb // // while still being fairly cheap. // const dr = r0*r0 - r1*r1; // const dg = g0*g0 - g1*g1; // const db = b0*b0 - b1*b1; // return (dr < 0 ? -dr : dr) + (dg < 0 ? -dg : dg) + (db < 0 ? -db : db); // Adapted from code in https://github.com/mapbox/pixelmatch // calculate color difference according to the paper "Measuring perceived color difference // using YIQ NTSC transmission color space in mobile applications" by Y. Kotsarenko and F. Ramos const r = r0 - r1; const g = g0 - g1; const b = b0 - b1; const y = r * 0.2124681075446384 + g * 0.4169973963260294 + b * 0.08137907133969426; const i = r * 0.3258860837850668 - g * 0.14992193838645426 - b * 0.17596414539861255; const q = r * 0.0935501584120867 - g * 0.23119531908149002 + b * 0.13764516066940333; return sqrt(y * y + i * i + q * q); } : function dPixel(p0, p1) { const abgr0 = inputData32[p0]; const abgr1 = inputData32[p1]; const r0 = abgr0 & 0xFF; const g0 = (abgr0 >> 8) & 0xFF; const b0 = (abgr0 >> 16) & 0xFF; const r1 = abgr1 & 0xFF; const g1 = (abgr1 >> 8) & 0xFF; const b1 = (abgr1 >> 16) & 0xFF; // const dr = r0*r0 - r1*r1; // const dg = g0*g0 - g1*g1; // const db = b0*b0 - b1*b1; // return (dr < 0 ? -dr : dr) + (dg < 0 ? -dg : dg) + (db < 0 ? -db : db); const r = r0 - r1; const g = g0 - g1; const b = b0 - b1; const y = r * 0.2124681075446384 + g * 0.4169973963260294 + b * 0.08137907133969426; const i = r * 0.3258860837850668 - g * 0.14992193838645426 - b * 0.17596414539861255; const q = r * 0.0935501584120867 - g * 0.23119531908149002 + b * 0.13764516066940333; return sqrt(y * y + i * i + q * q); }; // ================= // = Cost Map Init = // ================= // build up energy map // see https://avikdas.com/2019/07/29/improved-seam-carving-with-forward-energy.html const costMap = new Float64Array(w * h); function initCostMap() { // energy cost pixel top-left corner costMap[0] = dPixel(0, 1); // energy cost top row (except corners for (let x = 1; x < w - 1; x++) { costMap[x] = dPixel(x - 1, x + 1); } // energy cost pixel top-right corner costMap[w - 1] = dPixel(w - 2, w - 1); // build cost map for rest of the rows for (let y = 1; y < h; y++) { const line = y * w; const linem1 = line - w; // leftmost column { // const x = 0; // we can simply leave out the zero altogether here const Cu = dPixel(line, line + 1) * 2; const Cr = dPixel(linem1, line + 1); const Mu = costMap[linem1] + Cu; const Mr = costMap[linem1 + 1] + Cu + Cr; costMap[line] = Mu < Mr ? Mu : Mr; } // rest of the row for (let x = 1; x < w - 1; x++) { const Cl = dPixel(x + linem1, x - 1 + line); const Cu = dPixel(x - 1 + line, x + 1 + line); const Cr = dPixel(x + linem1, x + 1 + line); // Diagonals. We'll weight them sqrt(1/2) compared to the rest const Clr = SQRT1_2 * dPixel(x - 1 + line, x + 1 + linem1); // bottom left to top right const Crl = SQRT1_2 * dPixel(x + 1 + line, x - 1 + linem1); // bottom right to top left const Ml = costMap[x - 1 + linem1] + Cu + Cl + Crl; const Mu = costMap[x + linem1] + Cu + Clr + Crl; const Mr = costMap[x + 1 + linem1] + Cu + Cr + Clr; const M = Mr < Mu ? Mr : Mu; costMap[x + line] = M < Ml ? M : Ml; } // rightmost column { const x = w - 1; const Cl = dPixel(x + linem1, x - 1 + line); const Cu = dPixel(x - 1 + line, x + line); const Ml = costMap[x - 1 + linem1] + Cu + Cl; const Mu = costMap[x + linem1] + Cu * 2; costMap[x + line] = Ml < Mu ? Ml : Mu; } } } // =============== // = Carve Loop = // =============== while (true) { initBuffers(); initCostMap(); let inputEnd = w - 1; while (inputEnd > 0) { // find minimum of the bottom row let y = h - 1; let line = y * w; let xMin = 0, costMin = costMap[line]; for (let x = 1; x < inputEnd; x++) { const xCost = costMap[x + line]; if (xCost < costMin) { costMin = xCost; xMin = x; } } // remove seam - go from bottom minimum and carve up while (y >= 0) { // copy pixeldata + original location at xMin; // moving 32 bits at once is faster than three 8-bit chunks const rgbMin = inputData32[xMin + line]; const loc = originalColumns[xMin + line]; removedSeamLoc[y] = xMin; // shift all remaining pixels (in "original" image) one to the left, // including their locations let i = xMin; let idx = i + line; while (i < inputEnd) { inputData32[idx] = inputData32[idx + 1]; originalColumns[idx] = originalColumns[idx + 1]; idx++; i++; } // insert removed pixel at the appropriate location on the right while (i < w - 1 && loc > newColumns[idx + 1]) { inputData32[idx] = inputData32[idx + 1]; newColumns[idx] = newColumns[idx + 1]; idx++; i++; } // const idx4 = idx << 2; inputData32[idx] = rgbMin; newColumns[idx] = loc; insertedSeamLoc[y] = i; // go up by one row y--; line = y * w; // The minimum cost is guaranteed to be one of the // three pixels right above the pixel with the // minimum cost of the row below it costMin = costMap[xMin + line]; if (costMap[max(0, xMin - 1) + line] < costMin) { xMin = max(0, xMin - 1); costMin = costMap[xMin + line]; } if (costMap[Math.min(inputEnd, xMin + 1) + line] < costMin) { xMin = Math.min(inputEnd, xMin + 1); costMin = costMap[xMin + line]; } } // Move the right border of the original (carved) image to the left by one inputEnd--; // Update image and draw seams context.putImageData(input, 0, 0); // context.strokeStyle = 'red'; // context.beginPath(); // context.moveTo(removedSeamLoc[0], 0); // for (let y = 1; y < h; y++) { // context.lineTo(removedSeamLoc[y], y); // } // context.stroke(); // context.strokeStyle = 'blue'; // context.beginPath(); // context.moveTo(insertedSeamLoc[0], 0); // for (let y = 1; y < h; y++) { // context.lineTo(insertedSeamLoc[y], y); // } // context.stroke(); yield { context, input, inputData32, inputEnd, originalColumns, newColumns }; // Update cost map // Optimization: we only have to update the pyramid // under wherever xMin ends up being on the top row, // because pixels outside of this pyramid cannot be // affected by the change. // Handle top-left and top-right corners if (xMin === 0) { costMap[0] = dPixel(0, 1); costMap[1] = dPixel(0, 2); } else if (xMin >= inputEnd) { costMap[inputEnd - 1] = dPixel(inputEnd - 2, inputEnd); costMap[inputEnd] = dPixel(inputEnd - 1, inputEnd); } else { costMap[xMin - 1] = dPixel(max(0, xMin - 2), xMin); costMap[xMin] = dPixel(xMin - 1, xMin + 1); costMap[xMin + 1] = dPixel(xMin, min(inputEnd, xMin + 2)); } // update cost map for rest of the rows let xStart = max(0, xMin - 1), xEnd = min(xMin + 1, inputEnd); for (let y = 1; y < h; y++) { const line = y * w; const linem1 = line - w; let x = xStart; // handle case where pyramid overlaps with start if (x === 0) { const Cu = dPixel(line, line + 1) * 2; const Cr = dPixel(linem1, line + 1); const Mu = costMap[linem1] + Cu; const Mr = costMap[linem1 + 1] + Cu + Cr; costMap[line] = Mu < Mr ? Mu : Mr; } else { const Cl = dPixel(x + linem1, x - 1 + line); const Cu = dPixel(x - 1 + line, x + 1 + line); const Cr = dPixel(x + linem1, x + 1 + line); // Diagonals. We'll weight them sqrt(1/2) compared to the rest const Clr = SQRT1_2 * dPixel(x - 1 + line, x + 1 + linem1); // bottom left to top right const Crl = SQRT1_2 * dPixel(x + 1 + line, x - 1 + linem1); // bottom right to top left const Ml = costMap[x - 1 + linem1] + Cu + Cl + Crl; const Mu = costMap[x + linem1] + Cu + Clr + Crl; const Mr = costMap[x + 1 + linem1] + Cu + Cr + Clr; const M = Mr < Mu ? Mr : Mu; costMap[x + line] = M < Ml ? M : Ml; } // inner pyramid while (++x < xEnd) { const Cl = dPixel(x + linem1, x - 1 + line); const Cu = dPixel(x - 1 + line, x + 1 + line); const Cr = dPixel(x + linem1, x + 1 + line); const Ml = costMap[x - 1 + linem1] + Cu + Cl; const Mu = costMap[x + linem1] + Cu; const Mr = costMap[x + 1 + linem1] + Cu + Cr; const M = Mr < Mu ? Mr : Mu; costMap[x + line] = M < Ml ? M : Ml; } // handle case where pyramid overlaps with end if (x === inputEnd) { const Cl = dPixel(x + linem1, x - 1 + line); const Cu = dPixel(x - 1 + line, x + line); const Ml = costMap[x - 1 + linem1] + Cu + Cl; const Mu = costMap[x + linem1] + Cu * 2; costMap[x + line] = Ml < Mu ? Ml : Mu; } else { const Cl = dPixel(x + linem1, x - 1 + line); const Cu = dPixel(x - 1 + line, x + 1 + line); const Cr = dPixel(x + linem1, x + 1 + line); const Ml = costMap[x - 1 + linem1] + Cu + Cl; const Mu = costMap[x + linem1] + Cu; const Mr = costMap[x + 1 + linem1] + Cu + Cr; const M = Mr < Mu ? Mr : Mu; costMap[x + line] = M < Ml ? M : Ml; } xStart = max(0, xStart - 1); xEnd = min(inputEnd, xEnd + 1); } } context.putImageData(input, 0, 0); yield { context, input, inputData32, inputEnd, originalColumns, newColumns }; } }