diff --git a/API.md b/API.md index 7542027e..006329f1 100644 --- a/API.md +++ b/API.md @@ -72,6 +72,8 @@ Create an empty Paint layer with `addLayer("paint")`. `addPaintStroke(id, stroke For RAW sources, `setWhiteBalance({ temperature, tint })` sets absolute Kelvin and DNG tint, preserving unspecified values. Temperature accepts 2000–25000 K and tint accepts -150–150, extending either range to include the file's As Shot value. `setWhiteBalance()` restores that value (the decoder's daylight fallback if camera multipliers are unavailable). Non-RAW sources and invalid values throw. Incremental temperature/tint remain separate RGB adjustments. +`setNoiseReduction({ luminance, color })` reduces the photo's noise before its adjustments, keeping the strength it omits: a RAW photo's on its Bayer mosaic before demosaicing, any other image's in its encoded color. Each strength goes from 0, none, to 100; near 50 the result comes closest to the clean photo, and higher smooths further. Luminance reduces noise in light, color in color, so color alone leaves a monochrome grain. A photo's first reduction runs once on the GPU and takes a moment while the editor keeps showing the photo as it was; later strengths, white balance, and every other edit reuse it, so they apply at once. The original comparison shows the photo as decoded. RAW files without a 2 × 2 mosaic, such as Fujifilm X-Trans, and invalid values throw. + `autoWhiteBalance()` neutralizes the photo's color cast as one edit and resolves once it is made. The GPU measures the cast of the source before any edit, as its edges' average color, so calling it again gives the same balance. A RAW photo's white balance moves from As Shot; any other photo's incremental temperature and tint are replaced. When another photo opens first, it resolves without an edit. Edits update the scene and history synchronously. Rendering may finish later, particularly RAW development. Tests should wait for visible results; `exportImage()` renders and waits for its captured scene independently of the preview. @@ -185,7 +187,7 @@ The image renders at the document dimensions and downsamples to `longEdge` with `scene` is the `getState()` scene; each image layer's `source` names an entry in `sources`, and each Remove patch's `field` an entry in `fields`: where its texels sit, or `{}` for a field not yet synthesized. Version 1 files, which kept each Remove field inside its patch, open as version 2. Opening the file with `loadScene`, `openFile`, a drop, or the file picker decodes the stored source again and restores the frame and every layer as a new document with empty history. Preview settings and history are not saved. -Opening validates every value as the matching command does, and a file that fails leaves the workspace in its error state with a message naming the first invalid field. Fields OpenLight does not know are dropped. A parameter missing from `adjustments`, `details`, `vignette`, `grain`, `fill`, or `colorMixer` takes its default, so older files still open when a group gains a parameter; a RAW image without a white balance uses its As Shot value, and a frame without `perspective` has none. `version` increases only when older files can no longer open as written, or an older OpenLight would drop what a file holds: a scene that corrects perspective is version 3, which OpenLight before perspective correction rejects, and every other scene stays version 2. A newer version is rejected. +Opening validates every value as the matching command does, and a file that fails leaves the workspace in its error state with a message naming the first invalid field. Fields OpenLight does not know are dropped. A parameter missing from `adjustments`, `details`, `vignette`, `grain`, `fill`, or `colorMixer` takes its default, so older files still open when a group gains a parameter; a RAW image without a white balance uses its As Shot value, an image without `noiseReduction` has none, and a frame without `perspective` has none. `version` increases only when older files can no longer open as written, or an older OpenLight would drop what a file holds: a scene that corrects perspective is version 3, which OpenLight before perspective correction rejects, and every other scene stays version 2. A newer version is rejected. ## Drafts @@ -218,6 +220,7 @@ editor.run({ | `set-adjustments` | adjustments, `layerId?` | Like `setAdjustments`. | | `set-tone-curve` | `points?`, `layerId?` | Like `setToneCurve`. | | `set-white-balance` | `temperature?`, `tint?` | Like `setWhiteBalance` with a change. | +| `set-noise-reduction` | `luminance?`, `color?` | Like `setNoiseReduction`. | | `set-color-mixer` | `color`, `hue?`, `saturation?`, `luminance?`, `layerId?` | Like `setColorMixer`. | | `set-details` | `clarity?`, `sharpening?`, `sharpenRadius?`, `layerId?` | Like `setDetails`. | | `set-vignette` | `intensity?`, `softness?`, `layerId?` | Like `setVignette`. | @@ -225,7 +228,7 @@ editor.run({ | `add-mask` | `mask`, `adjustments?` | Adds a mask layer with those adjustments on top of the stack as one edit. | | `delete-layer` | `layerId` | Like `deleteLayer`. | | `set-crop` | `aspectRatio?`, `straighten?` | Replaces the frame with the largest centered crop of the source at `aspectRatio`, width over height, straightened by −45 to 45 degrees, keeping the frame's perspective. Without either, it removes the crop. | -| `reset` | none | Removes every layer and returns adjustments, the tone curve, white balance, and the frame, perspective included, to how the photo opened, as one edit. | +| `reset` | none | Removes every layer and returns adjustments, the tone curve, white balance, noise reduction, and the frame, perspective included, to how the photo opened, as one edit. | | `undo`, `redo` | none | Like `undo` and `redo`. | | `set-preview` | `comparison` | Shows `"edited"`, `"original"`, or `"split"`. | diff --git a/ARCHITECTURE.md b/ARCHITECTURE.md index 03fe88e3..68d9c6c5 100644 --- a/ARCHITECTURE.md +++ b/ARCHITECTURE.md @@ -94,4 +94,5 @@ Dependencies point downward. Features do not import each other; `app/` connects - A mask without child effects renders its adjustments, curve, and mix in one pass over the image below, so no full-size texture holds its edited image; while it is inspected, its curve input reads the adjusted image, so the three render separately. A mask that only gradients shape needs no texture: the mix pass computes it, children included. A mask with a painted brush or a range combines in the render graph: each child's coverage, a brush's cached raster, a gradient drawn by a node, or a range read from the image below the mask, folds into the group's in stored order through one combine pass, and the mix pass shares its add, subtract, and intersect. The group takes the image's resolution when a range takes part, since a range follows the photo's edges, and the brushes' otherwise. The renderer gives it out, with each range's own coverage, for the overlay and thumbnails. While a color is picked, renders also keep the image below the mask group, which the picker reads through one GPU pass and 16 bytes of readback per pick. - Renders reduce the source by the largest whole factor that still gives every device pixel a texel, a proxy, in gestures or not, so the photo's own size renders once the display shows more than half a device pixel per source pixel; the main canvas writes that density into the preview state, and the provider renders again when it changes. A proxy that lacks a Remove patch's field synthesizes it apart, outside a gesture: the layers below the patch render in full for the field alone, and the proxy composes again with it, so the display never shows a full render for it. The export view renders its own full snapshot, which its preview and the saved file both encode. Every render image carries `scale`, its source pixels per texel; shaders that take document coordinates or radii apply it. - The scene's frame maps output to source in one homogeneous transform, `frameTransform` in `core/image/frame.ts`: crop, rotation, flips, and scale, then the perspective correction about the source's center. Layers, paint, masks, and patches stay in source pixels, so correcting perspective leaves them in place; `transformImage` resamples the composite once, the display draws the comparison and the mask overlay from the source through the same transform, and a frame that changes nothing adds no pass. The proxy follows the frame's strongest magnification, perspective's included. `lib/projective.ts` holds the plane projective geometry the frame, the Crop tool's coverage, and the canvas mapping share. -- TIFF and camera RAW decode through `raw-webgpu`; OpenLight adapts its resources to document ownership. A RAW source keeps its as-shot development, which a renderer shows as long as that balance is selected; another balance develops into the renderer's own texture. +- TIFF and camera RAW decode through `raw-webgpu`; OpenLight adapts its resources to document ownership. A RAW source keeps its as-shot development, which a renderer shows as long as that balance is selected and noise reduction is off, and the original comparison always shows; another balance, or noise reduction, develops into the renderer's own texture. +- Noise reduction works before the image layer's adjustments, on a RAW photo's 2 × 2 mosaic before demosaicing or on any other image's encoded color, at the layer's luminance and color strengths. `features/noise-reduction` reduces the noise on the GPU once per source, the first time a renderer or control asks; the image source keeps the reduction until it closes and shares it with every renderer, and each renderer composes it at the layer's strengths in one pass into a texture of its own. A RAW photo develops the composed samples: the package develops only its own sensor texture, so `core/image/decode/raw.ts` moves the decoded samples to a copy the first time and, right before each development, writes the samples that development takes into that texture, in the same task, so no other renderer's write falls between. Any other image's composition takes the composed image in place of the source. Changing a strength composes and develops again, as a white balance does, without reducing noise again; a newer request need not wait for a first reduction. `app/editor/renderer.ts` gives the renderer the feature's reduction, so the core imports no feature. See [the algorithm](src/features/noise-reduction/algorithm.md). diff --git a/CONTEXT.md b/CONTEXT.md index 9822c021..71fabad8 100644 --- a/CONTEXT.md +++ b/CONTEXT.md @@ -10,6 +10,7 @@ OpenLight edits documents containing image content and adjustments. | **Stroke** | One brush drag: paint or erase, size, feather, flow, and points with pressure in source pixels; a Paint layer's strokes also carry a color. Strokes are content; the renderer rasterizes them into a cached texture at half the source's resolution, coverage for a mask and premultiplied color for a Paint layer, and the latest stroke into a shared stroke buffer until the next one starts. A Paint layer's or brush mask's strokes settle into pixels every 100 strokes, which it then paints over; undo still reaches every stroke. A Healing patch stores ordered hard paint/erase strokes, rasterized over only the tiles its paint reaches, and feathers its accumulated shape separately. | | **Range** | A mask that selects pixels of the image below its mask group by lightness, a luminance range, or by color, a color range, rather than by position. A new color range has no coverage until it takes a color from a click on the photo. | | **Healing patch** | One painted, non-destructive repair in a Healing layer. Its Heal mode copies a nearby donor, found automatically or set by hand, with a bounded boundary color correction. Clone mode copies the donor without color correction. Remove mode synthesizes many local donor patches from surrounding texture and has no donor offset; it keeps where each pixel copies from as its field, so the layers below change only the copied colors, and a stroke added or erased synthesizes its whole new shape. Patches replay in order and each sees the result of the preceding patches. Each gesture keeps its recorded size and separate path; Shift adds and Alt/Option subtracts on the selected patch. Feather and opacity remain editable on its accumulated shape. | +| **Noise reduction** | The image layer's luminance and color strengths, each 0 to 100, that remove noise before its adjustments: from a RAW photo's color filter mosaic before demosaicing, from any other image's color. The photo's first reduction takes a moment; changing a strength or editing anything else afterwards applies at once, and the original comparison shows the photo as decoded. RAW photos without a 2 × 2 Bayer mosaic have none. | | **Perspective correction** | Part of the frame: Vertical and Horizontal keystone that makes lines converging in the photo, such as a building's edges, parallel, before the crop, rotation, and flips. Layers keep their source pixels and follow it; tools draw and place them as it maps them. | | **Image source** | An image used as input to the document's composition. It is content within a document. | | **LUT** | A 3D color lookup table, such as a creative grade, read from a `.cube` file. A LUT layer holds one and applies it to the image below, at the layer's opacity. | diff --git a/src/app/commands.ts b/src/app/commands.ts index a70bd0ba..d9232426 100644 --- a/src/app/commands.ts +++ b/src/app/commands.ts @@ -16,6 +16,8 @@ import { setGrain } from "@/features/grain/edits"; import { grainSchema } from "@/features/grain/model"; import { addLayer, deleteLayer } from "@/features/layers/edits"; import { maskSchema } from "@/features/layers/model"; +import { setNoiseReduction } from "@/features/noise-reduction/edits"; +import { noiseReductionSchema } from "@/features/noise-reduction/model"; import { curveSchema, defaultCurve } from "@/features/tone-curves/curve"; import { setToneCurve } from "@/features/tone-curves/edits"; import { setVignette } from "@/features/vignette/edits"; @@ -65,6 +67,7 @@ function reset(document: EditorDocument) { adjustments: { ...defaultAdjustments }, toneCurve: defaultCurve, whiteBalance: source.raw?.asShot, + noiseReduction: undefined, }, ], }); @@ -104,6 +107,14 @@ export const commands = { return { layerId: image(document).id }; }, ), + "set-noise-reduction": command( + "Reduces the photo's noise before its adjustments, a RAW photo's before demosaicing: luminance and color from 0, none, to 100, keeping omitted values. 50 removes the noise the photo measures; more smooths further. The first reduction of a photo takes a moment; later strengths apply at once. RAW files without a 2 × 2 mosaic, such as Fujifilm X-Trans, have none.", + change(noiseReductionSchema), + (document, change) => { + setNoiseReduction(document, change); + return { layerId: image(document).id }; + }, + ), "set-color-mixer": command( "Shifts one color range's hue, saturation, or luminance from -100 to 100, keeping other colors and omitted values. Neutral grays are unaffected. Without layerId it edits the first Color Mixer layer, creating one when there is none.", z.extend(mixerChange, { color: mixerColor, layerId }), diff --git a/src/app/controls.ts b/src/app/controls.ts index ced6aa91..0897d3c0 100644 --- a/src/app/controls.ts +++ b/src/app/controls.ts @@ -28,7 +28,7 @@ import type { ToneCurve, Vignette, } from "@/core/document"; -import type { WhiteBalance } from "@/core/image"; +import type { NoiseReduction, WhiteBalance } from "@/core/image"; import { decode } from "@/core/image/decode"; import type { ImageFrame, Point } from "@/core/image/frame"; import { setAdjustments, setExposure } from "@/features/adjustments/edits"; @@ -62,6 +62,7 @@ import { setMaskOperation, } from "@/features/layers/edits"; import { defaultGradient } from "@/features/layers/gradient"; +import { setNoiseReduction } from "@/features/noise-reduction/edits"; import { addPaintStroke, setPaintBlend } from "@/features/paint/edits"; import { defaultCurve } from "@/features/tone-curves/curve"; import { setToneCurve } from "@/features/tone-curves/edits"; @@ -123,6 +124,8 @@ export function createControls( setWhiteBalance: (change?: Partial) => setWhiteBalance(workspace.getDocument(), change), autoWhiteBalance: () => autoWhiteBalance(workspace.getDocument(), gpu), + setNoiseReduction: (change: Partial) => + setNoiseReduction(workspace.getDocument(), change), setToneCurve: (curve?: ToneCurve, id?: string) => setToneCurve(workspace.getDocument(), curve, id), setColorMixer(color: MixerColor, change: MixerChange, id?: string) { diff --git a/src/app/editor/adjust/index.tsx b/src/app/editor/adjust/index.tsx index 84de6ff0..3a97fa12 100644 --- a/src/app/editor/adjust/index.tsx +++ b/src/app/editor/adjust/index.tsx @@ -1,5 +1,6 @@ import { Button } from "@roprgm/ui/button"; import { Slider } from "@roprgm/ui/slider"; +import { Spinner } from "@roprgm/ui/spinner"; import { type ReactNode, useCallback, useState } from "react"; import { useStore } from "zustand"; import { effectKinds } from "@/app/editor/layers"; @@ -41,6 +42,10 @@ import { Histogram } from "@/features/histogram"; import { setLayer } from "@/features/layers/edits"; import { OverlayToggle } from "@/features/layers/overlay-toggle"; import { LutCurves } from "@/features/lut/curves"; +import { + NoiseReductionControls, + useNoiseReduction, +} from "@/features/noise-reduction/controls"; import { PaintControls } from "@/features/paint/controls"; import { defaultCurve } from "@/features/tone-curves/curve"; import { setToneCurve } from "@/features/tone-curves/edits"; @@ -148,7 +153,9 @@ function SelectedControls({ layer }: { layer: Layer }) { id={layer.id} adjustments={layer.adjustments} whiteBalance={} - /> + > + {layer.kind === "image" && } +
@@ -240,6 +247,7 @@ function AdjustmentDials({ }) { const document = useDocument(); const whiteBalance = useWhiteBalance(layer); + const noiseReduction = useNoiseReduction(layer); const header = ( ); } + // Noise reduction follows the tone, with a spinner while the first reduction runs. return ( + {noiseReduction?.reducing && ( + + )} + {action} + + } + parameters={[...adjustments(tone), ...(noiseReduction?.parameters ?? [])]} /> ); } diff --git a/src/app/editor/renderer.ts b/src/app/editor/renderer.ts index bf6ccdbc..8dd43cca 100644 --- a/src/app/editor/renderer.ts +++ b/src/app/editor/renderer.ts @@ -24,6 +24,7 @@ import { grain } from "@/features/grain/pass"; import { heal, retainHealPatches } from "@/features/heal/pass"; import { histogramGrid } from "@/features/histogram/histogram"; import { lut } from "@/features/lut/pass"; +import { reduceNoise } from "@/features/noise-reduction/reduce"; import { paint } from "@/features/paint/pass"; import { toneCurves } from "@/features/tone-curves/pass"; import { vignette } from "@/features/vignette/pass"; @@ -181,7 +182,8 @@ function composeLayers( /** * Pure layer composition shares the same graph for preview, crop, and export. - * The image layer develops the photo first; layers above it process that floating-point result. + * The image layer develops the photo first, after its noise reduction; layers above it process that + * floating-point result. * Display and export encode only after composition, so later layers can recover HDR headroom. */ export function createEditorRenderer( @@ -220,6 +222,9 @@ export function createEditorRenderer( rangeSource: children.rangeSource ?? composite.rangeSource, }; }, - options, + { + ...options, + reduceNoise: (source, signal) => reduceNoise(gpu, source, signal), + }, ); } diff --git a/src/app/scene-file.ts b/src/app/scene-file.ts index 9fa17273..6513ec20 100644 --- a/src/app/scene-file.ts +++ b/src/app/scene-file.ts @@ -32,6 +32,7 @@ import { maskSchema, } from "@/features/layers/model"; import { lutSchema } from "@/features/lut/model"; +import { noiseReductionSchema } from "@/features/noise-reduction/model"; import { paintShape } from "@/features/paint/model"; import { curveSchema } from "@/features/tone-curves/curve"; import { defaultVignette, vignetteSchema } from "@/features/vignette/model"; @@ -225,6 +226,7 @@ const imageLayer = z.object({ whiteBalance: z.optional( z.object({ temperature: z.number(), tint: z.number() }), ), + noiseReduction: z.optional(noiseReductionSchema), adjustments, toneCurve: curveSchema, children: empty("The image layer cannot contain layers"), diff --git a/src/core/document/scene.ts b/src/core/document/scene.ts index f3eec98f..9fdeaf59 100644 --- a/src/core/document/scene.ts +++ b/src/core/document/scene.ts @@ -1,4 +1,4 @@ -import type { WhiteBalance } from "@/core/image"; +import type { NoiseReduction, WhiteBalance } from "@/core/image"; import type { ImageFrame, Point } from "@/core/image/frame"; export type Adjustments = { @@ -72,6 +72,8 @@ export type ImageLayer = { readonly name: string; readonly source: string; readonly whiteBalance?: Readonly; + /** Noise reduction on the photo before any adjustment, a RAW one's before demosaicing; none without it. */ + readonly noiseReduction?: Readonly; readonly adjustments: Readonly; readonly toneCurve: ToneCurve; readonly children: readonly ProcessingLayer[]; diff --git a/src/core/image/decode/raw-blend.wgsl b/src/core/image/decode/raw-blend.wgsl new file mode 100644 index 00000000..8e45f929 --- /dev/null +++ b/src/core/image/decode/raw-blend.wgsl @@ -0,0 +1,32 @@ +// Writes sensor samples for the package to develop from replacement samples, each 2 × 2 cell's four by +// position above black in one texel, dithered so smooth replacements don't band when rounded back to +// whole sensor units. +@group(0) @binding(0) var decoded: texture_2d; +@group(0) @binding(1) var replacement: texture_2d; +// Each 2 × 2 position's black level. +@group(0) @binding(2) var black: vec4f; + +// One triangle covering the whole destination. +@vertex fn vs_main(@builtin(vertex_index) index: u32) -> @builtin(position) vec4f { + let points = array(vec2f(-1, -1), vec2f(3, -1), vec2f(-1, 3)); + return vec4f(points[index], 0, 1); +} + +// Uniform in [-0.5, 0.5), from a hash of the sample's position (PCG, Jarzynski and Olano, 2020). +fn dither(p: vec2u) -> f32 { + var state = p.x * 747796405u + p.y * 2891336453u + 2891336453u; + state = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u; + return f32((state >> 22u) ^ state) / 4294967296.0 - 0.5; +} + +@fragment fn fs_main(@builtin(position) position: vec4f) -> @location(0) vec4u { + let p = vec2u(position.xy); + let cell = p / 2u; + // A last odd row or column has no cell and keeps its samples. + if (any(cell >= textureDimensions(replacement))) { + return textureLoad(decoded, p, 0); + } + let site = (p.y & 1u) * 2u + (p.x & 1u); + let value = textureLoad(replacement, cell, 0)[site] + black[site] + dither(p); + return vec4u(u32(clamp(round(value), 0.0, 65535.0))); +} diff --git a/src/core/image/decode/raw.ts b/src/core/image/decode/raw.ts index 6348be43..70a21372 100644 --- a/src/core/image/decode/raw.ts +++ b/src/core/image/decode/raw.ts @@ -1,11 +1,135 @@ -import { createRawDecoder } from "raw-webgpu"; -import { type Gpu, type Target, target } from "vgpu"; -import { createImageSource, type RawDevelopment } from "@/core/image"; +import { createRawDecoder, type RawSource } from "raw-webgpu"; +import { type Gpu, type Target, type Texture, target } from "vgpu"; +import { + createImageSource, + type Mosaic, + type RawDevelopment, +} from "@/core/image"; import { weakMemo } from "@/lib/weak-memo"; +import blendShader from "./raw-blend.wgsl"; /** One decoder per GPU, released with its device. */ const decoder = weakMemo((gpu: Gpu) => createRawDecoder(gpu.gpu)); +const blendPipeline = weakMemo((device: GPUDevice) => { + const module = device.createShaderModule({ code: blendShader.wgsl }); + return device.createRenderPipeline({ + layout: "auto", + vertex: { module, entryPoint: "vs_main" }, + fragment: { + module, + entryPoint: "fs_main", + targets: [{ format: "r16uint" }], + }, + }); +}); + +/** + * The samples the package develops, which live in its own `sensor` texture. Passes take turns: each + * writes the samples it develops just before developing, in the same task, so no other pass's write + * falls between. The decoded samples move to a copy the first time others replace them. + */ +function createSampleWriter( + gpu: Gpu, + sensor: GPUTexture, + blacks: readonly number[], +) { + const device = gpu.gpu; + const black = gpu.device.createBuffer({ + size: 16, + usage: ["uniform", "copy_dst"], + }); + black.write(new Float32Array(blacks)); + let decoded: Texture | undefined; + return { + /** The decoded samples, wherever they are. */ + decoded: () => decoded?.gpu ?? sensor, + /** Writes `samples` in place of the decoded ones, or the decoded ones back. */ + write(samples?: Target) { + if (!samples && !decoded) { + return; + } + const encoder = device.createCommandEncoder(); + const size: [number, number] = [sensor.width, sensor.height]; + if (!decoded) { + decoded = gpu.device.createTexture({ + size, + format: sensor.format, + usage: ["texture_binding", "copy_src", "copy_dst"], + }); + encoder.copyTextureToTexture( + { texture: sensor }, + { texture: decoded.gpu }, + size, + ); + } + if (samples) { + const pipeline = blendPipeline(device); + const pass = encoder.beginRenderPass({ + colorAttachments: [ + { view: sensor.createView(), loadOp: "clear", storeOp: "store" }, + ], + }); + pass.setPipeline(pipeline); + pass.setBindGroup( + 0, + device.createBindGroup({ + layout: pipeline.getBindGroupLayout(0), + entries: [ + { binding: 0, resource: decoded.gpu.createView() }, + { binding: 1, resource: samples.color.gpu.createView() }, + { binding: 2, resource: { buffer: black.gpu } }, + ], + }), + ); + pass.draw(3); + pass.end(); + } else { + encoder.copyTextureToTexture( + { texture: decoded.gpu }, + { texture: sensor }, + size, + ); + } + device.queue.submit([encoder.finish()]); + }, + dispose() { + black.dispose(); + decoded?.dispose(); + }, + }; +} + +/** + * The mosaic of a source the GPU demosaics from a 2 × 2 color filter of integer samples, and what + * writes the samples its passes develop. + */ +function mosaicOf(gpu: Gpu, source: RawSource) { + const { cfa, cfaSize = 2, demosaic, black, white } = source.metadata; + if ( + !cfa || + cfaSize !== 2 || + demosaic !== "gpu" || + source.texture.format !== "r16uint" + ) { + return undefined; + } + const writer = createSampleWriter( + gpu, + source.texture, + cfa.map((color) => black[color]), + ); + const mosaic: Mosaic = { + get samples() { + return writer.decoded(); + }, + pattern: cfa, + black, + white, + }; + return { mosaic, write: writer.write, dispose: writer.dispose }; +} + /** Adapts package-owned sensor sources to OpenLight's document and renderer lifetimes. */ export async function decodeRaw(gpu: Gpu, file: File) { const source = await decoder(gpu).load(file); @@ -23,26 +147,31 @@ export async function decodeRaw(gpu: Gpu, file: File) { } const { asShot } = source; const developed = original; + const mosaic = mosaicOf(gpu, source); const raw: RawDevelopment = { asShot, + mosaic: mosaic?.mosaic, createPass() { - // The as-shot development is the image the source keeps; another balance gets its own. + // The as-shot development of the decoded samples is the image the source keeps; another + // balance, or other samples, get their own. let output: Target | undefined; const pass = source.createDevelopPass(); let calibration = source.calibration; + let replaced: Target | undefined; let dirty = true; return { - async prepare(balance) { + async prepare(balance, samples) { const shot = balance.temperature === asShot.temperature && balance.tint === asShot.tint; calibration = shot ? source.calibration : await source.calibrate(balance); + replaced = samples; dirty = true; }, render() { - if (calibration === source.calibration) { + if (calibration === source.calibration && !replaced) { output?.color.dispose(); output = undefined; return developed; @@ -52,6 +181,7 @@ export async function decodeRaw(gpu: Gpu, file: File) { format: "rgba16float", }); if (dirty) { + mosaic?.write(replaced); pass.render({ destination: output.color.gpu, calibration }); dirty = false; } @@ -63,7 +193,10 @@ export async function decodeRaw(gpu: Gpu, file: File) { }, }; }, - dispose: source.dispose, + dispose() { + mosaic?.dispose(); + source.dispose(); + }, }; return createImageSource(original, { raw }); } catch (error) { diff --git a/src/core/image/index.ts b/src/core/image/index.ts index 02026edb..59a8987a 100644 --- a/src/core/image/index.ts +++ b/src/core/image/index.ts @@ -1,10 +1,43 @@ import type { Target } from "vgpu"; export type WhiteBalance = { temperature: number; tint: number }; +/** How much noise reduction removes from light and from color, each 0 to 100; 0 keeps it. */ +export type NoiseReduction = { luminance: number; color: number }; +export const noNoiseReduction: NoiseReduction = { luminance: 0, color: 0 }; +/** A source's noise, reduced once, which each renderer composes at any strengths at once. */ +export type Reduction = { + /** A composer with its own output, for one renderer. */ + compose(): { + /** + * The source reduced by `strengths`: for a RAW photo, samples to develop in place of the decoded + * ones; for any other image, the image itself. + */ + render(strengths: NoiseReduction): Target; + dispose(): void; + }; + dispose(): void; +}; +/** A RAW photo's 2 × 2 color filter mosaic before the GPU demosaics it. */ +export type Mosaic = { + /** One `r16uint` sample per sensor pixel, unrotated, as decoded. */ + readonly samples: GPUTexture; + /** Each 2 × 2 cell position's color, row-major: 0 red, 1 or 3 green, 2 blue. */ + pattern: readonly number[]; + /** Black level per color, in sensor units. */ + black: readonly number[]; + /** The level where samples clip, in sensor units. */ + white: number; +}; export type RawDevelopment = { asShot: WhiteBalance; + /** The mosaic, when the GPU demosaics a 2 × 2 one, for processing before demosaicing. */ + mosaic?: Mosaic; createPass(): { - prepare(balance: WhiteBalance): Promise; + /** + * Develops with `balance`, from the decoded samples or from `samples`: a half-size + * `rgba16float` image holding each 2 × 2 cell's four by position, above their black level. + */ + prepare(balance: WhiteBalance, samples?: Target): Promise; render(): Target; dispose(): void; }; @@ -39,6 +72,8 @@ export function createImageSource( { raw, primaries = "rec2020" }: SourceOptions = {}, ) { let references = 1; + const closed = new AbortController(); + let reduction: Promise | undefined; function release() { let active = true; return () => { @@ -47,6 +82,11 @@ export function createImageSource( } active = false; if (--references === 0) { + closed.abort(Error("Image source is closed.")); + void reduction?.then( + (reduced) => reduced.dispose(), + () => {}, + ); image.color.dispose(); raw?.dispose(); } @@ -56,6 +96,17 @@ export function createImageSource( image, raw, primaries, + /** + * The source's noise reduction, which `reduce` makes the first time a renderer or control asks, + * kept until the source closes; a failure lets the next ask try again. + */ + reduced(reduce: (signal: AbortSignal) => Promise) { + reduction ??= reduce(closed.signal).catch((error) => { + reduction = undefined; + throw error; + }); + return reduction; + }, retain() { if (!references) { throw Error("Image source is closed."); @@ -68,3 +119,8 @@ export function createImageSource( } export type ImageSource = ReturnType; + +/** Whether a source takes noise reduction: a RAW photo through its 2 × 2 mosaic, any other image itself. */ +export function reducible(source: ImageSource) { + return !source.raw || Boolean(source.raw.mosaic); +} diff --git a/src/core/renderer/index.ts b/src/core/renderer/index.ts index db3722ff..df657223 100644 --- a/src/core/renderer/index.ts +++ b/src/core/renderer/index.ts @@ -11,7 +11,14 @@ import { type Scene, walkLayers, } from "@/core/document"; -import type { ImageSource, WhiteBalance } from "@/core/image"; +import { + type ImageSource, + type NoiseReduction, + noNoiseReduction, + type Reduction, + reducible, + type WhiteBalance, +} from "@/core/image"; import type { Point } from "@/core/image/frame"; import { createFieldStore, type FieldLattice } from "./fields"; import { createRenderGraph } from "./graph"; @@ -129,6 +136,10 @@ function reductionFor(density: number) { /** How many reductions keep their proxy and intermediates, each a quarter of the photo at most. */ const keptReductions = 4; +function sameReduction(a: NoiseReduction, b: NoiseReduction) { + return a.luminance === b.luminance && a.color === b.color; +} + function sameBalance(a: WhiteBalance | undefined, b: WhiteBalance | undefined) { return a?.temperature === b?.temperature && a?.tint === b?.tint; } @@ -142,6 +153,11 @@ export type RendererOptions = { field?: (id: string) => RemoveField | undefined; /** Receives each Remove field the renderer synthesizes, read back; without it nothing is read back. */ saveField?: (id: string, field: RemoveField) => void; + /** Reduces a source's noise once, for the image layer's noise reduction. */ + reduceNoise?: ( + source: ImageSource, + signal: AbortSignal, + ) => Promise; }; /** @@ -159,6 +175,9 @@ export function createRenderer( }, field: savedField = () => undefined, saveField, + reduceNoise = () => { + throw Error("This renderer reduces no noise."); + }, }: RendererOptions = {}, ) { const source = resource.image; @@ -181,7 +200,15 @@ export function createRenderer( /** Mask coverage the graph rendered for the overlay and thumbnails, by layer ID. */ let shown = new Map(); let balance = resource.raw?.asShot; - /** Counts developments, so the proxy follows white-balance changes. */ + /** The noise reduction the development applies. */ + let applied = noNoiseReduction; + /** Renders the reduction at the applied strengths, for the reduction it composes. */ + let composer: + | (ReturnType & { of: Reduction }) + | undefined; + /** An image's reduced self, which compositions take in place of the source. */ + let reduced: Target | undefined; + /** Counts developments, so the proxy follows white-balance and noise-reduction changes. */ let version = 0; let last: RenderRequest | undefined; /** Whether the last render, a gesture's proxy, showed a Remove patch unfilled. */ @@ -195,6 +222,8 @@ export function createRenderer( let reductions: number[] = []; let next: RenderRequest | undefined; let pending: Promise | undefined; + /** Ends a wait for noise reduction once a newer request or closing makes it moot. */ + let interrupt: (() => void) | undefined; /** A raster being read back, a painting to settle or a Remove field to save; nothing renders meanwhile. */ let settling: Promise | undefined; /** Readbacks of synthesized fields in flight, by ID, which finish even once the renderer closes. */ @@ -270,7 +299,7 @@ export function createRenderer( if (unchanged && !(lacked && !request.interactive)) { return; } - const developed = raw?.render() ?? source; + const developed = raw?.render() ?? reduced ?? source; // Every brush and paint layer updates once, bypassed or not, so hidden layers keep their rasters. for (const { layer } of walkLayers(scene.layers)) { if (layer.kind === "mask" && layer.mask.kind === "brush") { @@ -414,9 +443,33 @@ export function createRenderer( } return stale; } + /** The noise reduction a scene asks of the image layer, none where the source takes none. */ + function wantedReduction(scene: Scene) { + return reducible(resource) + ? (scene.layers[0].noiseReduction ?? noNoiseReduction) + : noNoiseReduction; + } + /** The source reduced at `strengths`, through this renderer's composer for `reduction`. */ + function composeReduction(reduction: Reduction, strengths: NoiseReduction) { + if (composer?.of !== reduction) { + composer?.dispose(); + composer = { ...reduction.compose(), of: reduction }; + } + return composer.render(strengths); + } + /** Resolves with `work`, or with nothing once a newer request or closing interrupts the wait. */ + function interruptible(work: Promise) { + const interrupted = new Promise((resolve) => { + interrupt = () => resolve(undefined); + }); + return Promise.race([work, interrupted]).finally(() => { + interrupt = undefined; + }); + } /** * Renders the latest request once what it waits for is ready: a settle, a RAW development, whose - * calibration alone crosses the worker, settled paint, or Remove fields to load. + * calibration alone crosses the worker and whose noise reduction runs once per source, which a + * newer request need not wait for, settled paint, or Remove fields to load. */ async function develop() { while (next && !disposed) { @@ -429,12 +482,40 @@ export function createRenderer( if (disposed) return; if (next) continue; const selected = scene.layers[0].whiteBalance ?? resource.raw?.asShot; - if (raw && selected && !sameBalance(balance, selected)) { - await raw.prepare(selected); - if (disposed) { - return; + const wanted = wantedReduction(scene); + let failure: { error: unknown } | undefined; + if ( + (raw && !sameBalance(balance, selected)) || + !sameReduction(wanted, applied) + ) { + let output: Target | undefined; + if (wanted.luminance || wanted.color) { + try { + const reduction = await interruptible( + resource.reduced((signal) => reduceNoise(resource, signal)), + ); + output = reduction && composeReduction(reduction, wanted); + } catch (error) { + // The photo develops as decoded rather than not at all; the error shows once, and the + // reduction runs again when its strengths change. + failure = { error }; + } + } else { + composer?.dispose(); + composer = undefined; + } + if (disposed) return; + if (next) continue; + if (raw && selected) { + await raw.prepare(selected, output); + if (disposed) { + return; + } + balance = selected; + } else { + reduced = output; } - balance = selected; + applied = wanted; version++; } for (const { id, raster, rasters } of stalePaint(scene)) { @@ -452,6 +533,9 @@ export function createRenderer( if (!next) { render(request); } + if (failure) { + throw failure.error; + } } } /** @@ -477,11 +561,13 @@ export function createRenderer( factor: reductionFor(density), interactive, }; - // Renders wait, in order, for a settle, a RAW development, or settled paint and fields to load. + // Renders wait, in order, for a settle, a RAW development or noise reduction, or settled paint and + // fields to load. if ( !raw && !pending && !settling && + sameReduction(wantedReduction(scene), applied) && !stalePaint(scene).length && !staleFields(scene).length ) { @@ -489,6 +575,7 @@ export function createRenderer( return; } next = request; + interrupt?.(); pending ??= develop() .catch((error) => { if (!next) { @@ -529,6 +616,7 @@ export function createRenderer( strokes.dispose(); proxy.dispose(); raw?.dispose(); + composer?.dispose(); release(); } @@ -632,6 +720,7 @@ export function createRenderer( return; } disposed = true; + interrupt?.(); listeners.clear(); // Banded transfers keep their targets until their last GPU operation finishes. const finishing = pending ?? settling; diff --git a/src/features/adjustments/controls.tsx b/src/features/adjustments/controls.tsx index 07ba1432..d1bb99e1 100644 --- a/src/features/adjustments/controls.tsx +++ b/src/features/adjustments/controls.tsx @@ -66,21 +66,24 @@ function AdjustmentSliders({ ); } -/** Tone, then the white balance's controls, then color. */ +/** Tone, then the white balance's controls, then color, then any controls that follow. */ export function AdjustmentControls({ id, adjustments, whiteBalance, + children, }: { id: string; adjustments: Readonly; whiteBalance: ReactNode; + children?: ReactNode; }) { return (
{whiteBalance} + {children}
); } diff --git a/src/features/noise-reduction/algorithm.md b/src/features/noise-reduction/algorithm.md new file mode 100644 index 00000000..80abc07e --- /dev/null +++ b/src/features/noise-reduction/algorithm.md @@ -0,0 +1,28 @@ +# Noise reduction + +`reduceNoise(gpu, source, signal)` reduces a photo's noise once, before its adjustments: a RAW photo on its Bayer mosaic before the GPU demosaics it, any other image in its encoded color. It runs without React, model weights, or workers, entirely in WGSL but for small fits in TypeScript. The source keeps the result until it closes, and each renderer composes it at the image layer's luminance and color strengths in one pass, so moving either slider, or editing any layer above, never reduces noise again. + +1. **Noise model, RAW.** Every 16 × 16 block of sensor pixels gives each cell position its mean and the variance of its diagonal Haar detail, which flat blocks owe to noise alone. Over 40 equal-count bins along the signal, the quietest tenth of each bin gives a point, and a weighted line through them is the position's Poisson-Gaussian model, variance = gain · signal + floor. Texture only adds variance, so the fit runs again without points more than 30 % above the line, which follows the noise where a detailed photo has few flat blocks. Positions fit apart: some cameras scale a channel digitally, which scales its noise too. +2. **Spectrum.** For a RAW photo, the generalized Anscombe transform, measured from its value at black so half floats stay precise, makes each position's noise unit Gaussian, and an orthonormal transform turns red, the two greens, and blue into light, green against magenta, red against blue, and the greens' difference. Samples below the transform's lowest value, such as dead pixels, are clamped there. Any other image is filtered in the sRGB transfer's encoding, where noise varies far less with light than in linear, as light and two color differences, each at half size and apart: light as each 2 × 2 cell's Haar transform, its mean and three differences, and color as each cell's sum. Every spectrum then takes 8 bytes per 2 × 2 pixels, so the largest texture a device takes fits one storage binding. +3. **Pyramid.** Averaging 2 × 2 texels builds up to six levels, down to 32 texels. Each level measures each component's noise from its flattest blocks. For a RAW photo, the finest level's noise is the model's, which a lower measure corrects. For an image, never below 8-bit rounding, light takes the least of its four measures that find any noise, since white noise holds the same deviation in each and texture raises a coarse cell's mean most, and each color difference its own. A coarser level takes its own measure only up to 15 % over what white noise would leave there: lossy compression correlates noise and raises it at coarse levels, while texture, which a coarse level's flattest blocks still hold, raises a measure far more. +4. **Filter.** From the coarsest level up, sliding 8 × 8 DCT denoising runs in two steps: hard thresholding at 3 deviations, 3.5 for the greens' difference, then Wiener shrinkage by the first estimate, with color shrinking twice and the greens' difference eight times as hard as light. Patches start on every other texel, mirrored past the image's edges so edges and corners get as many as the rest; each output texel gathers the 16 patches over it in workgroup memory, weighted by how little of each survived and by a Kaiser window. That is a quarter of the patches of every texel's, within 0.1 dB of them on test photos. A patch's DC stays whole where a component holds a level, so every level stays unbiased, but not for the greens' difference: both greens see the same light, so its mean is noise or imbalance, which demosaicing would draw as a maze. After each step, the upsampled difference between the coarser estimate and the step's own 2 × 2 averages, blurred at the coarse size, replaces the step's low frequencies: broad noise the 8 × 8 patches cannot tell from the image goes, the coarser level's finest frequencies, where its patches ring at edges, stay out, and where the two levels agree nothing changes. +5. **Reduction.** The filter runs once, at the measured noise. For a RAW photo it keeps, through the closed-form inverse that is unbiased for Poisson-Gaussian noise (Mäkitalo and Foi), its change to the decoded samples; for an image, its light at full size and its color at half, which a reduced image's color needs no more of. +6. **Composing.** An amount from 0 to 100 sets how much of each component's difference between the photo and its reduction goes, in noise deviations: none at 0, all at 100. Between, a component loses nearly all of a difference whose 3 × 3 neighborhood varies by less than a limit, and about the limit's worth of one that varies more: noise, which varies by one deviation, goes evenly, while detail, which varies more, mostly stays. Weighing neighborhoods rather than single values keeps the few noise values that stand out from staying as dots. At an amount a, the limit is a / (100 − a) times its value at 50: one deviation for light, three for color, which holds little detail, and 0.7 for the greens' difference, which the filter shrinks hardest. Luminance sets light, color the two color differences, and the stronger of the two the greens' difference. For a RAW photo, differences count along the transform's slope where the reduction puts each sample, as the filter saw them. + +Patches are filtered alone. Grouping similar ones, as BM3D does, ranks candidates by a distance that noise dominates in dark areas; along edges it groups aligned patches whose average keeps a profile of noise, which showed as streaks. Hot pixels are left alone: telling one from a fine highlight by its neighbors removed far more detail from sharp photos than it saved. + +On synthetic RAW mosaics made from a sharp photo, at three noise levels, an amount near 50 comes closest to the clean photo, 2 to 3 dB nearer than the noisy one, within 0.3 dB of the best a single strength reaches; 100, the whole reduction, smooths further, which very noisy photos may want. An image with added sRGB noise gains 5 to 6 dB near 50. + +## Cost + +Reducing runs the pyramid once. At its peak, the finest level's noisy spectrum, its two step estimates, the coarser estimate, and the blur between them take about 8 bytes per pixel; coarser levels free as soon as they are filtered, and an image's color waits through light's at 1 byte per pixel. Afterwards a RAW photo keeps 2 bytes per pixel for its change and 2 for its decoded samples, and each renderer 2 for its composed samples, besides the development it already makes for a white balance. An image keeps 3 bytes per pixel for its light and color, and each renderer an image in the source's own format: 4 bytes per pixel for 8-bit sources, dithered so smooth reductions don't band. + +## Limits + +Detail with less signal than noise is gone for any filter: in the darkest areas of a high-ISO photo, strong amounts leave smooth rather than textured surfaces. Very fine repeated texture near the noise level softens. A detailed photo with little noise can still measure more noise than it has, which a lower amount corrects. The RAW model assumes Poisson-Gaussian noise; row or column banding, which some sensors add, passes through. An image's noise has no model, and JPEG block edges are not noise. Filtering an image's light through 2 × 2 Haar components, which keeps it at 2 bytes per pixel, comes about 0.6 dB short of filtering it at full size. X-Trans and linear DNG files have no noise reduction. + +## Verification + +`tests/noise-reduction.test.ts` checks the noise fit against blocks of known noise, a third of them textured. `tests/noise-reduction.e2e.ts` opens a RAW fixture with seeded Poisson-Gaussian noise and a noisy PNG, checks that their noise falls well below the original while their mean holds, that color reduces apart from light and keeps the image's color, and that 0 gives back the decoded pixels. + +References: [Yu and Sapiro, DCT image denoising](https://www.ipol.im/pub/art/2011/ys-dct/), [Dabov et al., BM3D](https://webpages.tuni.fi/foi/GCF-BM3D/), [Pierazzo, Morel, and Facciolo, multiscale DCT denoising](https://www.ipol.im/pub/art/2017/201/), [Foi et al., single-image Poisson-Gaussian noise fitting](https://webpages.tuni.fi/foi/papers/Foi-PoissonianGaussianClippedRaw-2007-IEEE_TIP.pdf), [Mäkitalo and Foi, the unbiased inverse](https://webpages.tuni.fi/foi/papers/OptGenAnscombeInverse-doublecolumn-preprint.pdf). diff --git a/src/features/noise-reduction/compose-image.wgsl b/src/features/noise-reduction/compose-image.wgsl new file mode 100644 index 00000000..3d417fb1 --- /dev/null +++ b/src/features/noise-reduction/compose-image.wgsl @@ -0,0 +1,60 @@ +// An image at chosen strengths: its encoded light and color differences less part of their difference +// from the reduction's, in noise deviations, which the limits of light and color and the difference +// over the 3 × 3 pixels around set; then decoded, and dithered when the output keeps 8 bits, so +// smooth reductions don't band. +import { Limits, removed } from "./limit.wgsl"; +import { decode, encode, fromOpponent, toOpponent } from "./opponent.wgsl"; + +struct Params { + // The noise deviation of each pixel's light and two color differences. + noise: vec4f, + // The output's step in the encoding, or 0 when it keeps floats. + quantum: f32, +} + +// Of light and the two color differences. +@group(0) @binding(0) var limits: Limits; +@group(0) @binding(1) var params: Params; +@group(0) @binding(2) var source: texture_2d; +@group(0) @binding(3) var linearSampler: sampler; +@group(0) @binding(4) var light: texture_2d; +// Color at half size, sampled between its texels. +@group(0) @binding(5) var color: texture_2d; + +// Uniform in [-0.5, 0.5), from a hash of the pixel's position (PCG, Jarzynski and Olano, 2020). +fn dither(p: vec2u) -> f32 { + var state = p.x * 747796405u + p.y * 2891336453u + 2891336453u; + state = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u; + return f32((state >> 22u) ^ state) / 4294967296.0 - 0.5; +} + +struct Pixel { + value: vec3f, + // The difference from the reduction, in noise deviations. + difference: vec3f, +} + +fn pixelAt(at: vec2i) -> Pixel { + let p = clamp(at, vec2i(0), vec2i(textureDimensions(source)) - 1); + let value = toOpponent(encode(textureLoad(source, p, 0).rgb)); + // Each half-size texel averages 2 × 2 pixels, whatever the image's parity. + let uv = (vec2f(p) + 0.5) * 0.5 / vec2f(textureDimensions(color)); + let reduced = vec3f(textureLoad(light, p, 0).x, textureSampleLevel(color, linearSampler, uv, 0.0).xy); + return Pixel(value, (value - reduced) / params.noise.xyz); +} + +@fragment fn fs_main(@builtin(position) position: vec4f) -> @location(0) vec4f { + let center = vec2i(position.xy); + let pixel = pixelAt(center); + var energy = pixel.difference * pixel.difference; + for (var i = 0; i < 9; i++) { + let offset = vec2i(i % 3, i / 3) - 1; + if (any(offset != vec2i(0))) { + let difference = pixelAt(center + offset).difference; + energy += difference * difference; + } + } + let part = removed(vec4f(pixel.difference, 0.0), vec4f(energy / 9.0, 0.0), limits).xyz; + let result = fromOpponent(pixel.value - params.noise.xyz * part); + return vec4f(decode(result + dither(vec2u(center)) * params.quantum), textureLoad(source, center, 0).a); +} diff --git a/src/features/noise-reduction/compose.wgsl b/src/features/noise-reduction/compose.wgsl new file mode 100644 index 00000000..71d1e081 --- /dev/null +++ b/src/features/noise-reduction/compose.wgsl @@ -0,0 +1,39 @@ +// Replacement samples for a mosaic at chosen strengths: the decoded samples less part of their +// difference from the reduced ones, in noise deviations per spectrum component, which the +// component's limits and the difference over the 3 × 3 cells around set. Shares of one give the +// reduction, and of zero the decoded samples, exactly. +import { Limits, removed } from "./limit.wgsl"; +import { fromSpectrum, Noise, toSpectrum } from "./stabilize.wgsl"; + +@group(0) @binding(0) var noise: Noise; +// Of light, the two color differences, and the greens' difference. +@group(0) @binding(1) var limits: Limits; +@group(0) @binding(2) var samples: texture_2d; +// The reduction's change to each 2 × 2 cell's samples, by position. +@group(0) @binding(3) var reduction: texture_2d; + +@fragment fn fs_main(@builtin(position) position: vec4f) -> @location(0) vec4f { + let cell = vec2i(position.xy); + let p = cell * 2; + let decoded = vec4f( + f32(textureLoad(samples, p, 0).r), + f32(textureLoad(samples, p + vec2i(1, 0), 0).r), + f32(textureLoad(samples, p + vec2i(0, 1), 0).r), + f32(textureLoad(samples, p + vec2i(1, 1), 0).r), + ) - noise.black; + let change = textureLoad(reduction, cell, 0); + // The transform's slope where the reduction puts each sample, so differences count as the filter + // saw them. Neighbors share it: light varies little between them except at edges, which the + // difference already marks. + let a = noise.gain; + let slope = sqrt(a * max(decoded + change, vec4f(0.0)) + 0.375 * a * a + noise.floor); + let last = vec2i(textureDimensions(reduction)) - 1; + var energy = vec4f(0.0); + for (var i = 0; i < 9; i++) { + let neighbor = clamp(cell + vec2i(i % 3, i / 3) - 1, vec2i(0), last); + let difference = toSpectrum(textureLoad(reduction, neighbor, 0) / slope, noise); + energy += difference * difference; + } + let part = removed(toSpectrum(-change / slope, noise), energy / 9.0, limits); + return decoded - slope * fromSpectrum(part, noise); +} diff --git a/src/features/noise-reduction/controls.tsx b/src/features/noise-reduction/controls.tsx new file mode 100644 index 00000000..3b70d6b3 --- /dev/null +++ b/src/features/noise-reduction/controls.tsx @@ -0,0 +1,93 @@ +import { Slider } from "@roprgm/ui/slider"; +import { Spinner } from "@roprgm/ui/spinner"; +import { useEffect, useState } from "react"; +import { useGpu } from "vgpu-react"; +import type { Parameter } from "@/components/editor/parameter"; +import { useDocument } from "@/components/editor/session"; +import type { ImageLayer, Layer } from "@/core/document"; +import { type ImageSource, noNoiseReduction, reducible } from "@/core/image"; +import { setNoiseReduction } from "./edits"; +import { reduceNoise } from "./reduce"; + +/** + * Whether the source's first reduction is running while some is asked for: it takes a moment, and + * the renderer shows the photo as it was until it is ready. + */ +function useReducing(source: ImageSource | undefined, wanted: boolean) { + const gpu = useGpu(); + const [settled, setSettled] = useState(); + useEffect(() => { + if (!source || !wanted) { + return; + } + let active = true; + const settle = () => active && setSettled(source); + source + .reduced((signal) => reduceNoise(gpu, source, signal)) + .then(settle, settle); + return () => { + active = false; + }; + }, [gpu, source, wanted]); + return wanted && settled !== source; +} + +const controls = [ + ["luminance", "Luminance noise"], + ["color", "Color noise"], +] as const; + +/** + * The image layer's noise reduction, light and color, and whether its first reduction is running; + * nothing for other layers, or a RAW photo without a 2 × 2 mosaic. + */ +export function useNoiseReduction( + layer: Layer, +): { parameters: Parameter[]; reducing: boolean } | undefined { + const document = useDocument(); + const image = layer.kind === "image" ? layer : undefined; + const resource = image && document.resources.get(image.source); + const source = resource && reducible(resource) ? resource : undefined; + const strengths = image?.noiseReduction ?? noNoiseReduction; + const reducing = useReducing( + source, + strengths.luminance > 0 || strengths.color > 0, + ); + if (!source) { + return undefined; + } + return { + parameters: controls.map(([name, label]) => ({ + id: name, + label, + value: strengths[name], + min: 0, + max: 100, + defaultValue: 0, + onChange: (value) => setNoiseReduction(document, { [name]: value }), + })), + reducing, + }; +} + +/** The noise reduction sliders, with a spinner beside the first's label while the first reduction runs. */ +export function NoiseReductionControls({ layer }: { layer: ImageLayer }) { + const noiseReduction = useNoiseReduction(layer); + if (!noiseReduction) { + return null; + } + return noiseReduction.parameters.map(({ id, ...parameter }, i) => ( +
+ + {i === 0 && noiseReduction.reducing && ( + // The slider's label row again, its text hidden, so the spinner follows the label. +
+ + {parameter.label} + + +
+ )} +
+ )); +} diff --git a/src/features/noise-reduction/down.wgsl b/src/features/noise-reduction/down.wgsl new file mode 100644 index 00000000..8413fb98 --- /dev/null +++ b/src/features/noise-reduction/down.wgsl @@ -0,0 +1,26 @@ +// Halves the spectrum by averaging each 2 × 2 block, which halves its noise; past an odd edge the last +// texel repeats. +import { pack, unpack } from "./texels.wgsl"; + +struct Sizes { + source: vec2u, + output: vec2u, +} + +@group(0) @binding(0) var sizes: Sizes; +@group(0) @binding(1) var source: array; +@group(0) @binding(2) var output: array; + +fn load(p: vec2u) -> vec4f { + let q = min(p, sizes.source - 1u); + return unpack(source[q.y * sizes.source.x + q.x]); +} + +@compute @workgroup_size(8, 8) fn down(@builtin(global_invocation_id) id: vec3u) { + if (any(id.xy >= sizes.output)) { + return; + } + let p = id.xy * 2u; + let sum = load(p) + load(p + vec2u(1u, 0u)) + load(p + vec2u(0u, 1u)) + load(p + 1u); + output[id.y * sizes.output.x + id.x] = pack(sum * 0.25); +} diff --git a/src/features/noise-reduction/edits.ts b/src/features/noise-reduction/edits.ts new file mode 100644 index 00000000..0ea75b10 --- /dev/null +++ b/src/features/noise-reduction/edits.ts @@ -0,0 +1,32 @@ +import { type EditorDocument, editLayer } from "@/core/document"; +import { type NoiseReduction, noNoiseReduction, reducible } from "@/core/image"; +import { change, parse } from "@/lib/parse"; +import { noiseReductionSchema } from "./model"; + +const noiseReductionChange = change(noiseReductionSchema); + +/** Sets the image layer's noise reduction, keeping the strengths `values` omits. */ +export function setNoiseReduction( + document: EditorDocument, + values: Partial, +) { + const [image] = document.scene.getState().layers; + if (!reducible(document.resources.get(image.source))) { + throw Error("Only RAW photos with a 2 × 2 mosaic take noise reduction."); + } + const strengths = parse( + noiseReductionChange, + values, + "Invalid noise reduction", + ); + editLayer(document, image.id, (layer) => { + if (layer.kind !== "image") { + throw Error("Select the image layer."); + } + const noiseReduction = { + ...(layer.noiseReduction ?? noNoiseReduction), + ...strengths, + }; + return { ...layer, noiseReduction }; + }); +} diff --git a/src/features/noise-reduction/finish-color.wgsl b/src/features/noise-reduction/finish-color.wgsl new file mode 100644 index 00000000..4e5f5675 --- /dev/null +++ b/src/features/noise-reduction/finish-color.wgsl @@ -0,0 +1,11 @@ +// A filtered color spectrum's two differences, as each 2 × 2 cell's mean: a reduced image's color +// holds no detail at that scale, and half its size takes a quarter of the memory. +import { unpack } from "./texels.wgsl"; + +@group(0) @binding(0) var size: vec2u; +@group(0) @binding(1) var spectrum: array; + +@fragment fn fs_main(@builtin(position) position: vec4f) -> @location(0) vec4f { + let cell = vec2u(position.xy); + return vec4f(unpack(spectrum[cell.y * size.x + cell.x]).xy * 0.5, 0.0, 1.0); +} diff --git a/src/features/noise-reduction/finish-light.wgsl b/src/features/noise-reduction/finish-light.wgsl new file mode 100644 index 00000000..00307891 --- /dev/null +++ b/src/features/noise-reduction/finish-light.wgsl @@ -0,0 +1,13 @@ +// A filtered light spectrum back to full size: each 2 × 2 cell's pixels from its Haar transform. +import { haar } from "./opponent.wgsl"; +import { unpack } from "./texels.wgsl"; + +@group(0) @binding(0) var size: vec2u; +@group(0) @binding(1) var spectrum: array; + +@fragment fn fs_main(@builtin(position) position: vec4f) -> @location(0) vec4f { + let p = vec2u(position.xy); + let cell = p / 2u; + let pixels = haar * unpack(spectrum[cell.y * size.x + cell.x]); + return vec4f(pixels[(p.y & 1u) * 2u + (p.x & 1u)], 0.0, 0.0, 1.0); +} diff --git a/src/features/noise-reduction/finish.wgsl b/src/features/noise-reduction/finish.wgsl new file mode 100644 index 00000000..85b70e10 --- /dev/null +++ b/src/features/noise-reduction/finish.wgsl @@ -0,0 +1,22 @@ +// How a filtered spectrum changes the mosaic: each 2 × 2 cell's four samples by position, as the +// unbiased inverse returns them, less the decoded ones. +import { fromSpectrum, Noise, unstabilize } from "./stabilize.wgsl"; +import { unpack } from "./texels.wgsl"; + +@group(0) @binding(0) var noise: Noise; +@group(0) @binding(1) var size: vec2u; +@group(0) @binding(2) var spectrum: array; +@group(0) @binding(3) var samples: texture_2d; + +@fragment fn fs_main(@builtin(position) position: vec4f) -> @location(0) vec4f { + let cell = vec2u(position.xy); + let p = vec2i(cell * 2u); + let decoded = vec4f( + f32(textureLoad(samples, p, 0).r), + f32(textureLoad(samples, p + vec2i(1, 0), 0).r), + f32(textureLoad(samples, p + vec2i(0, 1), 0).r), + f32(textureLoad(samples, p + vec2i(1, 1), 0).r), + ) - noise.black; + let z = fromSpectrum(unpack(spectrum[cell.y * size.x + cell.x]), noise); + return unstabilize(z, noise) - decoded; +} diff --git a/src/features/noise-reduction/fuse.wgsl b/src/features/noise-reduction/fuse.wgsl new file mode 100644 index 00000000..8fa4c72e --- /dev/null +++ b/src/features/noise-reduction/fuse.wgsl @@ -0,0 +1,77 @@ +// Replaces an estimate's low frequencies with those of the half-size spectrum's estimate, whose +// patches were large enough to tell broad noise from the image (after Pierazzo, Morel, and Facciolo, +// 2017): the difference between that coarser estimate and the estimate's own 2 × 2 averages, +// blurred at the coarse size so the coarser patches' finest frequencies, which ring at edges, stay +// out, then upsampled and added. Where the two levels agree, the estimate stays whole. +import { pack, unpack } from "./texels.wgsl"; + +struct Sizes { + fine: vec2u, + coarse: vec2u, +} + +// A Gaussian of 0.625 coarse texels, 1.25 fine ones, normalized over its five taps. +const taps = array(0.6377, 0.1773, 0.0038); + +@group(0) @binding(0) var sizes: Sizes; +@group(0) @binding(1) var fine: array; +// The fine estimate's 2 × 2 averages, and the coarser level's estimate. +@group(0) @binding(2) var average: array; +@group(0) @binding(3) var coarse: array; +// The difference blurred across rows, then down columns, at the coarse size. +@group(0) @binding(4) var across: array; +@group(0) @binding(5) var blurred: array; +@group(0) @binding(6) var output: array; + +fn coarseIndex(p: vec2i) -> u32 { + let q = vec2u(clamp(p, vec2i(0), vec2i(sizes.coarse) - 1)); + return q.y * sizes.coarse.x + q.x; +} + +fn difference(p: vec2i) -> vec4f { + let index = coarseIndex(p); + return unpack(coarse[index]) - unpack(average[index]); +} + +@compute @workgroup_size(8, 8) fn blurRows(@builtin(global_invocation_id) id: vec3u) { + if (any(id.xy >= sizes.coarse)) { + return; + } + let p = vec2i(id.xy); + var sum = taps[0] * difference(p); + for (var k = 1; k < 3; k++) { + sum += taps[k] * (difference(p - vec2i(k, 0)) + difference(p + vec2i(k, 0))); + } + across[coarseIndex(p)] = pack(sum); +} + +@compute @workgroup_size(8, 8) fn blurColumns(@builtin(global_invocation_id) id: vec3u) { + if (any(id.xy >= sizes.coarse)) { + return; + } + let p = vec2i(id.xy); + var sum = taps[0] * unpack(across[coarseIndex(p)]); + for (var k = 1; k < 3; k++) { + sum += taps[k] * (unpack(across[coarseIndex(p - vec2i(0, k))]) + unpack(across[coarseIndex(p + vec2i(0, k))])); + } + blurred[coarseIndex(p)] = pack(sum); +} + +fn blurredAt(p: vec2i) -> vec4f { + return unpack(blurred[coarseIndex(p)]); +} + +// The estimate plus the blurred difference, bilinear, each coarse texel centered on the 2 × 2 block +// it averages. +@compute @workgroup_size(8, 8) fn add(@builtin(global_invocation_id) id: vec3u) { + if (any(id.xy >= sizes.fine)) { + return; + } + let position = vec2f(id.xy) * 0.5 - 0.25; + let base = vec2i(floor(position)); + let f = position - floor(position); + let top = mix(blurredAt(base), blurredAt(base + vec2i(1, 0)), f.x); + let bottom = mix(blurredAt(base + vec2i(0, 1)), blurredAt(base + 1), f.x); + let index = id.y * sizes.fine.x + id.x; + output[index] = pack(unpack(fine[index]) + mix(top, bottom, f.y)); +} diff --git a/src/features/noise-reduction/image.ts b/src/features/noise-reduction/image.ts new file mode 100644 index 00000000..14b7deaf --- /dev/null +++ b/src/features/noise-reduction/image.ts @@ -0,0 +1,186 @@ +import { + type Compute, + compute, + type Effect, + effect, + frame, + type Gpu, + sampler, + type Target, + target, +} from "vgpu"; +import type { EncodedImage, NoiseReduction, Reduction } from "@/core/image"; +import { weakMemo } from "@/lib/weak-memo"; +import composeShader from "./compose-image.wgsl"; +import finishColorShader from "./finish-color.wgsl"; +import finishLightShader from "./finish-light.wgsl"; +import { componentLimits } from "./model"; +import prepareShader from "./prepare-image.wgsl"; +import { + block, + type Component, + groups, + half, + reduceSpectrum, + type Size, +} from "./pyramid"; + +/** The noise 8-bit samples hold at least, their rounding, in the encoding's units. */ +const rounding = 1 / 255 / Math.sqrt(12); + +/** + * One of the spectra an image reduces apart: how it is made and returns, at what size and format, how + * its components filter, and their noise from the finest level's measures. + */ +type Part = { + prepare: Compute; + finish: Effect; + size: (image: Size) => Size; + format: GPUTextureFormat; + components: readonly Component[]; + noise: (measured: readonly number[]) => readonly number[]; +}; + +/** + * Light returns at full size from its Haar components, which filter alike; white noise holds the + * same deviation in each, so the least measure, which texture raises least, gives all of theirs. + * Color returns at half size, which a reduced image's color needs no more of, and its two + * differences shrink twice as hard as light. + */ +const passes = weakMemo((gpu: Gpu) => ({ + light: { + prepare: compute(gpu, prepareShader, { entry: "light" }), + finish: effect(gpu, finishLightShader), + size: (image: Size) => image, + format: "r16float", + components: Array(4).fill({ + threshold: 3, + caution: 1, + mean: true, + }), + noise: (measured: readonly number[]) => { + // A component without detail anywhere, such as an enlarged image's differences, measures none. + const measures = measured.filter((value) => value > 0); + const least = measures.length ? Math.min(...measures) : 0; + return measured.map(() => Math.max(least, rounding)); + }, + } satisfies Part, + color: { + prepare: compute(gpu, prepareShader, { entry: "color" }), + finish: effect(gpu, finishColorShader), + size: half, + format: "rg16float", + components: Array(2).fill({ + threshold: 3, + caution: 2, + mean: true, + }), + noise: (measured: readonly number[]) => + measured.map((value) => Math.max(value, rounding)), + } satisfies Part, + compose: effect(gpu, composeShader), +})); + +/** + * Composes the image at any strengths from its reduced light and color, into a target of its own + * in the image's format: 8-bit images stay 8-bit, dithered. + */ +function createComposer( + gpu: Gpu, + { image }: EncodedImage, + light: Target, + color: Target, + noise: readonly number[], +) { + const output = target(gpu, { size: image.size, format: image.format }); + const quantum = image.format === "rgba8unorm-srgb" ? 1 / 255 : 0; + return { + render(strengths: NoiseReduction) { + const pass = passes(gpu).compose.set({ + limits: componentLimits(strengths), + params: { noise, quantum }, + source: image.color, + linearSampler: sampler(gpu, { + minFilter: "linear", + magFilter: "linear", + }), + light: light.color, + color: color.color, + }); + frame(gpu, (f) => f.pass(output, pass)); + return output; + }, + dispose: () => output.color.dispose(), + }; +} + +/** Reduces one of an image's spectra into a target of its own, with the noise it measured per component. */ +async function reducePart( + gpu: Gpu, + image: Target, + part: Part, + signal: AbortSignal, +) { + const size = half(image.size); + const { estimate, noise } = await reduceSpectrum( + gpu, + size, + (finest) => + part.prepare + .set({ source: image.color, size, output: finest.buffer }) + .dispatch(...groups(size)), + part.noise, + part.components, + signal, + ); + try { + const pass = part.finish.set({ size, spectrum: estimate.buffer }); + const output = target(gpu, { + size: part.size(image.size), + format: part.format, + }); + try { + frame(gpu, (f) => f.pass(output, pass)); + } catch (error) { + output.color.dispose(); + throw error; + } + return { output, noise }; + } finally { + estimate.buffer.dispose(); + } +} + +/** + * Reduces an image's noise in its encoded light and color, each at half size and apart, with no + * noise model: each level's noise comes from its flattest blocks. A composer then gives the image at + * any strengths at once. Color reduces first, so only its smaller result waits through light's. Stops + * between levels once `signal` aborts. + */ +export async function reduceImage( + gpu: Gpu, + encoded: EncodedImage, + signal: AbortSignal, +): Promise { + const { image } = encoded; + if (Math.min(...half(image.size)) < block) { + throw Error("The image is too small to reduce its noise."); + } + const parts = passes(gpu); + const color = await reducePart(gpu, image, parts.color, signal); + try { + const light = await reducePart(gpu, image, parts.light, signal); + const noise = [light.noise[0], color.noise[0], color.noise[1], 0]; + return { + compose: () => + createComposer(gpu, encoded, light.output, color.output, noise), + dispose() { + light.output.color.dispose(); + color.output.color.dispose(); + }, + }; + } catch (error) { + color.output.color.dispose(); + throw error; + } +} diff --git a/src/features/noise-reduction/limit.wgsl b/src/features/noise-reduction/limit.wgsl new file mode 100644 index 00000000..e8c870cb --- /dev/null +++ b/src/features/noise-reduction/limit.wgsl @@ -0,0 +1,20 @@ +// How much of the difference between an image and its reduction, in noise deviations, to remove: +// none at a share of 0 and all of it at 1. Between, a component loses nearly all of a difference +// whose neighborhood varies by less than a limit, mostly noise, and of one that varies more, mostly +// detail, about the limit's worth; the limit grows with the share, through its value at one half, +// and without bound toward 1. Weighing neighborhoods rather than single values keeps noise from +// standing out where it varies most. + +// Each component's share, and its limit at a share of one half. +export struct Limits { + shares: vec4f, + halves: vec4f, +} + +// `energy` is the difference's mean square around it. +export fn removed(difference: vec4f, energy: vec4f, limits: Limits) -> vec4f { + let limit = limits.halves * limits.shares; + let kept = 1.0 - limits.shares; + let part = difference * limit * inverseSqrt(limit * limit + kept * kept * energy); + return select(vec4f(0.0), part, limits.shares > vec4f(0.0)); +} diff --git a/src/features/noise-reduction/model.ts b/src/features/noise-reduction/model.ts new file mode 100644 index 00000000..9372b03f --- /dev/null +++ b/src/features/noise-reduction/model.ts @@ -0,0 +1,24 @@ +import { z } from "zod/mini"; +import type { NoiseReduction } from "@/core/image"; +import { range } from "@/lib/parse"; + +/** How much noise reduction removes from light and from color, each 0 to 100; 0 keeps it. */ +export const noiseReductionSchema = z.object({ + luminance: range(0, 100), + color: range(0, 100), +}) satisfies z.ZodMiniType; + +/** + * Per spectrum component, its share of the reduction, 0 to 1, and its limit at a share of one half, + * in noise deviations: less for detail the filter removes more of, more for color, which holds + * little. Light follows luminance, the two color differences color, and the greens' difference, which + * demosaicing draws as a maze, the stronger of the two. + */ +export function componentLimits({ luminance, color }: NoiseReduction) { + return { + shares: [luminance, color, color, Math.max(luminance, color)].map( + (amount) => amount / 100, + ), + halves: [1, 3, 3, 0.7], + }; +} diff --git a/src/features/noise-reduction/mosaic.ts b/src/features/noise-reduction/mosaic.ts new file mode 100644 index 00000000..476ff2d2 --- /dev/null +++ b/src/features/noise-reduction/mosaic.ts @@ -0,0 +1,148 @@ +import { compute, effect, frame, type Gpu, type Target, target } from "vgpu"; +import type { Mosaic, NoiseReduction, Reduction } from "@/core/image"; +import { weakMemo } from "@/lib/weak-memo"; +import composeShader from "./compose.wgsl"; +import finishShader from "./finish.wgsl"; +import { componentLimits } from "./model"; +import { measureNoise, type NoiseModel } from "./noise"; +import prepareShader from "./prepare.wgsl"; +import { + block, + type Component, + groups, + modeled, + reduceSpectrum, + type Size, +} from "./pyramid"; + +const passes = weakMemo((gpu: Gpu) => ({ + prepare: compute(gpu, prepareShader), + finish: effect(gpu, finishShader), + compose: effect(gpu, composeShader), +})); + +/** + * How the filter treats light, the two color differences, and the greens' difference: color shrinks + * twice as hard as light. The greens' difference holds less detail, so a coefficient needs more to + * stand out and shrinks eight times as hard; both greens see the same light, so its patches' means + * are noise or imbalance, which demosaicing would draw as a maze. + */ +const components: Component[] = [ + { threshold: 3, caution: 1, mean: true }, + { threshold: 3, caution: 2, mean: true }, + { threshold: 3, caution: 2, mean: true }, + { threshold: 3.5, caution: 8, mean: false }, +]; + +/** + * The transform's largest value, in noise deviations: half floats keep it, and the filter's + * fixed-point sums of its patches stay in range. + */ +const largest = 1024; + +/** + * The least floor that keeps the transform of a sample `range` above black within `largest`: solving + * 2r / (√(a·r + c) + √c) = L for √c, with c = 3/8 · a² + b. A photo with little noise, which would + * stretch the transform further, is filtered as if it had that much. + */ +function leastFloor(gain: number, range: number) { + const root = range / largest - (gain * largest) / 4; + return root > 0 ? root * root - 0.375 * gain * gain : 0; +} + +/** The filter's noise model per 2 × 2 position, and the positions of red, the two greens, and blue. */ +function noiseParams( + { pattern, black, white }: Mosaic, + { gain, floor }: NoiseModel, +) { + const greens = [0, 1, 2, 3].filter((p) => pattern[p] % 2 === 1); + const blacks = pattern.map((color) => black[color]); + return { + order: [pattern.indexOf(0), greens[0], greens[1], pattern.indexOf(2)], + black: blacks, + gain, + floor: floor.map((b, p) => + Math.max(b, leastFloor(gain[p], white - blacks[p])), + ), + }; +} + +type NoiseParams = ReturnType; + +/** Composes replacement samples at any strengths from the reduction's change, into a target of its own. */ +function createComposer( + gpu: Gpu, + mosaic: Mosaic, + noise: NoiseParams, + reduction: Target, +) { + const output = target(gpu, { size: reduction.size, format: "rgba16float" }); + return { + render(strengths: NoiseReduction) { + const pass = passes(gpu).compose.set({ + noise, + limits: componentLimits(strengths), + samples: mosaic.samples, + reduction: reduction.color, + }); + frame(gpu, (f) => f.pass(output, pass)); + return output; + }, + dispose: () => output.color.dispose(), + }; +} + +/** + * Reduces a mosaic's noise before demosaicing: its own noise model makes the noise unit Gaussian, a + * pyramid of patch filtering removes it from fine detail to broad color blotches, and the unbiased + * inverse returns samples, kept as their change to the decoded ones. A composer then gives + * replacement samples at any strengths at once: a half-size `rgba16float` target holding each 2 × 2 + * cell's four samples by position, above black. Stops between levels once `signal` aborts. + */ +export async function reduceMosaic( + gpu: Gpu, + mosaic: Mosaic, + signal: AbortSignal, +): Promise { + const { prepare, finish } = passes(gpu); + const size: Size = [ + Math.floor(mosaic.samples.width / 2), + Math.floor(mosaic.samples.height / 2), + ]; + if (Math.min(...size) < block) { + throw Error("The photo is too small to reduce its noise."); + } + const noise = noiseParams(mosaic, await measureNoise(gpu, mosaic)); + const { estimate } = await reduceSpectrum( + gpu, + size, + (finest) => + prepare + .set({ samples: mosaic.samples, noise, size, output: finest.buffer }) + .dispatch(...groups(size)), + modeled, + components, + signal, + ); + try { + const pass = finish.set({ + noise, + size, + spectrum: estimate.buffer, + samples: mosaic.samples, + }); + const change = target(gpu, { size, format: "rgba16float" }); + try { + frame(gpu, (f) => f.pass(change, pass)); + } catch (error) { + change.color.dispose(); + throw error; + } + return { + compose: () => createComposer(gpu, mosaic, noise, change), + dispose: () => change.color.dispose(), + }; + } finally { + estimate.buffer.dispose(); + } +} diff --git a/src/features/noise-reduction/noise.ts b/src/features/noise-reduction/noise.ts new file mode 100644 index 00000000..649d64ef --- /dev/null +++ b/src/features/noise-reduction/noise.ts @@ -0,0 +1,132 @@ +import { compute, type Gpu } from "vgpu"; +import type { Mosaic } from "@/core/image"; +import { weakMemo } from "@/lib/weak-memo"; +import shader from "./noise.wgsl"; + +type Vec4 = [number, number, number, number]; + +/** + * Sensor noise per 2 × 2 cell position: the variance of a sample `x` above black is + * `gain · x + floor`, photon noise growing with the signal over the read noise. + */ +export type NoiseModel = { gain: Vec4; floor: Vec4 }; + +const measure = weakMemo((gpu: Gpu) => + compute(gpu, shader, { entry: "measure" }), +); + +/** Sensor pixels along each block's side. */ +const block = 16; +/** Equal-count bins of blocks along the signal, each of which yields one point to fit. */ +const bins = 40; +/** + * The variance of the quietest tenth of a bin's blocks comes from its flattest ones, which texture + * reaches least; a block's estimate over 16 details is chi-squared, whose tenth percentile lies at + * 0.582 of its mean. + */ +const quiet = 0.1; +const quietBias = 0.582; +/** + * Texture only adds variance, so a bin this far above the fitted line holds texture rather than + * noise; fitting again without such bins follows the noise where a detailed photo has few flat blocks. + */ +const envelope = 1.3; + +function fitLine(points: readonly [number, number][]) { + // Weighted by the inverse variance squared, so dark and bright bins count by relative error. + let sw = 0; + let sx = 0; + let sy = 0; + let sxx = 0; + let sxy = 0; + for (const [x, y] of points) { + const w = 1 / (y * y); + sw += w; + sx += w * x; + sy += w * y; + sxx += w * x * x; + sxy += w * x * y; + } + const gain = (sw * sxy - sx * sy) / (sw * sxx - sx * sx); + const floor = (sy - gain * sx) / sw; + if (!(gain > 0)) return { gain: 0, floor: sy / sw }; + if (!(floor > 0)) return { gain: sxy / sxx, floor: 0 }; + return { gain, floor }; +} + +/** + * Fits one position's model to its blocks' means and variances, skipping clipped blocks; without + * any, the photo shows no noise. + */ +export function fitNoise( + means: readonly number[], + variances: readonly number[], +) { + const usable = means + .map((mean, i) => [mean, variances[i]] as [number, number]) + .filter(([, variance]) => variance > 0) + .sort(([a], [b]) => a - b); + if (!usable.length) { + return { gain: 0, floor: 0 }; + } + // A small photo fits fewer points, down to one, which leaves read noise alone. + const count = Math.min(bins, Math.ceil(usable.length / 16)); + const points: [number, number][] = []; + for (let i = 0; i < count; i++) { + const bin = usable.slice( + Math.floor((i * usable.length) / count), + Math.floor(((i + 1) * usable.length) / count), + ); + const variances = bin.map(([, variance]) => variance).sort((a, b) => a - b); + const mean = bin[Math.floor(bin.length / 2)][0]; + const variance = variances[Math.floor(bin.length * quiet)]; + points.push([mean, variance / quietBias]); + } + let fit = fitLine(points); + for (let pass = 0; pass < 5; pass++) { + const { gain, floor } = fit; + const below = points.filter(([x, y]) => y <= envelope * (gain * x + floor)); + if (below.length < 3) break; + fit = fitLine(below); + } + return fit; +} + +/** Measures the mosaic's noise on the GPU from its flat blocks, then fits each position's model. */ +export async function measureNoise( + gpu: Gpu, + { samples, pattern, black, white }: Mosaic, +): Promise { + const columns = Math.floor(samples.width / block); + const rows = Math.floor(samples.height / block); + const stats = gpu.device.createBuffer({ + size: columns * rows * 32, + usage: ["storage", "copy_src"], + }); + try { + measure(gpu) + .set({ + samples, + params: { black: pattern.map((color) => black[color]), white }, + blocks: stats, + }) + .dispatch(columns, rows); + const values = new Float32Array(await stats.read(columns * rows * 32)); + const gain: Vec4 = [0, 0, 0, 0]; + const floor: Vec4 = [0, 0, 0, 0]; + for (let position = 0; position < 4; position++) { + const means: number[] = []; + const variances: number[] = []; + for (let i = 0; i < values.length; i += 8) { + means.push(values[i + position]); + variances.push(values[i + 4 + position]); + } + const fit = fitNoise(means, variances); + gain[position] = fit.gain; + floor[position] = fit.floor; + } + return { gain, floor }; + } finally { + stats.dispose(); + } +} diff --git a/src/features/noise-reduction/noise.wgsl b/src/features/noise-reduction/noise.wgsl new file mode 100644 index 00000000..14a3a4f4 --- /dev/null +++ b/src/features/noise-reduction/noise.wgsl @@ -0,0 +1,64 @@ +// Noise statistics of a Bayer mosaic: for each block of 16 × 16 sensor pixels, each 2 × 2 cell +// position's mean above black and the variance of its diagonal Haar detail, which flat blocks owe to +// noise alone; small blocks leave more of them flat in a detailed photo. A position with a clipped +// sample in the block reports a negative variance, which the fit skips. +const side = 4u; +const threads = side * side; + +struct Params { + black: vec4f, + white: f32, +} + +@group(0) @binding(0) var samples: texture_2d; +@group(0) @binding(1) var params: Params; +// Per block: the four positions' means, then their variances. +@group(0) @binding(2) var blocks: array; + +var means: array; +var details: array; +var clipped: array; + +// The 2 × 2 cell at `cell`, one sample per position, above black. +fn cell(cell: vec2u) -> vec4f { + let p = vec2i(cell * 2u); + return vec4f( + f32(textureLoad(samples, p, 0).r), + f32(textureLoad(samples, p + vec2i(1, 0), 0).r), + f32(textureLoad(samples, p + vec2i(0, 1), 0).r), + f32(textureLoad(samples, p + vec2i(1, 1), 0).r), + ) - params.black; +} + +@compute @workgroup_size(side, side) fn measure( + @builtin(workgroup_id) block: vec3u, + @builtin(local_invocation_id) local: vec3u, + @builtin(local_invocation_index) thread: u32, + @builtin(num_workgroups) blocks_: vec3u, +) { + // Each thread takes 2 × 2 cells: their mean and their diagonal detail, (a - b - c + d) / 2. + let origin = (block.xy * side + local.xy) * 2u; + let a = cell(origin); + let b = cell(origin + vec2u(1u, 0u)); + let c = cell(origin + vec2u(0u, 1u)); + let d = cell(origin + vec2u(1u, 1u)); + let detail = (a - b - c + d) * 0.5; + means[thread] = (a + b + c + d) * 0.25; + details[thread] = detail * detail; + let top = max(max(a, b), max(c, d)) + params.black; + clipped[thread] = vec4u(top >= vec4f(params.white)); + workgroupBarrier(); + for (var stride = threads / 2u; stride > 0u; stride /= 2u) { + if (thread < stride) { + means[thread] += means[thread + stride]; + details[thread] += details[thread + stride]; + clipped[thread] = max(clipped[thread], clipped[thread + stride]); + } + workgroupBarrier(); + } + if (thread == 0u) { + let index = (block.y * blocks_.x + block.x) * 2u; + blocks[index] = means[0] / f32(threads); + blocks[index + 1u] = select(details[0] / f32(threads), vec4f(-1.0), clipped[0] != vec4u(0u)); + } +} diff --git a/src/features/noise-reduction/opponent.wgsl b/src/features/noise-reduction/opponent.wgsl new file mode 100644 index 00000000..8e3b7c29 --- /dev/null +++ b/src/features/noise-reduction/opponent.wgsl @@ -0,0 +1,37 @@ +// An image is filtered in the sRGB transfer's encoding, where noise varies far less with light than +// in linear, extended past 0 and 1 by symmetry and its own curve; and in an orthonormal basis of +// light and two color differences, so noise keeps its deviation in each. + +export fn encode(x: vec3f) -> vec3f { + let a = abs(x); + return sign(x) * select(1.055 * pow(a, vec3f(1.0 / 2.4)) - 0.055, 12.92 * a, a <= vec3f(0.0031308)); +} + +export fn decode(e: vec3f) -> vec3f { + let a = abs(e); + return sign(e) * select(pow((a + 0.055) / 1.055, vec3f(2.4)), a / 12.92, a <= vec3f(0.04045)); +} + +// Light, red against blue, and green against magenta. +const opponent = mat3x3f( + 0.57735027, 0.70710678, 0.40824829, + 0.57735027, 0.0, -0.81649658, + 0.57735027, -0.70710678, 0.40824829, +); + +export fn toOpponent(rgb: vec3f) -> vec3f { + return opponent * rgb; +} + +export fn fromOpponent(value: vec3f) -> vec3f { + return value * opponent; +} + +// Each 2 × 2 cell's four values, row-major, and its mean, horizontal, vertical, and diagonal +// differences, one to the other: orthonormal, so noise keeps its deviation, and its own inverse. +export const haar = mat4x4f( + 0.5, 0.5, 0.5, 0.5, + 0.5, -0.5, 0.5, -0.5, + 0.5, 0.5, -0.5, -0.5, + 0.5, -0.5, -0.5, 0.5, +); diff --git a/src/features/noise-reduction/prepare-image.wgsl b/src/features/noise-reduction/prepare-image.wgsl new file mode 100644 index 00000000..4f3232fe --- /dev/null +++ b/src/features/noise-reduction/prepare-image.wgsl @@ -0,0 +1,37 @@ +// An image's pixels as the filter's spectra at half size, from the light and two color differences of +// its encoded color: each 2 × 2 cell's light by its Haar transform, or its color differences, summed +// orthonormally so noise keeps its deviation. A last odd row or column repeats the one before. +import { encode, haar, toOpponent } from "./opponent.wgsl"; +import { pack } from "./texels.wgsl"; + +@group(0) @binding(0) var source: texture_2d; +@group(0) @binding(1) var size: vec2u; +@group(0) @binding(2) var output: array; + +fn opponentAt(p: vec2u) -> vec3f { + let q = min(p, textureDimensions(source) - 1u); + return toOpponent(encode(textureLoad(source, q, 0).rgb)); +} + +@compute @workgroup_size(8, 8) fn light(@builtin(global_invocation_id) id: vec3u) { + if (any(id.xy >= size)) { + return; + } + let p = id.xy * 2u; + let cell = vec4f( + opponentAt(p).x, + opponentAt(p + vec2u(1u, 0u)).x, + opponentAt(p + vec2u(0u, 1u)).x, + opponentAt(p + 1u).x, + ); + output[id.y * size.x + id.x] = pack(haar * cell); +} + +@compute @workgroup_size(8, 8) fn color(@builtin(global_invocation_id) id: vec3u) { + if (any(id.xy >= size)) { + return; + } + let p = id.xy * 2u; + let sum = opponentAt(p).yz + opponentAt(p + vec2u(1u, 0u)).yz + opponentAt(p + vec2u(0u, 1u)).yz + opponentAt(p + 1u).yz; + output[id.y * size.x + id.x] = pack(vec4f(sum * 0.5, 0.0, 0.0)); +} diff --git a/src/features/noise-reduction/prepare.wgsl b/src/features/noise-reduction/prepare.wgsl new file mode 100644 index 00000000..bff73ba1 --- /dev/null +++ b/src/features/noise-reduction/prepare.wgsl @@ -0,0 +1,23 @@ +// Packs each 2 × 2 cell of the mosaic into one texel of the filter's spectrum, where noise is unit +// Gaussian everywhere. +import { Noise, stabilize, toSpectrum } from "./stabilize.wgsl"; +import { pack } from "./texels.wgsl"; + +@group(0) @binding(0) var samples: texture_2d; +@group(0) @binding(1) var noise: Noise; +@group(0) @binding(2) var size: vec2u; +@group(0) @binding(3) var output: array; + +@compute @workgroup_size(8, 8) fn prepare(@builtin(global_invocation_id) id: vec3u) { + if (any(id.xy >= size)) { + return; + } + let p = vec2i(id.xy * 2u); + let cell = vec4f( + f32(textureLoad(samples, p, 0).r), + f32(textureLoad(samples, p + vec2i(1, 0), 0).r), + f32(textureLoad(samples, p + vec2i(0, 1), 0).r), + f32(textureLoad(samples, p + vec2i(1, 1), 0).r), + ); + output[id.y * size.x + id.x] = pack(toSpectrum(stabilize(cell - noise.black, noise), noise)); +} diff --git a/src/features/noise-reduction/pyramid.ts b/src/features/noise-reduction/pyramid.ts new file mode 100644 index 00000000..0f454900 --- /dev/null +++ b/src/features/noise-reduction/pyramid.ts @@ -0,0 +1,291 @@ +import { type Buffer, compute, type Gpu } from "vgpu"; +import { weakMemo } from "@/lib/weak-memo"; +import downShader from "./down.wgsl"; +import fuseShader from "./fuse.wgsl"; +import shrinkShader from "./shrink.wgsl"; +import spreadShader from "./spread.wgsl"; + +export type Size = readonly [number, number]; +/** One level of the spectrum: four half floats per texel, row after row. */ +export type Spectrum = { buffer: Buffer; size: Size }; +/** How the filter treats one of a spectrum's components. */ +export type Component = { + /** Noise deviations past which the first step keeps a coefficient. */ + threshold: number; + /** How much more its noise counts in the second step's Wiener weights. */ + caution: number; + /** Whether patches keep their means whole, as a component that holds a level must to stay unbiased. */ + mean: boolean; +}; + +const passes = weakMemo((gpu: Gpu) => ({ + down: compute(gpu, downShader), + spread: compute(gpu, spreadShader), + shrink: compute(gpu, shrinkShader), + blurRows: compute(gpu, fuseShader, { entry: "blurRows" }), + blurColumns: compute(gpu, fuseShader, { entry: "blurColumns" }), + add: compute(gpu, fuseShader, { entry: "add" }), +})); + +/** Levels of the pyramid at most, and the shortest side a level may have. */ +const levels = 6; +const smallest = 32; +/** Spectrum texels along the side of the tile each workgroup of the filter takes. */ +const tile = 32; +/** Spectrum texels along the side of each block that measures a level's noise. */ +export const block = 16; +/** The flattest quarter of blocks gives the noise; their chi-squared estimates sit at 0.876 of it. */ +const quartileBias = 0.876; +/** + * How far a level's noise may exceed what white noise would leave there: lossy compression + * correlates noise, raising it at coarse levels; texture raises a measure far more, and is not noise. + */ +const margin = 1.15; + +function spectrum(gpu: Gpu, size: Size): Spectrum { + const buffer = gpu.device.createBuffer({ + size: size[0] * size[1] * 8, + usage: ["storage"], + }); + return { buffer, size }; +} + +export function groups([width, height]: Size, side = 8) { + return [Math.ceil(width / side), Math.ceil(height / side)] as const; +} + +export function half([width, height]: Size): Size { + return [Math.ceil(width / 2), Math.ceil(height / 2)]; +} + +/** Averages 2 × 2 texels of `source` into `output`, half its size. */ +function down(gpu: Gpu, source: Spectrum, output: Spectrum) { + passes(gpu) + .down.set({ + sizes: { source: source.size, output: output.size }, + source: source.buffer, + output: output.buffer, + }) + .dispatch(...groups(output.size)); +} + +/** Adds averages of `pyramid`'s last level, each half the one before, down to `smallest` texels. */ +function extendPyramid(gpu: Gpu, pyramid: Spectrum[]) { + while (pyramid.length < levels) { + const finer = pyramid[pyramid.length - 1]; + const size = half(finer.size); + if (Math.min(...size) < smallest) return; + const level = spectrum(gpu, size); + pyramid.push(level); + down(gpu, finer, level); + } +} + +/** + * Each component's noise deviation in one level, from its flattest blocks, which texture can only + * raise. Blocks without any detail, such as clipped highlights, hold no noise to measure; without + * others, a component measures none. + */ +async function measureSpread(gpu: Gpu, { buffer, size }: Spectrum) { + const columns = Math.floor(size[0] / block); + const rows = Math.floor(size[1] / block); + const blocks = gpu.device.createBuffer({ + size: columns * rows * 16, + usage: ["storage", "copy_src"], + }); + try { + passes(gpu) + .spread.set({ size, spectrum: buffer, blocks }) + .dispatch(columns, rows); + const values = new Float32Array(await blocks.read(columns * rows * 16)); + return [0, 1, 2, 3].map((component) => { + const variances = values + .filter((variance, i) => i % 4 === component && variance > 0) + .sort((a, b) => a - b); + const quartile = variances[Math.floor(variances.length / 4)]; + return quartile ? Math.sqrt(quartile / quartileBias) : 0; + }); + } finally { + blocks.dispose(); + } +} + +/** + * Each level's noise deviation per component, finest first. `finest` turns the finest level's + * measures into its noise. A coarser level's flattest blocks give its noise, up to a margin over what + * white noise would leave there, since they still hold texture. + */ +async function measureLevels( + gpu: Gpu, + pyramid: readonly Spectrum[], + finest: (measured: readonly number[]) => readonly number[], +) { + const measured = []; + for (const level of pyramid) measured.push(await measureSpread(gpu, level)); + const base = finest(measured[0]); + return measured.map((level, l) => + level.map((value, c) => { + const white = base[c] * 2 ** -l; + if (!l) return base[c]; + return value ? Math.min(value, margin * white) : white; + }), + ); +} + +/** The finest level's noise for a noise model that predicts unit noise, which lower measures correct. */ +export function modeled(measured: readonly number[]) { + return measured.map((value) => Math.min(value || 1, 1)); +} + +/** Each of four components' value, 0 past the components a spectrum holds. */ +const padded = (values: readonly number[]) => + [0, 1, 2, 3].map((c) => values[c] ?? 0); + +/** One step of the patch filter over `noisy` into `output`, shrinking by `pilot`'s Wiener weights if given. */ +function shrink( + gpu: Gpu, + noisy: Spectrum, + output: Spectrum, + sigma: readonly number[], + components: readonly Component[], + pilot?: Spectrum, +) { + const { size } = noisy; + passes(gpu) + .shrink.set({ + params: { + size, + channels: components.length, + wiener: Number(Boolean(pilot)), + sigma, + threshold: padded(components.map((c) => c.threshold)), + caution: padded(components.map((c) => c.caution)), + mean: padded(components.map((c) => Number(c.mean))), + }, + noisy: noisy.buffer, + pilot: (pilot ?? noisy).buffer, + output: output.buffer, + }) + .dispatch(...groups(size, tile)); +} + +/** Replaces `estimate`'s low frequencies with those of `coarse`, the estimate of the level below, into `output`. */ +function fuse( + gpu: Gpu, + estimate: Spectrum, + coarse: Spectrum, + output: Spectrum, +) { + const { blurRows, blurColumns, add } = passes(gpu); + const average = spectrum(gpu, coarse.size); + const across = spectrum(gpu, coarse.size); + const blurred = spectrum(gpu, coarse.size); + const sizes = { fine: estimate.size, coarse: coarse.size }; + try { + down(gpu, estimate, average); + blurRows + .set({ + sizes, + average: average.buffer, + coarse: coarse.buffer, + across: across.buffer, + }) + .dispatch(...groups(coarse.size)); + blurColumns + .set({ sizes, across: across.buffer, blurred: blurred.buffer }) + .dispatch(...groups(coarse.size)); + add + .set({ + sizes, + fine: estimate.buffer, + blurred: blurred.buffer, + output: output.buffer, + }) + .dispatch(...groups(estimate.size)); + } finally { + average.buffer.dispose(); + across.buffer.dispose(); + blurred.buffer.dispose(); + } +} + +/** + * Denoises one level, taking its low frequencies from the estimate of the level below when there is + * one: thresholding first, then Wiener shrinkage guided by that first estimate. + */ +function filterLevel( + gpu: Gpu, + noisy: Spectrum, + coarse: Spectrum | undefined, + sigma: readonly number[], + components: readonly Component[], +) { + const first = spectrum(gpu, noisy.size); + try { + const second = spectrum(gpu, noisy.size); + try { + shrink(gpu, noisy, first, sigma, components); + if (!coarse) { + shrink(gpu, noisy, second, sigma, components, first); + return second; + } + fuse(gpu, first, coarse, second); + shrink(gpu, noisy, first, sigma, components, second); + fuse(gpu, first, coarse, second); + return second; + } catch (error) { + second.buffer.dispose(); + throw error; + } + } finally { + first.buffer.dispose(); + } +} + +/** + * Denoises a spectrum of `size`, which `prepare` fills, through a pyramid of its averages from the + * coarsest level up: each level's `components` are filtered at the noise they measure, which + * `noiseOf` turns into the finest level's, and each level takes its low frequencies from the + * estimate of the level below. Resolves to the finest estimate and the finest level's noise per + * component; stops between steps once `signal` aborts. + */ +export async function reduceSpectrum( + gpu: Gpu, + size: Size, + prepare: (finest: Spectrum) => void, + noiseOf: (measured: readonly number[]) => readonly number[], + components: readonly Component[], + signal: AbortSignal, +) { + signal.throwIfAborted(); + const pyramid = [spectrum(gpu, size)]; + let estimate: Spectrum | undefined; + try { + prepare(pyramid[0]); + extendPyramid(gpu, pyramid); + const noise = await measureLevels(gpu, pyramid, noiseOf); + signal.throwIfAborted(); + for (let l = pyramid.length - 1; ; l--) { + const coarse = estimate; + const filtered = filterLevel( + gpu, + pyramid[l], + coarse, + noise[l], + components, + ); + estimate = filtered; + coarse?.buffer.dispose(); + pyramid[l].buffer.dispose(); + // Each level submits before the next, so other GPU work, such as the display, runs between. + await gpu.gpu.queue.onSubmittedWorkDone(); + signal.throwIfAborted(); + if (!l) return { estimate: filtered, noise: noise[0] }; + } + } catch (error) { + estimate?.buffer.dispose(); + throw error; + } finally { + for (const level of pyramid) level.buffer.dispose(); + } +} diff --git a/src/features/noise-reduction/reduce.ts b/src/features/noise-reduction/reduce.ts new file mode 100644 index 00000000..60651f67 --- /dev/null +++ b/src/features/noise-reduction/reduce.ts @@ -0,0 +1,23 @@ +import type { Gpu } from "vgpu"; +import type { ImageSource, Reduction } from "@/core/image"; +import { reduceImage } from "./image"; +import { reduceMosaic } from "./mosaic"; + +/** + * Reduces a source's noise once: a RAW photo on its mosaic before demosaicing, any other image in + * its encoded color. + */ +export async function reduceNoise( + gpu: Gpu, + source: ImageSource, + signal: AbortSignal, +): Promise { + const mosaic = source.raw?.mosaic; + if (mosaic) { + return reduceMosaic(gpu, mosaic, signal); + } + if (source.raw) { + throw Error("Only RAW photos with a 2 × 2 mosaic take noise reduction."); + } + return reduceImage(gpu, source, signal); +} diff --git a/src/features/noise-reduction/shrink.wgsl b/src/features/noise-reduction/shrink.wgsl new file mode 100644 index 00000000..417d27f5 --- /dev/null +++ b/src/features/noise-reduction/shrink.wgsl @@ -0,0 +1,153 @@ +// One step of sliding DCT denoising (Yu and Sapiro, 2011), the patch filter BM3D builds on (Dabov, +// Foi, Katkovnik, and Egiazarian, 2007), over a tile of the spectrum per workgroup: every 8 × 8 patch +// that starts on every other texel and reaches the tile is transformed, shrunk, and added back to the +// 64 texels it covers, weighted by how little of it survived, so 16 patches cover each texel. The +// first step keeps the coefficients that stand out of the noise; the second shrinks each by the +// Wiener weight the first step's estimate gives it. Patches are filtered alone rather than grouped +// with similar ones, which would mistake aligned noise for likeness. +import { pack, unpack } from "./texels.wgsl"; + +const side = 8u; +const tile = 32u; +const stride = 2u; +// Patches whose texels reach the tile, along each side. +const span = 19u; +const threads = 256u; +// Fixed-point scales of the sums. Light, the largest component, reaches twice the transform's largest +// value, 1024, and ringing at most doubles it; the 16 patches over a texel weigh 9 at most, in their +// Kaiser windows. Both sums stay in range. +const valueScale = 8192.0; +const weightScale = 16777216.0; +const kaiser = array(0.438676, 0.681324, 0.87684, 0.985823, 0.985823, 0.87684, 0.681324, 0.438676); + +struct Params { + size: vec2u, + // Components the spectrum holds, up to four. + channels: u32, + // 0 thresholds; 1 shrinks by the Wiener weights of `pilot`, the first step's estimate. + wiener: u32, + // Each component's noise deviation. + sigma: vec4f, + // Per component, the noise deviations past which the first step keeps a coefficient, how much more + // its noise counts in the second step's Wiener weights, and 1 where patches keep their means whole. + threshold: vec4f, + caution: vec4f, + mean: vec4u, +} + +@group(0) @binding(0) var params: Params; +@group(0) @binding(1) var noisy: array; +@group(0) @binding(2) var pilot: array; +@group(0) @binding(3) var output: array; + +// The 8-point DCT basis, [frequency][position]. +var basis: array; +var sums: array, 1024>; +var weights: array, 1024>; + +// A texel past the edge reads the image mirrored there, so patches cover the edge as fully as the rest. +fn mirrored(p: vec2i) -> u32 { + let last = vec2i(params.size) - 1; + let q = min(abs(p), 2 * last - p); + return u32(q.y) * params.size.x + u32(q.x); +} + +// The 2D DCT of a block, or its inverse, separably along rows and then columns. +fn transform(block: ptr>, inverse: bool) { + var line: array; + for (var step = 0u; step < 2u; step++) { + for (var a = 0u; a < 8u; a++) { + for (var k = 0u; k < 8u; k++) { + var sum = 0.0; + for (var i = 0u; i < 8u; i++) { + let b = select(basis[k * 8u + i], basis[i * 8u + k], inverse); + sum += b * (*block)[select(a * 8u + i, i * 8u + a, step == 1u)]; + } + line[k] = sum; + } + for (var k = 0u; k < 8u; k++) { + (*block)[select(a * 8u + k, k * 8u + a, step == 1u)] = line[k]; + } + } + } +} + +@compute @workgroup_size(16, 16) fn shrink( + @builtin(workgroup_id) cell: vec3u, + @builtin(local_invocation_index) t: u32, +) { + let frequency = (t / 8u) % 8u; + let position = t % 8u; + if (t < 64u) { + basis[t] = select(0.5, 0.35355339, frequency == 0u) * cos(f32((2u * position + 1u) * frequency) * 0.19634954); + } + let origin = vec2i(cell.xy * tile); + var results: array; + for (var channel = 0u; channel < params.channels; channel++) { + for (var i = t; i < tile * tile; i += threads) { + atomicStore(&sums[i], 0); + atomicStore(&weights[i], 0u); + } + workgroupBarrier(); + let sigma = params.sigma[channel]; + for (var i = t; i < span * span; i += threads) { + let corner = origin - i32(side - stride) + i32(stride) * vec2i(vec2u(i % span, i / span)); + var values: array; + var guide: array; + for (var j = 0u; j < 64u; j++) { + let index = mirrored(corner + vec2i(vec2u(j % 8u, j / 8u))); + values[j] = unpack(noisy[index])[channel]; + if (params.wiener == 1u) { + guide[j] = unpack(pilot[index])[channel]; + } + } + transform(&values, false); + if (params.wiener == 1u) { + transform(&guide, false); + } + // A patch's mean, its DC, stays whole where a component holds a level: the inverse transform + // expects every level unbiased. + var kept = 0.0; + for (var j = 0u; j < 64u; j++) { + if (j == 0u && params.mean[channel] == 1u) { + kept += 1.0; + } else if (params.wiener == 1u) { + let g = guide[j] * guide[j]; + let w = g / (g + params.caution[channel] * sigma * sigma); + values[j] *= w; + kept += w * w; + } else if (abs(values[j]) > params.threshold[channel] * sigma) { + kept += 1.0; + } else { + values[j] = 0.0; + } + } + transform(&values, true); + let weight = 1.0 / max(kept, 1.0); + for (var j = 0u; j < 64u; j++) { + let p = corner + vec2i(vec2u(j % 8u, j / 8u)) - origin; + if (all(p >= vec2i(0)) && all(p < vec2i(i32(tile)))) { + let w = weight * kaiser[j % 8u] * kaiser[j / 8u]; + let index = u32(p.y) * tile + u32(p.x); + atomicAdd(&sums[index], i32(round(values[j] * w * valueScale))); + atomicAdd(&weights[index], u32(round(w * weightScale))); + } + } + } + workgroupBarrier(); + for (var k = 0u; k < 4u; k++) { + let i = t + k * threads; + let value = f32(atomicLoad(&sums[i])) / valueScale; + let weight = f32(atomicLoad(&weights[i])) / weightScale; + results[k][channel] = value / weight; + } + workgroupBarrier(); + } + for (var k = 0u; k < 4u; k++) { + let i = t + k * threads; + let p = vec2u(origin) + vec2u(i % tile, i / tile); + if (all(p < params.size)) { + output[p.y * params.size.x + p.x] = pack(results[k]); + } + } +} diff --git a/src/features/noise-reduction/spread.wgsl b/src/features/noise-reduction/spread.wgsl new file mode 100644 index 00000000..80ee2ccc --- /dev/null +++ b/src/features/noise-reduction/spread.wgsl @@ -0,0 +1,36 @@ +// The noise left in one level of the spectrum: for each block of 16 × 16 texels, the variance of each +// component's diagonal Haar detail, which flat blocks owe to noise alone. +import { unpack } from "./texels.wgsl"; + +const side = 8u; + +@group(0) @binding(0) var size: vec2u; +@group(0) @binding(1) var spectrum: array; +@group(0) @binding(2) var blocks: array; + +var details: array; + +fn load(p: vec2u) -> vec4f { + return unpack(spectrum[p.y * size.x + p.x]); +} + +@compute @workgroup_size(side, side) fn measure( + @builtin(workgroup_id) block: vec3u, + @builtin(local_invocation_id) local: vec3u, + @builtin(local_invocation_index) thread: u32, + @builtin(num_workgroups) count: vec3u, +) { + let p = (block.xy * side + local.xy) * 2u; + let detail = (load(p) - load(p + vec2u(1u, 0u)) - load(p + vec2u(0u, 1u)) + load(p + 1u)) * 0.5; + details[thread] = detail * detail; + workgroupBarrier(); + for (var stride = 32u; stride > 0u; stride /= 2u) { + if (thread < stride) { + details[thread] += details[thread + stride]; + } + workgroupBarrier(); + } + if (thread == 0u) { + blocks[block.y * count.x + block.x] = details[0] / 64.0; + } +} diff --git a/src/features/noise-reduction/stabilize.wgsl b/src/features/noise-reduction/stabilize.wgsl new file mode 100644 index 00000000..4464ca1c --- /dev/null +++ b/src/features/noise-reduction/stabilize.wgsl @@ -0,0 +1,65 @@ +// The mosaic's noise model per 2 × 2 cell position, and the transforms between sensor samples and the +// filter's spectrum, where noise is unit Gaussian everywhere. + +export struct Noise { + // The positions of red, the two greens, and blue. + order: vec4u, + black: vec4f, + // A sample x above black varies by gain · x + floor. + gain: vec4f, + floor: vec4f, +} + +// Orthonormal, so unit noise stays unit noise: light, then green against magenta, red against blue, +// and the two greens' difference, which holds little more than noise. +const spectrum = mat4x4f( + 0.5, 0.5, 0.70710678, 0.0, + 0.5, -0.5, 0.0, 0.70710678, + 0.5, -0.5, 0.0, -0.70710678, + 0.5, 0.5, -0.70710678, 0.0, +); + +// The lowest value the transform takes, at x = −c/a, below which the model allows no sample: dead +// pixels lie there. Unbounded without photon noise. +fn lowest(noise: Noise) -> vec4f { + let root = sqrt(0.375 * noise.gain * noise.gain + noise.floor); + return select(vec4f(-3.4e38), -2.0 * root / noise.gain, noise.gain > vec4f(0.0)); +} + +// The generalized Anscombe transform 2/a · (√(a·x + c) − √c), with c = 3/8 · a² + b, of samples +// above black, written so it keeps its precision as the gain goes to zero. Zero at black, which +// keeps half floats precise. +export fn stabilize(x: vec4f, noise: Noise) -> vec4f { + let c = 0.375 * noise.gain * noise.gain + noise.floor; + return max(2.0 * x / (sqrt(max(noise.gain * x + c, vec4f(0.0))) + sqrt(c)), lowest(noise)); +} + +// The inverse that is unbiased for Poisson-Gaussian noise, in the closed form of Mäkitalo and Foi +// (2013): the algebraic inverse plus a correction in powers of 1 / D, where D is the transform +// measured from its own zero. +export fn unstabilize(value: vec4f, noise: Noise) -> vec4f { + let a = noise.gain; + let root = sqrt(0.375 * a * a + noise.floor); + let z = max(value, lowest(noise)); + let u = clamp(a / (a * z + 2.0 * root), vec4f(0.0), vec4f(1.0)); + let correction = 0.25 + u * (0.30618622 + u * (-1.375 + u * 0.76546554)); + return root * z + 0.25 * a * z * z + a * correction; +} + +// Stabilized samples by position into the spectrum. +export fn toSpectrum(z: vec4f, noise: Noise) -> vec4f { + let o = noise.order; + return spectrum * vec4f(z[o.x], z[o.y], z[o.z], z[o.w]); +} + +// The spectrum back to stabilized samples by position. +export fn fromSpectrum(value: vec4f, noise: Noise) -> vec4f { + let colors = value * spectrum; + let o = noise.order; + var z: vec4f; + z[o.x] = colors.x; + z[o.y] = colors.y; + z[o.z] = colors.z; + z[o.w] = colors.w; + return z; +} diff --git a/src/features/noise-reduction/texels.wgsl b/src/features/noise-reduction/texels.wgsl new file mode 100644 index 00000000..db7454fd --- /dev/null +++ b/src/features/noise-reduction/texels.wgsl @@ -0,0 +1,8 @@ +// The spectrum lives in storage buffers: four half floats per texel, row after row. +export fn unpack(texel: vec2u) -> vec4f { + return vec4f(unpack2x16float(texel.x), unpack2x16float(texel.y)); +} + +export fn pack(value: vec4f) -> vec2u { + return vec2u(pack2x16float(value.xy), pack2x16float(value.zw)); +} diff --git a/tests/fixtures/noisy.png b/tests/fixtures/noisy.png new file mode 100644 index 00000000..ab3c4d21 Binary files /dev/null and b/tests/fixtures/noisy.png differ diff --git a/tests/fixtures/raw/README.md b/tests/fixtures/raw/README.md index be95c2d4..d39443b6 100644 --- a/tests/fixtures/raw/README.md +++ b/tests/fixtures/raw/README.md @@ -5,3 +5,5 @@ These two synthetic MIT-licensed fixtures are copied unchanged from [raw-webgpu] `bayer.dng` stores RGGB samples; `linear-jxl.dng` stores equivalent linear RGB with JPEG XL compression. Both use an sRGB-like camera profile, samples `[32768, 16384, 8192]`, white level 65535 and orientation 6. Expected output is 96×128 with sRGB pixels approximately `[188, 137, 99]`. OpenLight tests loading, white-balance edits, preview/export agreement, history and resource lifetime with these files. + +`noisy-bayer.dng` is `bayer.dng` with seeded Poisson-Gaussian noise in its samples, at 64 sensor units per electron over a read noise of 300, for noise reduction. diff --git a/tests/fixtures/raw/noisy-bayer.dng b/tests/fixtures/raw/noisy-bayer.dng new file mode 100644 index 00000000..0d751341 Binary files /dev/null and b/tests/fixtures/raw/noisy-bayer.dng differ diff --git a/tests/noise-reduction.e2e.ts b/tests/noise-reduction.e2e.ts new file mode 100644 index 00000000..d1fb5f33 --- /dev/null +++ b/tests/noise-reduction.e2e.ts @@ -0,0 +1,104 @@ +import { expect, test } from "./fixtures"; +import { readImage } from "./images"; + +/** A grid of points inside a `width` × `height` image, clear of its edges. */ +function grid(width: number, height: number) { + return Array.from( + { length: 48 }, + (_, i) => + [ + Math.round(((i % 6) + 1) * (width / 7)), + Math.round((Math.floor(i / 6) + 1) * (height / 9)), + ] as const, + ); +} + +const mean = (values: number[]) => + values.reduce((sum, value) => sum + value, 0) / values.length; +const deviation = (values: number[]) => + Math.sqrt(mean(values.map((value) => (value - mean(values)) ** 2))); + +/** Light, the channels' mean, and color, red less blue, across samples of a flat photo. */ +function measure(samples: number[][]) { + const light = samples.map(([r, g, b]) => (r + g + b) / 3); + const color = samples.map(([r, , b]) => r - b); + return { + light: mean(light), + color: mean(color), + lightNoise: deviation(light), + colorNoise: deviation(color), + }; +} + +test("Noise reduction cleans a noisy RAW before demosaicing and returns it as decoded at 0", async ({ + page, +}) => { + await page.goto("/"); + await page.waitForFunction(() => window.openlight); + await page.evaluate( + (url) => window.openlight.loadUrl(url), + "/tests/fixtures/raw/noisy-bayer.dng", + ); + await expect( + page.getByRole("slider", { name: "Luminance noise" }), + ).toBeVisible(); + // Off its As Shot balance, a RAW develops into its own texture, from samples it writes back. + await page.evaluate(() => + window.openlight.setWhiteBalance({ temperature: 4000, tint: 20 }), + ); + const points = grid(96, 128); + const noisy = await readImage(page, undefined, points); + await page.evaluate(() => + window.openlight.setNoiseReduction({ luminance: 75, color: 75 }), + ); + const reduced = await readImage(page, undefined, points); + const before = measure(noisy.samples ?? []); + const after = measure(reduced.samples ?? []); + expect(after.lightNoise).toBeLessThan(before.lightNoise / 3); + expect(after.colorNoise).toBeLessThan(before.colorNoise / 3); + expect(Math.abs(after.light - before.light)).toBeLessThan(2); + // The reductions stay, so changing strengths needs no other, and 0 develops the decoded samples. + await page.evaluate(() => + window.openlight.setNoiseReduction({ luminance: 0, color: 0 }), + ); + expect((await readImage(page, undefined, points)).samples).toEqual( + noisy.samples, + ); +}); + +test("Noise reduction cleans an image's color and light apart", async ({ + page, +}) => { + await page.goto("/"); + await page.waitForFunction(() => window.openlight); + await page.evaluate( + (url) => window.openlight.loadUrl(url), + "/tests/fixtures/noisy.png", + ); + const points = grid(64, 96); + const noisy = await readImage(page, undefined, points); + const before = measure(noisy.samples ?? []); + await page.evaluate(() => + window.openlight.setNoiseReduction({ luminance: 0, color: 75 }), + ); + const colorOnly = measure( + (await readImage(page, undefined, points)).samples ?? [], + ); + expect(colorOnly.colorNoise).toBeLessThan(before.colorNoise / 2); + expect(colorOnly.lightNoise).toBeGreaterThan(before.lightNoise * 0.6); + expect(Math.abs(colorOnly.color - before.color)).toBeLessThan(3); + await page.evaluate(() => + window.openlight.setNoiseReduction({ luminance: 75 }), + ); + const both = measure( + (await readImage(page, undefined, points)).samples ?? [], + ); + expect(both.lightNoise).toBeLessThan(before.lightNoise / 2); + expect(Math.abs(both.light - before.light)).toBeLessThan(2); + await page.evaluate(() => + window.openlight.setNoiseReduction({ luminance: 0, color: 0 }), + ); + expect((await readImage(page, undefined, points)).samples).toEqual( + noisy.samples, + ); +}); diff --git a/tests/noise-reduction.test.ts b/tests/noise-reduction.test.ts new file mode 100644 index 00000000..a9c9dbe7 --- /dev/null +++ b/tests/noise-reduction.test.ts @@ -0,0 +1,30 @@ +import { expect, test } from "bun:test"; +import { fitNoise } from "@/features/noise-reduction/noise"; + +test("noise fits photon and read noise from the flattest blocks, past texture", () => { + // A seeded generator, so the blocks are the same every run. + let seed = 7; + const random = () => { + seed = (seed * 1664525 + 1013904223) >>> 0; + return seed / 2 ** 32; + }; + const normal = () => + Math.sqrt(-2 * Math.log(1 - random())) * Math.cos(2 * Math.PI * random()); + const means: number[] = []; + const variances: number[] = []; + for (let i = 0; i < 20000; i++) { + const mean = random() ** 2 * 4000; + // Each block averages 16 squared details: chi-squared with 16 degrees of freedom around the noise. + let sum = 0; + for (let k = 0; k < 16; k++) sum += normal() ** 2; + const noise = (32 * mean + 1500) * (sum / 16); + // A third of the blocks hold texture as well as noise. + const texture = i % 3 === 0 ? (1 + 10 * random()) * noise : 0; + means.push(mean); + variances.push(noise + texture); + } + // Texture in every bin lifts its quietest tenth a little, never to the textured blocks' level. + const { gain, floor } = fitNoise(means, variances); + expect(Math.abs(gain / 32 - 1)).toBeLessThan(0.15); + expect(Math.abs(floor / 1500 - 1)).toBeLessThan(0.15); +});