31:Tfc9,md`# Variations on Symmetric Nearest Neighbour smoothing [The 2.10 release of GIMP](https://www.gimp.org) added a *Symmetric Nearest Neighbour* filter:

I had never heard of this before, so expected it to be some novel fancy complicated algorithm, but it turns out to be quite simple *and* old, which made its obscurity surprising. Beyond the GIMP documentation itself, the only bit of useful information that I could find at first was [a brief but sufficient description on the SubSurf Wiki](http://subsurfwiki.org/wiki/Symmetric_nearest_neighbour_filter) ([archive.org mirror](https://web.archive.org/web/20190706044806/http://subsurfwiki.org/wiki/Symmetric_nearest_neighbour_filter)), a wiki for subsurface science. Looking for the sources mentioned there on Google Scholar showed that this filter was first described by Harwood et al. in a 1987 paper called [*"A new class of edge-preserving smoothing filters"*](https://scholar.google.se/scholar?cluster=11619776633817084182&hl=en&as_sdt=0,5). That makes it quite old, and even more surprising that it was not in any photo-editing software packages before. Perhaps it was encumbered by software patents? Or just somehow not picked up (I guess things were different in the pre-internet era). Whatever the cause, it's in GIMP now, and we can start playing around with it! ## How it works. The description of the filter [from gimp.org itself](https://pippin.gimp.org/image-processing/chap_area.html) is quite good: > *An improvement over the kuwahara filter is the symmetric nearest neighbour filter that compares symmetric pairs/quadruples of pixels with the center pixel and only takes into account the one from this set which is closest in value to the center pixel.* In other words: we blur, but only with the pixels most similar to the pixel being blurred over. The image below illustrates selecting between a symmetric pair of points:

Given a center point \`C\`, we look at symmetrically opposed pairs like \`P1\` and \`P2\`. The point that has the least color difference to \`C\` is selected (in this pairing that would be \`P2\`). This is done for the whole neighbourhood. Then the selected values will be averaged. The illustration shows the pair-wise selection method, but this notebook will mostly work with the cross-wise version, as it turns out to generally perform better visually speaking. Conceptually speaking this is quite a simple approach. One tricky aspect is defining a distance function for two colors (or, alternatively and not explored here, doing this separately for each color channel), but we can build on other people's work there. I decided to use the Kotsarenko/Ramos algorithm. ## Variations Thinking about the SNN algorithm some more, it becomes clear these are essentially two composed algorithms: - the pair-wise or cross-wise selection of pixels - averaging the selected pixels (which in the GIMP version as well as the original paper seems to be a simple box blur) That also means we can decouple the two, and experiment with using the first part as a form of "preselection" with a whole bunch of other image filters and see what the results of that will be! It turns out that the authors of the original paper already realised this, and created both a "SNN mean" and "SNN median" variant, effectively composing the selection algorithm with a box blur and a median blur, respectively. But let's try replacing the box blur with something generally considered a better image filter, like a [Gaussian blur](https://en.wikipedia.org/wiki/Gaussian_blur), see what the effects are on the image quality, and explore from there. As test images I picked two photographs I took about a decade ago: `32:T7a1,// parameters are passed as an object, so that we can turn boxblur into a worker function boxblur({source, width, height, radius}) { if (!(radius > 0)) return source.slice(); let target = new Uint8ClampedArray(source.length); /* * total number of summed pixels is two times the radius plus * the central pixel, squared. See for example this neighbourhood * of radius 3 around a pixel: * * . . . . . . . * . . . . . . . * . . . . . . . * . . . # . . . * . . . . . . . * . . . . . . . * . . . . . . . * * At the end we have to divide by the total area, but it's quicker * to multiply by the inverse instead. */ const factor = 1 / ((2*radius + 1) * (2*radius + 1)); const {min, max} = Math; /* * Tangent: this is a really slow way to implement a blur function * We loop over the entire area for every pixel, which makes it * effectively O(n * r²). That gets slow really quickly for larger * radii (not to mention memory access patterns and such). There * are faster ways methods - we explore a two-pass version below * that is effectively 2 * O(n*r) near the end of this notebook. * For educational purposes we stick to these slow versions first. */ for(let y = 0; y < height; y++) { const line = y * width * 4; for (let x = 0; x < width; x++) { let r = 0; let g = 0; let b = 0; for (let dy = -radius; dy <= radius; dy++) { const line_dy = min(max(0, y + dy), height-1) * width; for (let dx = -radius; dx <= radius; dx++) { const idx = (min(max(0, x + dx), width-1) + line_dy) * 4; r += source[idx]; g += source[idx + 1]; b += source[idx + 2]; } } const cIdx = x*4 + line; target[cIdx] = r * factor; target[cIdx+1] = g * factor; target[cIdx+2] = b * factor; target[cIdx+3] = 255; } } return target; }33:T63d,// This is not very efficient, so don't use this in filters with big // radii. On the other hand, it's easier to follow code-wise function selectNearest(cx, cy, dx, dy, source, width, height) { // assuming that source is an RGBA array, the central index should be at: const cIdx = (cx + cy * width) * 4; // Side-note: interesting dilemma: do we pick "per color" or "per color channel"? // I suspect the latter will perform better. However, we first try this with a // closest-color formula using Kotsarenko/Ramos. // 35215 is the maximum possible value for the Kotsarenko/Ramos YIQ difference metric let minDiff = 35215 ; let pick = cIdx; // we only want to test pixels that are inside the boundaries if (cx + dx < width && cy + dy < height) { const pos = cIdx + (dx + dy * width) * 4; const diff = colorDelta(source, cIdx, pos); if (diff < minDiff) { pick = pos; minDiff = diff; } } if (cx - dx >= 0 && cy - dy >= 0) { const pos = cIdx - (dx + dy * width) * 4; const diff = colorDelta(source, cIdx, pos); if (diff < minDiff) { pick = pos; minDiff = diff; } } if (cx + dy < width && cy - dx >= 0) { const pos = cIdx + (dy - dy * width) * 4; const diff = colorDelta(source, cIdx, pos); if (diff < minDiff) { pick = pos; minDiff = diff; } } if (cx - dy >= 0 && cy + dx < height) { const pos = cIdx - (dy - dy * width) * 4; const diff = colorDelta(source, cIdx, pos); if (diff < minDiff) { pick = pos; minDiff = diff; } } return pick; }34:T911,// Adapted from code in https://github.com/mapbox/pixelmatch /* Inlined for better performance // function rgb2y(r, g, b) { return r * 0.29889531 + g * 0.58662247 + b * 0.11448223; } // function rgb2i(r, g, b) { return r * 0.59597799 - g * 0.27417610 - b * 0.32180189; } // function rgb2q(r, g, b) { return r * 0.21147017 - g * 0.52261711 + b * 0.31114694; } // // blend semi-transparent color with white // function blend(c, a) { // return 255 + (c - 255) * a; // } */ // 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 function colorDelta(img, pos1, pos2) { let r1 = img[pos1 + 0]; let g1 = img[pos1 + 1]; let b1 = img[pos1 + 2]; // let a1 = img[pos1 + 3]; // we ignore alpha in this demo, for perf reasons let r2 = img[pos2 + 0]; let g2 = img[pos2 + 1]; let b2 = img[pos2 + 2]; // let a2 = img[pos2 + 3]; if (r1 === r2 && g1 === g2 && b1 === b2) return 0; // if (a1 < 255) { // a1 /= 255; // r1 = 255 + (r1 - 255) * a1; // blend(r1, a1); // g1 = 255 + (g1 - 255) * a1; // blend(g1, a1); // b1 = 255 + (b1 - 255) * a1; // blend(b1, a1); // } // if (a2 < 255) { // a2 /= 255; // r2 = 255 + (r2 - 255) * a2; // blend(r2, a2); // g2 = 255 + (g2 - 255) * a2; // blend(g2, a2); // b2 = 255 + (b2 - 255) * a2; // blend(b2, a2); // } // const y = rgb2y(r1, g1, b1) - rgb2y(r2, g2, b2); // const i = rgb2i(r1, g1, b1) - rgb2i(r2, g2, b2); // const q = rgb2q(r1, g1, b1) - rgb2q(r2, g2, b2); const r = r1 - r2; const g = g1 - g2; const b = b1 - b2; // Without strength-reduction // const y = r * 0.29889531 + g * 0.58662247 + b * 0.11448223 // const i = r * 0.59597799 - g * 0.27417610 - b * 0.32180189 // const q = r * 0.21147017 - g * 0.52261711 + b * 0.31114694 // return 0.5053 * y * y + 0.299 * i * i + 0.1957 * q * q; // With strength-reduced constants 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 y*y + i*i + q*q }35:T4ba,md`And with a simple tweak we created a new filter: Gaussian Symmetric Nearest Neighbour Blurring. There are definitely some differences with the previous filter, although it might not be very obvious if you don't see the pictures side-by-side. Specifically, it looks like our edge-effects are less severe than before: the leaves don't have a darkened center, and the bright spot in the second photo remains bright. This makes sense: pixels further away from the central pixel have a lower weight, so even if our radius might be too large, the central pixels are weighed more heavily. As a result, the inverted-center effect that we had with the unweighted box blur is significantly reduced. Another thing to notice is that the image is less smoothed than before for the same radius value. The dust in the top-left corner isn't blurred away as strongly as it was with the SNN box blur. This is again a result of central pixels having more weight than outer ones. Still, in overall image quality I think this wins over our previous filter. ## Symmetrical *Furthest* Neighbour? So far we have been picking the colors that are the most similar. What would be the effect of picking the *least* similar color?`36:T555,function boxblurSNN({source, width, height, radius}) { if (!(radius > 0)) return source.slice(); let target = new Uint8ClampedArray(source.length); /* * total number of summed pixels is the central pixel, * plus one quadrant, which turns out to be radius * (radius + 1) * See for example this neighbourhood of radius 3 around a pixel: * * | | | | - - - * | | | | - - - * | | | | - - - * - - - # - - - * - - - | | | | * - - - | | | | * - - - | | | | * At the end we have to divide by that number, but it's quicker * to multiply by the inverse instead. */ const factor = 1 / (radius * (radius + 1) + 1); for(let y = 0; y < height; y++) { const line = y * width * 4; for (let x = 0; x < width; x++) { const cIdx = x*4 + line; let r = source[cIdx]; let g = source[cIdx + 1]; let b = source[cIdx + 2]; for (let dx = 1; dx <= radius; dx++) { for (let dy = 0; dy <= radius; dy++) { const pick = selectNearest(x, y, dx, dy, source, width, height); r += source[pick]; g += source[pick + 1]; b += source[pick + 2]; } } target[cIdx] = r * factor | 0; target[cIdx+1] = g * factor | 0; target[cIdx+2] = b * factor | 0; target[cIdx+3] = 255; } } return target; }37:T916,md`And just like that our general blur turns into a filter that both smooths and (somewhat) enhances edges! For example, the dirt in the upper left of the tulip photo is mostly smoothed away. At the same time, the edge of the outer ring of the fisheye lens is more defined. For other examples of the edge-enhancing effect, look at the blurry buildings in the horizon. While there is no fine detail to recover, what used to be a blurry transition to the sky changed into a sharp edge. Similarly, the couple in the second photo has a sharper silhouette, while the overal amount of grain is greatly reduced. However, it's not all good. If you look closely there are some strange "inverted area" effects happening. For example, some of the small leaves in the bottom background of the first photo: their edges are enhanced, but their centers are actually blurred with the pixels *outside* the leaf. A similar effect can be seen with the bright spot in the upper-left of the kissing couple. It is much more muted than before. What causes these "inversion" effects? Well, I set the blur radius to 10 pixels, so the total pixel area being averaged over is 21 by 21 pixels. If a leaf has a smaller area than that, say 15 by 15 pixels, then a pixel in the middle of that leaf will have many instances where all four selectable neighbours are outside of the leaf. Meanwhile, for a pixel on the edge of the leaf there will always be one selectable neighbour inside the leaf. So counter-intuitively, pixels at the edge of an area have a better chance of only selecting neighbours inside that area! So that is our basic SNN filter. Let's see what else we can build with this! Maybe we can even deal with the problem described above? Let's "invent" a Gaussian Symmetric Nearest Neighbour Blur. A Gaussian blur is basically just a blur with weights. So all we need to do is determine a Gaussian kernel of weights and add it into our existing algorithms. Although [I wrote a notebook for determining a proper Gaussian kernel a while ago](https://observablehq.com/@jobleonard/gaussian-kernel-calculater), I don't really know how to automatically find the right sigma-value for a given radius. So instead we'll use an approximation based on [the quadratic stackblur](https://observablehq.com/@jobleonard/mario-klingemans-stackblur) method:`38:T5cb,function gaussianBlurSNN({source, width, height, radius}) { if (!(radius > 0)) return source.slice(); let target = new Uint8ClampedArray(source.length); let totalWeight = 0; const kernel = bellcurvish2D(radius); // because we're dealing with quadrants again, we only want // 1/4 of the kernel weights plus the central pixel. for(let x = radius ; x < kernel.length; x++) { const row = kernel[x]; for (let y = radius + 1; y < row.length; y++) { totalWeight += row[y]; } } const centerWeight = kernel[radius][radius]; totalWeight += centerWeight; const norm = 1 / totalWeight const {min, max} = Math; for(let y = 0; y < height; y++) { const line = y * width * 4; for (let x = 0; x < width; x++) { const cIdx = x*4 + line; let r = source[cIdx] * centerWeight; let g = source[cIdx + 1] * centerWeight; let b = source[cIdx + 2] * centerWeight; for (let dx = 1; dx <= radius; dx++) { const row = kernel[radius + dx]; for (let dy = 0; dy <= radius; dy++) { const pick = selectNearest(x, y, dx, dy, source, width, height); const weight = row[radius + dy]; r += source[pick] * weight; g += source[pick + 1] * weight; b += source[pick + 2] * weight; } } target[cIdx] = r * norm | 0; target[cIdx+1] = g * norm | 0; target[cIdx+2] = b * norm | 0; target[cIdx+3] = 255; } } return target; }39:T5cc,function gaussianBlurSFN({source, width, height, radius}) { if (!(radius > 0)) return source.slice(); let target = new Uint8ClampedArray(source.length); let totalWeight = 0; const kernel = bellcurvish2D(radius); // because we're dealing with quadrants again, we only want // 1/4 of the kernel weights plus the central pixel. for(let x = radius ; x < kernel.length; x++) { const row = kernel[x]; for (let y = radius + 1; y < row.length; y++) { totalWeight += row[y]; } } const centerWeight = kernel[radius][radius]; totalWeight += centerWeight; const norm = 1 / totalWeight const {min, max} = Math; for(let y = 0; y < height; y++) { const line = y * width * 4; for (let x = 0; x < width; x++) { const cIdx = x*4 + line; let r = source[cIdx] * centerWeight; let g = source[cIdx + 1] * centerWeight; let b = source[cIdx + 2] * centerWeight; for (let dx = 1; dx <= radius; dx++) { const row = kernel[radius + dx]; for (let dy = 0; dy <= radius; dy++) { const pick = selectFurthest(x, y, dx, dy, source, width, height); const weight = row[radius + dy]; r += source[pick] * weight; g += source[pick + 1] * weight; b += source[pick + 2] * weight; } } target[cIdx] = r * norm | 0; target[cIdx+1] = g * norm | 0; target[cIdx+2] = b * norm | 0; target[cIdx+3] = 255; } } return target; }3a:T656,function selectDipole(cx, cy, dx, dy, source, width, height) { // assuming that source is an RGBA array, the central index should be at: const cIdx = (cx + cy * width) * 4; // 35215 is the maximum possible value for the Kotsarenko/Ramos YIQ difference metric let minDiff = 35215, maxDiff = 0; let minPick = cIdx, maxPick = cIdx; // we only want to test pixels that are within the boundaries if (cx + dx < width && cy + dy < height) { const pos = cIdx + (dx + dy * width) * 4; const diff = colorDelta(source, cIdx, pos); if (diff < minDiff) { minPick = pos; minDiff = diff; } if (diff > maxDiff) { maxPick = pos; maxDiff = diff; } } if (cx - dx >= 0 && cy - dy >= 0) { const pos = cIdx - (dx + dy * width) * 4; const diff = colorDelta(source, cIdx, pos); if (diff < minDiff) { minPick = pos; minDiff = diff; } if (diff > maxDiff) { maxPick = pos; maxDiff = diff; } } if (cx + dy < width && cy - dx >= 0) { const pos = cIdx + (dy - dy * width) * 4; const diff = colorDelta(source, cIdx, pos); if (diff < minDiff) { minPick = pos; minDiff = diff; } if (diff > maxDiff) { maxPick = pos; maxDiff = diff; } } if (cx - dy >= 0 && cy + dx < height) { const pos = cIdx - (dy - dy * width) * 4; const diff = colorDelta(source, cIdx, pos); if (diff < minDiff) { minPick = pos; minDiff = diff; } if (diff > maxDiff) { maxPick = pos; maxDiff = diff; } } return {minPick, minDiff, maxPick, maxDiff}; }3b:T775,function symmetricalNeighbourEdgeDetection({source, width, height, radius}){ if (!(radius > 0)) return source.slice(); let target = new Uint8ClampedArray(source.length); let totalWeight = 0; const kernel = bellcurvish2D(radius); // because we're dealing with quadrants again, we only want // 1/4 of the kernel weights plus the central pixel. for(let x = radius ; x < kernel.length; x++) { const row = kernel[x]; for (let y = radius + 1; y < row.length; y++) { totalWeight += row[y]; } } const centerWeight = kernel[radius][radius]; totalWeight += centerWeight; const norm = 1 / totalWeight const {min, max, abs, round} = Math; for(let y = 0; y < height; y++) { const line = y * width * 4; for (let x = 0; x < width; x++) { const cIdx = x*4 + line; let rf = source[cIdx] * centerWeight; let gf = source[cIdx + 1] * centerWeight; let bf = source[cIdx + 2] * centerWeight; let rn = rf; let gn = gf; let bn = bf; for (let dx = 1; dx <= radius; dx++) { const row = kernel[radius + dx]; for (let dy = 0; dy <= radius; dy++) { const {minPick, maxPick} = selectDipole(x, y, dx, dy, source, width, height); const weight = row[radius + dy]; rf += source[maxPick] * weight; gf += source[maxPick + 1] * weight; bf += source[maxPick + 2] * weight; rn += source[minPick] * weight; gn += source[minPick + 1] * weight; bn += source[minPick + 2] * weight; } } // In order to show whether the edge is brighter or darker, // we subtract the difference from 127.5, or neutral grey target[cIdx] = round(127.5 - (rn - rf) * norm); target[cIdx+1] = round(127.5 - (gn - gf) * norm); target[cIdx+2] = round(127.5 - (bn - bf) * norm); target[cIdx+3] = 255; } } return target; }3c:T439,function selectFurthest(cx, cy, dx, dy, source, width, height) { const cIdx = (cx + cy * width) * 4; let maxDiff = 0 ; let pick = cIdx; // we only want to test pixels that are inside the boundaries if (cx + dx < width && cy + dy < height) { const pos = cIdx + (dx + dy * width) * 4; const diff = colorDelta(source, cIdx, pos); if (diff > maxDiff) { pick = pos; maxDiff = diff; } } if (cx - dx >= 0 && cy - dy >= 0) { const pos = cIdx - (dx + dy * width) * 4; const diff = colorDelta(source, cIdx, pos); if (diff > maxDiff) { pick = pos; maxDiff = diff; } } if (cx + dy < width && cy - dx >= 0) { const pos = cIdx + (dy - dy * width) * 4; const diff = colorDelta(source, cIdx, pos); if (diff > maxDiff) { pick = pos; maxDiff = diff; } } if (cx - dy >= 0 && cy + dx < height) { const pos = cIdx - (dy - dy * width) * 4; const diff = colorDelta(source, cIdx, pos); if (diff > maxDiff) { pick = pos; maxDiff = diff; } } return pick; }3d:T768,// A version that shows the absolute value of the change - highlights which color is // close to an edge without including "direction" of change function symmetricalNeighbourEdgeDetectionAbsolute({source, width, height, radius}){ if (!(radius > 0)) return source.slice(); let target = new Uint8ClampedArray(source.length); let totalWeight = 0; const kernel = bellcurvish2D(radius); // because we're dealing with quadrants again, we only want // 1/4 of the kernel weights plus the central pixel. for(let x = radius ; x < kernel.length; x++) { const row = kernel[x]; for (let y = radius + 1; y < row.length; y++) { totalWeight += row[y]; } } const centerWeight = kernel[radius][radius]; totalWeight += centerWeight; const norm = 1 / totalWeight const {min, max, abs, round} = Math; for(let y = 0; y < height; y++) { const line = y * width * 4; for (let x = 0; x < width; x++) { const cIdx = x*4 + line; let rf = source[cIdx] * centerWeight; let gf = source[cIdx + 1] * centerWeight; let bf = source[cIdx + 2] * centerWeight; let rn = rf; let gn = gf; let bn = bf; for (let dx = 1; dx <= radius; dx++) { const row = kernel[radius + dx]; for (let dy = 0; dy <= radius; dy++) { const {minPick, maxPick} = selectDipole(x, y, dx, dy, source, width, height); const weight = row[radius + dy]; rf += source[maxPick] * weight; gf += source[maxPick + 1] * weight; bf += source[maxPick + 2] * weight; rn += source[minPick] * weight; gn += source[minPick + 1] * weight; bn += source[minPick + 2] * weight; } } target[cIdx] = abs(rn - rf) * norm; target[cIdx+1] = abs(gn - gf) * norm; target[cIdx+2] = abs(bn - bf) * norm; target[cIdx+3] = 255; } } return target; }3e:T5ed,function unsharpMask({source, width, height, radius}) { if (!(radius > 0)) return source.slice(); let target = new Uint8ClampedArray(source.length); let totalWeight = 0; const kernel = bellcurvish2D(radius); // we need a negative weight to get an unsharp-mask effect for(let i = 0; i < kernel.length; i++) { const row = kernel[i]; for(let j = 0; j < kernel.length; j++) { row[j] = -row[j]; } } // the central kernel has to be offset by the negative weight const centerWeight = 2 + kernel[radius][radius]; kernel[radius][radius] = centerWeight; const {min, max} = Math; for(let y = 0; y < height; y++) { const line = y * width * 4; for (let x = 0; x < width; x++) { let r = 0; let g = 0; let b = 0; for (let dy = -radius; dy <= radius; dy++) { const line_dy = min(max(0, y + dy), height-1) * width; // funny enough, because our kernel is symmetrical we can // interchangably do kernel[x][y] or kernel[y][x] const column = kernel[radius + dy]; for (let dx = -radius; dx <= radius; dx++) { const idx = (min(max(0, x + dx), width-1) + line_dy) * 4; const weight = column[radius + dx]; r += source[idx] * weight; g += source[idx + 1] * weight b += source[idx + 2] * weight; } } const cIdx = x*4 + line; target[cIdx] = r; target[cIdx+1] = g; target[cIdx+2] = b; target[cIdx+3] = 255; } } return target; }3f:T780,function symmetricalNeighbourSharpen({source, width, height, radius}){ if (!(radius > 0)) return source.slice(); let target = new Uint8ClampedArray(source.length); let totalWeight = 0; const kernel = bellcurvish2D(radius); // because we're dealing with quadrants again, we only want // 1/4 of the kernel weights plus the central pixel. for(let x = radius ; x < kernel.length; x++) { const row = kernel[x]; for (let y = radius + 1; y < row.length; y++) { totalWeight += row[y]; } } const centerWeight = kernel[radius][radius]; totalWeight += centerWeight; const norm = 1 / totalWeight const {min, max, abs, round} = Math; for(let y = 0; y < height; y++) { const line = y * width * 4; for (let x = 0; x < width; x++) { const cIdx = x*4 + line; const r0 = source[cIdx]; const g0 = source[cIdx + 1]; const b0 = source[cIdx + 2]; let rf = r0 * centerWeight; let gf = g0 * centerWeight; let bf = b0 * centerWeight; let rn = rf; let gn = gf; let bn = bf; for (let dx = 1; dx <= radius; dx++) { const row = kernel[radius + dx]; for (let dy = 0; dy <= radius; dy++) { const {minPick, maxPick} = selectDipole(x, y, dx, dy, source, width, height); const weight = row[radius + dy]; rf += source[maxPick] * weight; gf += source[maxPick + 1] * weight; bf += source[maxPick + 2] * weight; rn += source[minPick] * weight; gn += source[minPick + 1] * weight; bn += source[minPick + 2] * weight; } } // instead of neutral grey, we now apply // the difference to the original pixel target[cIdx] = round(r0 + (rn - rf) * norm); target[cIdx+1] = round(g0 + (gn - gf) * norm); target[cIdx+2] = round(b0 + (bn - bf) * norm); target[cIdx+3] = 255; } } return target; }40:T444,md`And out pops a sharpening filter. The halos are even more pronounced than with the regular unsharp mask though. We'll try to fix that in a bit, but notice too that there are some enhancements: with the plain unsharp mask filter, the grain of the kissing couple photo got more pronounced (more extreme black and white pixels). This version of sharpening seems to avoid that problem, while still keeping the grain. So about those halos. One way to think about what causes this is that we push the pixels away from the furthest neighbours, and pull them towards the closest neighbours. However, we do this push by subtracting the differences. Which means that bigger differences will have a proportionally larger effect. If you then add the realization that the *furthest* neighbour will likely have a bigger difference than the *nearest* neighbour, kind of by definition, it should not be surprising any more that we end up overshooting our target. But that also means we have an easy fix: we have to scale the "push" of the furthest difference down to the "pull" of the nearest neighbour!`41:T8f7,function symmetricalNeighbourSharpenNoHalo({source, width, height, radius}){ if (!(radius > 0)) return source.slice(); let target = new Uint8ClampedArray(source.length); let totalWeight = 0; const kernel = bellcurvish2D(radius); // because we're dealing with quadrants again, we only want // 1/4 of the kernel weights plus the central pixel. for(let x = radius ; x < kernel.length; x++) { const row = kernel[x]; for (let y = radius + 1; y < row.length; y++) { totalWeight += row[y]; } } const centerWeight = kernel[radius][radius]; totalWeight += centerWeight; const norm = 1 / totalWeight const {min, max, abs, round} = Math; for(let y = 0; y < height; y++) { const line = y * width * 4; for (let x = 0; x < width; x++) { const cIdx = x*4 + line; const r0 = source[cIdx]; const g0 = source[cIdx + 1]; const b0 = source[cIdx + 2]; let rf = r0 * centerWeight; let gf = g0 * centerWeight; let bf = b0 * centerWeight; let rn = rf; let gn = gf; let bn = bf; let diffF = 0; let diffN = 0; for (let dx = 1; dx <= radius; dx++) { const row = kernel[radius + dx]; for (let dy = 0; dy <= radius; dy++) { const {minPick, minDiff, maxPick, maxDiff} = selectDipole(x, y, dx, dy, source, width, height); diffF += maxDiff; diffN += minDiff; const weight = row[radius + dy]; rf += source[maxPick] * weight; gf += source[maxPick + 1] * weight; bf += source[maxPick + 2] * weight; rn += source[minPick] * weight; gn += source[minPick + 1] * weight; bn += source[minPick + 2] * weight; } } const diffscale = (diffN > 0 && diffF > 0) ? diffN / diffF : 1; rf = (rf * norm - r0) * diffscale; gf = (gf * norm - g0) * diffscale; bf = (bf * norm - b0) * diffscale; rn = rn * norm - r0; gn = gn * norm - g0; bn = bn * norm - b0; // instead of neutral grey, we now apply // the difference to the original pixel target[cIdx] = round(r0 + (rn - rf)); target[cIdx+1] = round(g0 + (gn - gf)); target[cIdx+2] = round(b0 + (bn - bf)); target[cIdx+3] = 255; } } return target; }42:T402,md`## Two-Pass Symmetric Nearest Neighbour Filters? Regular blur functions have the annoying property that a single-pass filter have to average over a surface area, which scales badly. To be exact: if our blurring radius is ${tex`n`} pixels, our total number of pixels being averaged over is actually ${tex`(2n+1)^2`} pixels. One optimization is to blur in two passes instead: one horizontal, and then one vertical. Because blurring is [commutative](https://en.wikipedia.org/wiki/Commutative_property) as long as we avoid severe rounding errors, this is visually almost the same, but instead only needs to process ${tex`2(2n+1)`} pixels per pixel. However, applying the SNN/SFN approaches in this case would not work in the exact same way: we would have to first apply it to the horizontal blur, and then to the vertical blur. I have no idea what the expected effects on the results would be, so let's see what it does! I'll skip rendering the plain blur versions, because it has (almost) identical output to the originals.`43:T50e,function boxblurTwoPass({source, width, height, radius}) { let target = new Uint8ClampedArray(source.length), buffer = new Uint32Array(source.length); const factor = 1 / ((2*radius + 1)*(2*radius + 1)); const {min, max} = Math; // Horizontal pass for(let y = 0; y < height; y++) { const line = y * width * 4; for (let x = 0; x < width; x++) { let r = 0; let g = 0; let b = 0; for (let dx = -radius; dx <= radius; dx++) { const idx = min(max(0, x + dx), width-1) * 4 + line; r += source[idx]; g += source[idx + 1]; b += source[idx + 2]; } const cIdx = x*4 + line; buffer[cIdx] = r; buffer[cIdx+1] = g; buffer[cIdx+2] = b; } } for(let y = 0; y < height; y++) { for (let x = 0; x < width; x++) { let r = 0; let g = 0; let b = 0; const cIdx = (x + y * width) * 4; for (let dy = -radius; dy <= radius; dy++) { const idx = (x + min(max(0, y + dy), height-1)*width) * 4; r += buffer[idx]; g += buffer[idx + 1]; b += buffer[idx + 2]; } target[cIdx] = r * factor | 0; target[cIdx+1] = g * factor | 0; target[cIdx+2] = b * factor | 0; target[cIdx+3] = 255; } } return target; }44:T56a,function boxblurSNNTwoPass({source, width, height, radius}) { let target = new Uint8ClampedArray(source.length), buffer = new Float64Array(source.length); const factor = 1 / (radius + 1); const {min, max} = Math; // Horizontal pass for(let y = 0; y < height; y++) { const line = y * width * 4; for (let x = 0; x < width; x++) { const cIdx = x*4 + line; let r = source[cIdx]; let g = source[cIdx]; let b = source[cIdx]; for (let dx = 1; dx <= radius; dx++) { const pick = selectNearest(x, y, dx, 0, source, width, height); r += source[pick]; g += source[pick + 1]; b += source[pick + 2]; } buffer[cIdx] = r * factor; buffer[cIdx+1] = g * factor; buffer[cIdx+2] = b * factor; } } for(let y = 0; y < height; y++) { for (let x = 0; x < width; x++) { const cIdx = (x + y * width) * 4; let r = buffer[cIdx]; let g = buffer[cIdx+1]; let b = buffer[cIdx+2]; for (let dy = 1; dy <= radius; dy++) { const pick = selectNearest(x, y, 0, dy, buffer, width, height); r += buffer[pick]; g += buffer[pick + 1]; b += buffer[pick + 2]; } target[cIdx] = r * factor | 0; target[cIdx+1] = g * factor | 0; target[cIdx+2] = b * factor | 0; target[cIdx+3] = 255; } } return target; }45:T71c,function gaussianBlurSNNTwoPass({source, width, height, radius}) { let target = new Uint8ClampedArray(source.length), buffer = new Float64Array(source.length); const kernel = bellcurvish1D(radius); let totalWeight = 0; for(let i = radius; i < kernel.length; i++) totalWeight += kernel[i]; const norm = 1 / totalWeight; const {min, max} = Math; // Horizontal pass const centerWeight = kernel[radius]; for(let y = 0; y < height; y++) { const line = y * width * 4; for (let x = 0; x < width; x++) { const cIdx = x*4 + line; let r = source[cIdx] * centerWeight; let g = source[cIdx+1] * centerWeight; let b = source[cIdx+2] * centerWeight; for (let dx = 1; dx <= radius; dx++) { const pick = selectNearest(x, y, dx, 0, source, width, height); const weight = kernel[radius + dx] r += source[pick] * weight; g += source[pick + 1] * weight; b += source[pick + 2] * weight; } buffer[cIdx] = r * norm; buffer[cIdx+1] = g * norm; buffer[cIdx+2] = b * norm; } } // Vertical pass for(let y = 0; y < height; y++) { const line = y * width * 4; for (let x = 0; x < width; x++) { const cIdx = x*4 + line; let r = buffer[cIdx] * centerWeight; let g = buffer[cIdx+1] * centerWeight; let b = buffer[cIdx+2] * centerWeight; for (let dy = 1; dy <= radius; dy++) { const pick = selectNearest(x, y, 0, dy, buffer, width, height); const weight = kernel[radius + dy]; r += buffer[pick] * weight; g += buffer[pick + 1] * weight; b += buffer[pick + 2] * weight; } target[cIdx] = r * norm; target[cIdx+1] = g * norm; target[cIdx+2] = b * norm; target[cIdx+3] = 255; } } return target; }46:Td4b,// This is becoming quite the mouthful function symmetricalNeighbourSharpenNoHaloTwoPass({source, width, height, radius}){ if (!(radius > 0)) return source.slice(); let target = new Uint8ClampedArray(source.length), buffer = new Float64Array(source.length); const kernel = bellcurvish1D(radius); let totalWeight = 0; const centerWeight = kernel[radius]; totalWeight += centerWeight; for(let i = radius+1; i < kernel.length; i++) { // we apply the unsharp mask in two passes, so we should // reduce the weights a bit in order to avoid double-sharpening const v = kernel[i] * 0.5; totalWeight += v; kernel[i] = v; } const norm = 1 / totalWeight; const {min, max, abs, round} = Math; // Horizontal pass for(let y = 0; y < height; y++) { const line = y * width * 4; for (let x = 0; x < width; x++) { const cIdx = x*4 + line; const r0 = source[cIdx]; const g0 = source[cIdx + 1]; const b0 = source[cIdx + 2]; let rf = r0 * centerWeight; let gf = g0 * centerWeight; let bf = b0 * centerWeight; let rn = rf; let gn = gf; let bn = bf; let diffF = 0; let diffN = 0; for (let dx = 1; dx <= radius; dx++) { const {minPick, minDiff, maxPick, maxDiff} = selectDipole(x, y, dx, 0, source, width, height); diffF += maxDiff; diffN += minDiff; const weight = kernel[radius + dx]; rf += source[maxPick] * weight; gf += source[maxPick + 1] * weight; bf += source[maxPick + 2] * weight; rn += source[minPick] * weight; gn += source[minPick + 1] * weight; bn += source[minPick + 2] * weight; } const diffscale = (diffN > 0 && diffF > 0) ? diffN / diffF : 1; rf = (rf * norm - r0) * diffscale; gf = (gf * norm - g0) * diffscale; bf = (bf * norm - b0) * diffscale; rn = rn * norm - r0; gn = gn * norm - g0; bn = bn * norm - b0; buffer[cIdx] = rn - rf; buffer[cIdx+1] = gn - gf; buffer[cIdx+2] = bn - bf; } } // Vertical pass for(let y = 0; y < height; y++) { const line = y * width * 4; for (let x = 0; x < width; x++) { const cIdx = x*4 + line; const r0 = source[cIdx]; const g0 = source[cIdx + 1]; const b0 = source[cIdx + 2]; let rf = buffer[cIdx] * centerWeight; let gf = buffer[cIdx + 1] * centerWeight; let bf = buffer[cIdx + 2] * centerWeight; let rn = rf; let gn = gf; let bn = bf; let diffF = 0; let diffN = 0; for (let dy = 1; dy <= radius; dy++) { const {minPick, minDiff, maxPick, maxDiff} = selectDipole(x, y, dy, 0, source, width, height); diffF += maxDiff; diffN += minDiff; const weight = kernel[radius + dy]; rf += buffer[maxPick] * weight; gf += buffer[maxPick + 1] * weight; bf += buffer[maxPick + 2] * weight; rn += buffer[minPick] * weight; gn += buffer[minPick + 1] * weight; bn += buffer[minPick + 2] * weight; } const diffscale = (diffN > 0 && diffF > 0) ? diffN / diffF : 1; target[cIdx] = round(r0 + (rn - rf * diffscale) * norm); target[cIdx+1] = round(g0 + (gn - gf * diffscale) * norm); target[cIdx+2] = round(b0 + (bn - bf * diffscale) * norm); target[cIdx+3] = 255; } } return target; }47:T4ab,md`There might be some way to formulate the ideal sharpening strength automatically, but right now I don't see it. ## Conclusions So, as we've seen, with the simple addition of a pixel selection method (nearest or furthest neighbour in the color space) we can turn simple blurring filters into edge-enhancing smoothing filters. One issue with the original algorithm is the area-inversion effect when the radius covers larger areas than the fine detail. By using a Gaussian blur kernel we can partially address this issue. Furthermore, by combining nearest and furthest neighbour selection, we can perform edge detection and sharpening, and even come up with a way to auto-scale an unsharp mask to mostly avoid halos. Finally, we tried converting the filters into two-pass versions, which scale significantly better. It seems to work fairly well for blurring filters, but has mixed results for the sharpening version. There might be some ways to address the latter issues though. Also, if anyone responsible for implementing GIMP filters is reading this: maybe we can expand the existing SNN filter to include a few of the ones explored here? ;) Thanks for reading, hope it was interesting!`48:T480,// Assumes the filter will return a typed array function makeFilterWorker(filter, params) { const script = makeFilterWorkerScript(filter); const context = DOM.context2d(params.width, params.height, 1); const imageData = context.getImageData(0, 0, params.width, params.height); // use original image to start with imageData.data.set(params.source); context.putImageData(imageData, 0, 0); // emphasize that this is the unrendered image context.font = '40px serif'; context.fillStyle = '#00000088'; context.fillRect(0, 0, params.width, params.height); context.fillText('Rendering...', 24, 64); context.fillStyle = 'white'; context.fillText('Rendering...', 20, 60); context.strokeStyle = 'black'; context.strokeText('Rendering...', 20, 60); const worker = new Worker(script); const messaged = ({data: {rgba}}) => { imageData.data.set(rgba); context.putImageData(imageData, 0, 0); }; invalidation.then(() => worker.terminate()); worker.onmessage = messaged; worker.postMessage({params}); context.canvas.style = "image-rendering: crisp-edges; image-rendering: pixelated" return context.canvas; }49:T40a,// Assumes the filter will return a typed array function makeFilterWorkerScript(filter) { const filterScript = filter.toString(); const script = ` // Inline all the generic functions used by the filter ${colorDelta.toString()} ${filterScript.indexOf("selectNearest") !== -1 ? selectNearest.toString() : ""} ${filterScript.indexOf("selectFurthest") !== -1 ? selectFurthest.toString() : ""} ${filterScript.indexOf("selectDipole") !== -1 ? selectDipole.toString() : ""} ${filterScript.indexOf("bellcurvish2D") !== -1 ? bellcurvish2D.toString() : ""} ${filterScript.indexOf("bellcurvish1D") !== -1 ? bellcurvish1D.toString() : ""} // Inline the filter ${filterScript} onmessage = event => { const {params} = event.data; const rgba = ${filter.name}(params); postMessage({rgba}, [rgba.buffer]); close(); } `; const url = URL.createObjectURL(new Blob([script], {type: 'text/javascript'})) invalidation.then(() => URL.revokeObjectURL(url)); return url; }4a:[] 2f:{"title":"Variations on Symmetric Nearest Neighbour smoothing","cells":[{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":35,"value":"tulip = await new Promise((resolve, reject) => {\n const image = new Image();\n image.crossOrigin = \"anonymous\";\n image.onerror = reject;\n image.onload = () => {\n resolve(image);\n };\n image.src = \"https://blindedcyclops.neocities.org/ObservableData/SNN/tulip.jpg\";\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":999,"value":"kiss = await new Promise((resolve, reject) => {\n const image = new Image();\n image.crossOrigin = \"anonymous\";\n image.onerror = reject;\n image.onload = () => {\n // image is 1200 x 600 pixels,\n // if width is less than that we crop \n // rather than resize to avoid artefacts\n if (width >= 1200) resolve(image);\n const ctx = DOM.context2d(width, 600);\n ctx.drawImage(image, 0, 0);\n resolve(ctx.canvas);\n };\n image.src = \"https://blindedcyclops.neocities.org/ObservableData/SNN/kiss.jpg\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":115,"value":"md`The tulip photo was picked because it is blurry at the edges, having nice examples of gradients in real-world photos, but also has elements that are in sharp focus. The dust in the corner will be a good point of reference too. The kissing couple combines fog with grain to create a very noisy signal that human eyes can extract information from quite well, but naive image filters tend to struggle with (for example, the subtle brick road texture that fades away).\n\nFirst we should recreate the original Symmetric Nearest Neighbour filter. This is done by combining a symmetric nearest neighbour selection algorithm with a box blur. Let's start with implementing the box blur to see if that works as expected:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":93,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":96,"value":"boxblurTulip = {\n const blurParams = makeFilterParams(tulipData, 10);\n return makeFilterWorker(boxblur, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1010,"value":"boxblurKiss = {\n const blurParams = makeFilterParams(kissData, 10);\n return makeFilterWorker(boxblur, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":81,"value":"md`Next we write a generic \\`selectNearest\\` selection function. To make it easier to plug into our filters, it will be applied directly to the \\`Uint8Array\\` of the image data, and return the index of the picked color in that array:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":22,"value":"$33","pinned":true,"hidden":false,"mode":"ojs"},{"id":1021,"value":"gaussianBlurSNNKiss = {\n const blurParams = makeFilterParams(kissData, 10);\n return makeFilterWorker(gaussianBlurSNN, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":23,"value":"$34","pinned":true,"hidden":false,"mode":"ojs"},{"id":223,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":32,"value":"$36","pinned":true,"hidden":false,"mode":"ojs"},{"id":194,"value":"md`Works like a charm! Ok, onto turning that into a symmetric nearest neighbour variant:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":46,"value":"boxblurSNNTulip = {\n const blurParams = makeFilterParams(tulipData, 10);\n return makeFilterWorker(boxblurSNN, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1013,"value":"boxblurSNNKiss = {\n const blurParams = makeFilterParams(kissData, 10);\n return makeFilterWorker(boxblurSNN, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":134,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":1289,"value":"function bellcurvish2D(radius) {\n const p = new Uint32Array(radius+1);\n for (let i = 0; i <= (radius+1)/2; i++){\n for (let j = 0; j <= radius/2; j++) {\n p[i+j]++;\n }\n }\n const q = new Uint32Array(p.length*2 - 1);\n for (let i = 0; i < p.length; i++){\n for (let j = 0; j < p.length; j++) {\n q[i+j] += p[j];\n }\n }\n \n let total = 0;\n const quadratic2D = [];\n for (let i = 0; i < q.length; i++){\n const c = q[i];\n const q2 = [];\n for(let j = 0; j < q.length; j++){\n const val = c * q[j];\n total += val;\n q2.push(val);\n }\n quadratic2D.push(q2);\n }\n const norm = 1/total;\n for (let i = 0; i < q.length; i++){\n for(let j = 0; j < q.length; j++){\n quadratic2D[i][j] *= norm;\n }\n }\n return quadratic2D;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":163,"value":"md`Now let's write a Gaussian(ish) blur:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":169,"value":"function gaussianBlur({source, width, height, radius}) {\n if (!(radius > 0)) return source.slice();\n let target = new Uint8ClampedArray(source.length);\n const kernel2 = bellcurvish2D(radius);\n \n const {min, max} = Math;\n for(let y = 0; y < height; y++) {\n const line = y * width * 4;\n for (let x = 0; x < width; x++) {\n let r = 0;\n let g = 0;\n let b = 0;\n for (let dy = -radius; dy <= radius; dy++) {\n const line_dy = min(max(0, y + dy), height-1) * width;\n for (let dx = -radius; dx <= radius; dx++) {\n const idx = (min(max(0, x + dx), width-1) + line_dy) * 4;\n const weight = kernel2[radius + dx][radius + dy];\n r += source[idx] * weight;\n g += source[idx + 1] * weight;\n b += source[idx + 2] * weight;\n }\n }\n const cIdx = x*4 + line;\n target[cIdx] = r;\n target[cIdx+1] = g;\n target[cIdx+2] = b;\n target[cIdx+3] = 255;\n }\n }\n\n return target;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":87,"value":"md`Now we can plug that \\`selectNearest\\` function into our simple box blur:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":181,"value":"gaussianBlurTulip = {\n const blurParams = makeFilterParams(tulipData, 10);\n return makeFilterWorker(gaussianBlur, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1017,"value":"gaussianBlurKiss = {\n const blurParams = makeFilterParams(kissData, 10);\n return makeFilterWorker(gaussianBlur, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":208,"value":"gaussianBlurSNNTulip = {\n const blurParams = makeFilterParams(tulipData, 10);\n return makeFilterWorker(gaussianBlurSNN, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":203,"value":"$38","pinned":true,"hidden":false,"mode":"ojs"},{"id":1685,"value":"$39","pinned":false,"hidden":false,"mode":"ojs"},{"id":1687,"value":"gaussianBlurSFNTulip = {\n const blurParams = makeFilterParams(tulipData, 10);\n return makeFilterWorker(gaussianBlurSFN, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2037,"value":"gaussianBlurSFNKiss = {\n const blurParams = makeFilterParams(kissData, 10);\n return makeFilterWorker(gaussianBlurSFN, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1694,"value":"md`So our edge-enhancing blur turned into an edge-*inverting* blur. That might seem kind of useless at first beyond intentional glitch art, but think about what basically happened here: if we pick the color that is the least like the central one, then pixels near an edge will likely pick a color outside of the edge. So the pixels near an edge will \"flip\", as we see in the above images, and the closer to the edge the stronger the flip.\n\nThis means that our Symmetrical Furthest Neighbour selection method responds most strongly to edges. So perhaps by combining nearest and furthers neighbour selection we could create a form of edge detection?`","pinned":false,"hidden":false,"mode":"ojs"},{"id":918,"value":"$3a","pinned":true,"hidden":false,"mode":"ojs"},{"id":1721,"value":"$3b","pinned":true,"hidden":false,"mode":"ojs"},{"id":1742,"value":"symmetricalNeighbourEdgeDetectTulip = {\n const blurParams = makeFilterParams(tulipData, 10);\n return makeFilterWorker(symmetricalNeighbourEdgeDetection, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":252,"value":"$3c","pinned":false,"hidden":false,"mode":"ojs"},{"id":1759,"value":"symmetricalNeighbourEdgeDetectKiss = {\n const blurParams = makeFilterParams(kissData, 10);\n return makeFilterWorker(symmetricalNeighbourEdgeDetection, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2234,"value":"$3d","pinned":false,"hidden":false,"mode":"ojs"},{"id":2241,"value":"symmetricalNeighbourEdgeDetectAbsoluteTulip = {\n const blurParams = makeFilterParams(tulipData, 10);\n return makeFilterWorker(symmetricalNeighbourEdgeDetectionAbsolute, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2238,"value":"symmetricalNeighbourEdgeDetectAbsoluteKiss = {\n const blurParams = makeFilterParams(kissData, 10);\n return makeFilterWorker(symmetricalNeighbourEdgeDetectionAbsolute, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1715,"value":"md`Looks like it works! Even the subtle brickwork in the kissing couple photo is slightly highlighted.\n\nNow that we have edge-detection we can use it to sharpen images. As a point of comparison, let's first look at a more traditional method: the [unsharp mask](https://en.wikipedia.org/wiki/Unsharp_masking). It turns out that once you have a Gaussian blur, you can quite easily implement a basic unsharp mask by subtracting the inverted version of the blurred image from the original image. The linked wikipedia article explains it in more detail, but here is an implementation:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":231,"value":"$3e","pinned":true,"hidden":false,"mode":"ojs"},{"id":235,"value":"unsharpMaskTulip = {\n const blurParams = makeFilterParams(tulipData, 10);\n return makeFilterWorker(unsharpMask, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1025,"value":"unsharpMaskKiss = {\n const blurParams = makeFilterParams(kissData, 10);\n return makeFilterWorker(unsharpMask, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":302,"value":"md`There we go! One unsharp mask, with all the familiar halos and over-sharpened noise. \n\nSo how to turn our earlier symmetric neighbour edge-detection method into a sharpening filter? Well, the output of the first edge-detection function showed darker or lighter pixels near edges, with lightness depending on which \"side\" of the edge our filter was on - a bit like a positive or negative sign. We did that by subtracting the difference of the sum of the nearest and furthest neighbors from a neutral gray value. All we have to do now is replace this neutral grey value with the value of the original pixel:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1807,"value":"$3f","pinned":true,"hidden":false,"mode":"ojs"},{"id":1810,"value":"symmetricalNeighbourSharpenTulip = {\n const blurParams = makeFilterParams(tulipData, 10);\n return makeFilterWorker(symmetricalNeighbourSharpen, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1815,"value":"symmetricalNeighbourSharpenKiss = {\n const blurParams = makeFilterParams(kissData, 10);\n return makeFilterWorker(symmetricalNeighbourSharpen, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2035,"value":"$40","pinned":false,"hidden":false,"mode":"ojs"},{"id":2054,"value":"$41","pinned":true,"hidden":false,"mode":"ojs"},{"id":2056,"value":"symmetricalNeighbourSharpenNoHaloTulip = {\n const blurParams = makeFilterParams(tulipData, 10);\n return makeFilterWorker(symmetricalNeighbourSharpenNoHalo, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2065,"value":"symmetricalNeighbourSharpenNoHaloKiss = {\n const blurParams = makeFilterParams(kissData, 10);\n return makeFilterWorker(symmetricalNeighbourSharpenNoHalo, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2078,"value":"md`Well, the halos are indeed gone! But isn't this just like our earlier Gaussian SNN blur now? It does look similar, but upon closer inspection we can see that contrast in fine details has been enhanced in this version: the dust in the upper-left corner of the tulip photo was smoothed away in the GSNNB filter, while it survives here. Similarly, small highlights in the kissing couple photo, like the lights in the background and the reflections on the bike were slightly blurred out with the GSNNB, while their contrast is higher in these outputs.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":911,"value":"$42","pinned":false,"hidden":false,"mode":"ojs"},{"id":669,"value":"$43","pinned":false,"hidden":false,"mode":"ojs"},{"id":680,"value":"$44","pinned":false,"hidden":false,"mode":"ojs"},{"id":684,"value":"boxblurSNNTwoPassTulip = {\n const blurParams = makeFilterParams(tulipData, 10);\n return makeFilterWorker(boxblurSNNTwoPass, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1053,"value":"boxblurSNNTwoPassKiss = {\n const blurParams = makeFilterParams(kissData, 10);\n return makeFilterWorker(boxblurSNNTwoPass, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":690,"value":"md`That turned out... surprisingly good actually? And since it now scales a lot better we can try using a big radius and see what the results are:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":778,"value":"bigboxblurTwoPassTulipSNN = {\n const blurParams = makeFilterParams(tulipData, 50);\n return makeFilterWorker(boxblurSNNTwoPass, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1064,"value":"bigboxblurTwoPassKissSNN = {\n const blurParams = makeFilterParams(kissData, 50);\n return makeFilterWorker(boxblurSNNTwoPass, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1837,"value":"md`... I'm not sure what I was expecting: I mentioned before that when the radius is too big, we get \"inversion\" problems. There is nothing about the two-pass version that would avoid this problem, so obviously we should see that effect too. \n\nSo if this works for box blur, we should expect the Gaussian blur to work as well. We need the one-dimensional version of our bellcurve approximation, but otherwise the principle is the same:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1363,"value":"function bellcurvish1D(radius) {\n const p = new Uint32Array(radius+1);\n for (let i = 0; i <= (radius+1)/2; i++){\n for (let j = 0; j <= radius/2; j++) {\n p[i+j]++;\n }\n }\n const q = new Uint32Array(p.length*2 - 1);\n for (let i = 0; i < p.length; i++){\n for (let j = 0; j < p.length; j++) {\n q[i+j] += p[j];\n }\n }\n\n const quadratic = [];\n for (let v of q){\n quadratic.push(v*v);\n }\n return quadratic;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":701,"value":"$45","pinned":true,"hidden":false,"mode":"ojs"},{"id":704,"value":"gaussianBlurSNNTwoPassTulip = {\n const blurParams = makeFilterParams(tulipData, 10);\n return makeFilterWorker(gaussianBlurSNNTwoPass, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1056,"value":"gaussianBlurSNNTwoPassKiss = {\n const blurParams = makeFilterParams(kissData, 10);\n return makeFilterWorker(gaussianBlurSNNTwoPass, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2109,"value":"md`Also working fine. And finally, let's try adapting our no-halo sharpening:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2084,"value":"$46","pinned":true,"hidden":false,"mode":"ojs"},{"id":2090,"value":"symmetricalNeighbourSharpenNoHaloTwoPassTulip = {\n const blurParams = makeFilterParams(tulipData, 10);\n return makeFilterWorker(symmetricalNeighbourSharpenNoHaloTwoPass, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2091,"value":"\nsymmetricalNeighbourSharpenNoHaloTwoPassKiss = {\n const blurParams = makeFilterParams(kissData, 10);\n return makeFilterWorker(symmetricalNeighbourSharpenNoHaloTwoPass, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2142,"value":"md`Well, our two-pass version unsharpens, but it also re-introduces the halos and amplefies the grain noise again. Less so than a naive unsharp mask though. It may also be a matter of tweaking the strength of the filter through other means, because if I set the radius to 5 we seem to hit a sweet spot in sharpening without visible halos:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2164,"value":"symmetricalNeighbourSharpenNoHaloTwoPassTulip5 = {\n const blurParams = makeFilterParams(tulipData, 5);\n return makeFilterWorker(symmetricalNeighbourSharpenNoHaloTwoPass, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2161,"value":"\nsymmetricalNeighbourSharpenNoHaloTwoPassKissR5 = {\n const blurParams = makeFilterParams(kissData, 5);\n return makeFilterWorker(symmetricalNeighbourSharpenNoHaloTwoPass, blurParams);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2171,"value":"$47","pinned":false,"hidden":false,"mode":"ojs"},{"id":166,"value":"html`

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":38,"value":"tulipData = {\n width; // trigger resize when going to fullscreen, or landscape on mobile\n const ctx = DOM.context2d(tulip.width, tulip.height, 1);\n ctx.drawImage(tulip, 0, 0, ctx.canvas.width, ctx.canvas.height);\n const imageData = ctx.getImageData(0, 0, ctx.canvas.width, ctx.canvas.height);\n return imageData;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1005,"value":"kissData = {\n width;\n const ctx = DOM.context2d(kiss.width, kiss.height, 1);\n ctx.drawImage(kiss, 0, 0, ctx.canvas.width, ctx.canvas.height);\n const imageData = ctx.getImageData(0, 0, ctx.canvas.width, ctx.canvas.height);\n return imageData;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1028,"value":"function makeFilterParams(imageData, radius) {\n return {\n source: imageData.data,\n width: imageData.width,\n height: imageData.height,\n radius,\n };\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":544,"value":"$48","pinned":true,"hidden":false,"mode":"ojs"},{"id":1112,"value":"function identityFilter({source}){\n return source;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":568,"value":"$49","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W4a","stdlib":"1","titling":"auto"}