Seam Carving (with forward energy), but we make a new image with the discarded seams / Job van der Zwan

images = ({

"De Stellingwerven": {

image: await FileAttachment("IMG_20180615_134540_HDR-PANO.jpg").image(),

credit: md`*Photo Credit:* Job van der Zwan (me)`

},

"Blue Beach Surf": {

image: await FileAttachment("blue.jpg").image(),

credit: md`*Photo Credit:* [Blue Beach Surf](https://pixabay.com/photos/blue-beach-surf-travel-surfer-4145659/) by [Kiril Dobrev](https://pixabay.com/users/kirildobrev-12266114/)`

},

"North Windows Arch": {

image: await FileAttachment("North Windows Arch.jpg").image(),

credit: md`*Photo Credit:* [North Windows Arch, Arches National Park - 9/24/2007](https://flic.kr/p/4hxxz50) by

Mike Goad`

},

})

height = Math.round(width / image.naturalWidth * image.naturalHeight)

function createContext(img) {

Math.round(width / img.naturalWidth * img.naturalHeight)

const context = DOM.context2d(width, height, 1);

context.canvas.style = "max-width: 100%;";

context.drawImage(img, 0, 0, width, height);

return context;

}

loadImage = {

refresh;

return function loadImage(src) {

const img = DOM.element('img', {crossorigin: 'anonymous', src});

return img.decode().then(() => img);

};

}

// Note: Multiple requests made within a 3 seconds might return the same URL.

async function randomImage(sizeOrWidth = null, height = sizeOrWidth) {

const size = sizeOrWidth ? `/${+sizeOrWidth}x${+height}` : '';

return (await fetch(`https://source.unsplash.com/random${size}`, {

method: 'head',

cache: 'no-store',

})).url;

}

carved = {

refresh;

refresh2;

if (!image) {

yield md`*Awaiting image selection*`;

} else {

const context = createContext(image);

context.lineCap = "round";

context.lineJoin = "round";

context.lineWidth = 3;

const input = context.getImageData(0, 0, width, height);

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(width * height);

for (let y = 0; y < height; y++) {

for (let x = 0; x < width; x++) {

originalColumns[x + y * width] = x;

}

const newColumns = new Uint16Array(width * height);

const removedSeamLoc = new Uint16Array(height);

const insertedSeamLoc = new Uint16Array(height);

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(width * height);

// energy cost pixel top-left corner

costMap[0] = dPixel(0, 1);

// energy cost top row (except corners

for (let x = 1; x < width - 1; x++) {

costMap[x] = dPixel(x - 1, x + 1);

}

// energy cost pixel top-right corner

costMap[width - 1] = dPixel(width - 2, width - 1);

// build cost map for rest of the rows

for (let y = 1; y < height; y++) {

const line = y * width;

const linem1 = line - width;

// 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 < width - 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 = width - 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 =

// ===============

let inputEnd = width - 1;

while (inputEnd > width / 2) {

// find minimum of the bottom row

let y = height - 1;

let line = y * width;

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 < width - 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 * width;

// 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 < height; y++) {

context.lineTo(removedSeamLoc[y], y);

}

context.stroke();

context.strokeStyle = 'blue';

context.beginPath();

context.moveTo(insertedSeamLoc[0], 0);

for (let y = 1; y < height; 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 < height; y++) {

const line = y * width;

const linem1 = line - width;

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

};

}

// Restores (an approximation of) the original image from either the retargeted pixels or the discarded seams

function restoreOriginal(targetImage, columns, colStart, colEnd) {

const {inputData32, inputEnd} = carved;

const {context, output32, output} = targetImage;

for (let y = 0; y < height; y++) {

const line = y * width;

let x = 0;

for (let i = colStart; i < colEnd; i++) {

const xEnd = columns[i + line];

const rgba = inputData32[i + line];

while(x < xEnd) {

output32[x + line] = rgba;

x++;

}

const rgba = inputData32[inputEnd + line];

while (x < width) {

output32[x + line] = rgba;

x++;

}

context.putImageData(output, 0, 0);

}

// Same as above, but lerping

function restoreOriginalLerp(targetImage, columns, colStart, colEnd) {

const {

input: { data },

inputEnd

} = carved;

const { context, output32, output } = targetImage;

for (let y = 0; y < height; y++) {

const line = y * width;

const l4 = line * 4;

let x = 0;

let idx = (colStart + line) * 4;

let r0 = data[idx];

let g0 = data[idx + 1];

let b0 = data[idx + 2];

for (let i = colStart + 1; i < colEnd - 1; i++) {

const xEnd = columns[i + line];

idx = (i + line) * 4;

const r1 = data[idx];

const g1 = data[idx + 1];

const b1 = data[idx + 2];

const dx = xEnd - x;

for (let j = 0; x + j < xEnd; j++) {

const r = ((r0 * (dx - j) + r1 * j) / dx) & 0xFF;

const g = ((g0 * (dx - j) + g1 * j) / dx) & 0xFF;

const b = ((b0 * (dx - j) + b1 * j) / dx) & 0xFF;

output[j * 4 + idx] = r;

output[j * 4 + idx + 1] = g;

output[j * 4 + idx + 2] = b;

output[j * 4 + idx + 3] = 0xFF;

}

x = xEnd;

r0 = r1;

g0 = g1;

b0 = b1;

}

idx = (inputEnd + line) * 4;

const r1 = data[idx];

const g1 = data[idx + 1];

const b1 = data[idx + 2];

const dx = width - x;

for (let j = 0; j + x < width; j++) {

const r = ((r0 * (dx - j) + r1 * j) / dx) & 0xFF;

const g = ((g0 * (dx - j) + g1 * j) / dx) & 0xFF;

const b = ((b0 * (dx - j) + b1 * j) / dx) & 0xFF;

output[j * 4 + idx] = r;

output[j * 4 + idx + 1] = g;

output[j * 4 + idx + 2] = b;

output[j * 4 + idx + 3] = 0xFF;

}

context.putImageData(output, 0, 0);

}

// The restored image from the retargeted pixels

// Initialized in a separate cell to avoid pointless re-allocation

originalRestored = {

const context = createContext(image);

const output = context.createImageData(width, height);

const output32 = new Uint32Array(output.data.buffer);

return { context, output, output32 };

}

// The restored image from the discarded pixels

removedRestored = {

const context = createContext(image);

const output = context.createImageData(width, height);

const output32 = new Uint32Array(output.data.buffer);

return { context, output, output32 };

}

restoreOriginal(originalRestored, carved.originalColumns, 0, carved.inputEnd)

restoreOriginal(removedRestored, carved.newColumns, carved.inputEnd, width)

Purpose-built for displays of data