Imaging and computer vision

afield.media.image can look at a picture as well as make one. Its CPU imaging and computer-vision family has a single value type, Img, and around it a set of operations that turn a raster into structure, from convolution and gradients up through edges, corners, lines, blobs, pyramids, warps and motion. It is the raster counterpart to the time-domain timeseries toolkit and the frequency-domain spectral family, and it is built the same way: pure Kotlin, no external dependency, works headless in a box.

Every detector here reduces pixels to geometry, and geometry in Field is an FLine, so a Canny result, a set of corners, a Hough line or an optical-flow field lands in the same drawing, animation and export machinery as anything you drew by hand. Scalar responses colour-map through the same Colormaps tables RectSL uses. And Img implements ToHTML, so evaluating one in a box shows the picture, not Img@3f1a2b.

For learned estimates of depth and image content, see depth and CLIP, which runs neural models over the same RGBA bytes an Img can produce.

The value type

An image is a size, a channel count and a flat array of floats:

Img(width, height, channels, data)

width * height * channels floats, interleaved, row-major, origin top-left, the same layout BufferedImage, PNG and a GPU upload all want. Channels are 1 (scalar / gray / mask / gradient / distance), 3 (RGB) or 4 (RGBA); nothing else.

Three choices change what the numbers mean:

• Float, not bytes. Every operation is float math. Operations keep floating-point values through the pipeline instead of quantising to bytes at every stage; clamping happens on output.
• Signed where it should be. Gradients, Laplacians and differences legitimately go negative and are left that way. colormapDiverging(...) is the right way to look at them.
• Gradients are in true units-per-pixel. The Sobel derivative factor is [-1, 0, 1]/2 and the smoothing factor is normalised, so sobel().dx on a ramp of slope s really is s, and a unit step edge integrates to a total gradient of exactly 1. A threshold like 0.05 therefore means the same thing on every gradient operator here.

Build one from a function, a file, a BufferedImage, PNG bytes or a plain FloatArray:

val ramp = Img.generate(64, 64) { x, _ -> x / 63f }
val photo = imgFromFile("cat.png") // javax.imageio under the hood
val fromGl = floatsToImg(heights, w, h) // a heightfield, a spectrogram, any scalar grid

Basics — resample, luma, HSV, threshold

resize area-averages when it shrinks (the anti-aliased answer) and interpolates bilinearly when it grows. luma() is BT.709. otsuThreshold() computes the best two-way split of the histogram rather than making you guess a number.

:::static
// An Img is width x height x channels floats, interleaved, origin top-left.
// Build one straight from a function of (x, y) — the usual way to make a test image.
val card =
    Img.generateRGB(240, 180) { x, y ->
        val u = x / 239f
        val v = y / 179f
        val disc = if (hypot(x - 60.0, y - 60.0) <= 34.0) 1f else 0f
        val bar = if (y in 130..142 && x in 20..219) 1f else 0f
        Triple(
            (0.15f + 0.7f * u * (1f - disc) + disc * 0.9f + bar).coerceIn(0f, 1f),
            (0.15f + 0.6f * v + disc * 0.25f + bar).coerceIn(0f, 1f),
            (0.55f - 0.4f * u + disc * 0.1f + bar).coerceIn(0f, 1f),
        )
    }

val gray = card.luma() // BT.709 perceptual grayscale
val hsv = card.rgbToHsv() // h, s, v each in 0..1
val small = card.resize(60, 45) // area-averaged (anti-aliased) downsample
val backUp = small.resize(240, 180) // bilinear upsample: the detail is really gone
val t = gray.otsuThreshold() // computed, not guessed

showImgGrid(
    "source RGB (${card.width}x${card.height})" to card,
    "luma()" to gray,
    "hue channel" to hsv.channel(0),
    "saturation channel" to hsv.channel(1),
    "resize(60,45) then back up" to backUp,
    "thresholdOtsu() @ ${t.f3}" to gray.threshold(t),
    columns = 3,
    cellSide = 260,
)
println(
    "Img ${card.width}x${card.height} x${card.channels} ch, mean ${card.mean().f3}, " +
        "otsu ${t.f3}, foreground ${gray.threshold(t).data.count { it > 0.5f }} px",
)

showImgGrid(...) is the fastest way to see a whole pipeline: it lays captioned panels out side by side, maps every panel through one shared display range, and pushes one PNG to the box output through Img.toPngBytes() → PngEncoder → display(...). Ordinary images use 0..1; out-of-range panels choose one common automatic range, or you can pass displayRange = -1f..1f (for example). No AWT is involved, so it works anywhere.

Choosing a filter

This box runs four smoothers over the same noisy image and prints each one's residual against the clean source. On impulse noise the median is two orders of magnitude better than the Gaussian. The bilateral filter, which is excellent at preserving edges, is bad at impulses: a salt pixel is far away in range and so keeps its own weight.

:::static
// A clean synthetic source, and the same source with 4% salt-and-pepper corruption.
val clean =
    Img.generate(240, 180) { x, y ->
        var v = 0.12f + 0.10f * (y / 180f)
        if (hypot(x - 60.0, y - 70.0) <= 34.0) v = 0.85f
        if (x in 120..185 && y in 30..95) v = 0.65f
        if (y in 135..142 && x in 20..219) v = 0.95f
        v
    }
val rng = kotlin.random.Random(11)
val noisy = clean.copy()
repeat((noisy.pixelCount * 0.04).toInt()) {
    noisy.data[rng.nextInt(noisy.pixelCount)] = if (rng.nextBoolean()) 0f else 1f
}

val gauss = noisy.gaussianBlur(2.0) // separable, normalised: preserves the image mean
val med = noisy.median(2) // the right tool for impulse noise
val bil = noisy.bilateral(sigmaSpace = 3.0, sigmaRange = 0.12)
val sharp = clean.gaussianBlur(2.0).unsharpMask(sigma = 2.0, amount = 1.4)

showImgGrid(
    "clean source" to clean,
    "+ 4% salt-and-pepper" to noisy,
    "gaussianBlur(2.0)" to gauss,
    "median(2)" to med,
    "bilateral(3.0, 0.12)" to bil,
    "blur then unsharpMask(2.0, 1.4)" to sharp,
    columns = 3,
    cellSide = 260,
)
println(
    "residual vs clean — gaussian ${gauss.meanAbsDiff(clean).f4}, " +
        "median ${med.meanAbsDiff(clean).f4}, " +
        "bilateral ${bil.meanAbsDiff(clean).f4}",
)

boxBlur(radius) is also there, computed with running sums so its cost is independent of the radius, and boxBlurApproxGaussian(sigma) chains three of them when sigma is large and you care more about speed than about the exact kernel shape.

Gradients, edges, and the FLines they become

sobel() returns a Gradient: dx, dy, and derived magnitude(), orientation() and a hue-encoded toColorImg(). canny() runs the four classic stages (smooth, gradient, non-maximum suppression, hysteresis) and returns a thin binary mask. toContourFLines(...) traces that mask and hands back drawable FLines, simplified with Ramer–Douglas–Peucker so a staircase pixel boundary becomes a handful of straight runs.

:::static
val src =
    Img.generate(240, 180) { x, y ->
        var v = 0.12f + 0.08f * (y / 180f)
        if (hypot(x - 62.0, y - 70.0) <= 34.0) v = 0.85f
        if (x in 120..185 && y in 30..95) v = 0.68f
        val tx = x - 210.0
        val ty = y - 40.0
        if (ty in 0.0..60.0 && kotlin.math.abs(tx) <= ty * 0.5) v = 0.95f
        v
    }

// Sobel here is scaled to true units-per-pixel: a unit step edge integrates to 1.
val grad = src.sobel()
val edges = src.canny(low = 0.05, high = 0.11, sigma = 1.2)
val auto = src.cannyAuto() // thresholds derived from this image's own gradients

// Edges are geometry, and geometry in Field is an FLine.
val flines = edges.toContourFLines(minPerimeter = 12, simplifyTolerance = 1.2)

showImgGrid(
    "source" to src,
    "sobel().dx (signed, coolwarm)" to grad.dx.colormapDiverging("coolwarm"),
    "sobel().dy (signed, coolwarm)" to grad.dy.colormapDiverging("coolwarm"),
    "magnitude (inferno)" to grad.magnitude().normalized().colormap("inferno"),
    "orientation as hue x magnitude" to grad.toColorImg(),
    "canny(0.05, 0.11)" to edges,
    columns = 3,
    cellSide = 260,
)
println(
    "|grad| max ${grad.magnitude().max().f3}; " +
        "canny ${edges.data.count { it > 0.5f }} px, cannyAuto ${auto.data.count { it > 0.5f }} px; " +
        "${flines.size} FLines, ${flines.sumOf { it.nodes.size }} nodes total",
)

dx fires on vertical edges and dy on horizontal ones; the orientation panel wheels its hue right around the disc, because that is what the gradient direction does.

Canny on low-contrast images

A (low, high) pair tuned on a bright synthetic card finds almost nothing in a dark, softly-lit photograph. The thresholds are in gradient-magnitude units, and a dark image does not reach them. cannyAuto(...) derives them from the image's own distribution: it takes a percentile of the non-maximum-suppressed magnitudes as the strong seed threshold and Canny's own recommended ratio for the weak one.

The figure below is a real repo PNG run through the whole pipeline. The fixed pair finds 252 edge pixels, barely the car at all. cannyAuto() finds 1952, and traces the whole silhouette; the last panel is those edges as 57 FLines drawn back over the source.

Corners

Three detectors, all on a square whose corners are at exactly known coordinates, so each detector's answer can be measured against the truth.

• harrisResponse() — \( \det(M) - k\,\operatorname{trace}(M)^2 \) over the structure tensor. Large and positive at a corner, negative along an edge, near zero on flat ground; the diverging colour map in the figure is showing exactly that sign structure.
• shiTomasiResponse() — the smaller eigenvalue of the same tensor. No magic k, and directly interpretable: it is the weakest intensity variation over any direction, so a large value means the point is localisable in both axes. A straight edge scores ~0, which is the aperture problem stated numerically. This is the response an optical-flow tracker wants, and goodFeaturesToTrack(...) is its friendlier name.
• fastCorners() — no gradients and no eigenvalues, just a wrap-around run-length test on a 16-pixel circle. Much faster, and it scores each corner by the largest threshold at which it would still qualify.

The two response maps are drawn from a 240 px copy of the square, shown 1:1, and raised to the 0.25 power so the weak values show.

Run headless, Harris and Shi–Tomasi each return exactly four corners, and the worst of the four sits 0.94 px from the true corner of the square (sub-pixel refinement is on by default, via a three-tap parabola fit in x and y). cornersToFLines(...) turns them into drawable crosses or boxes.

Hough lines

Every edge pixel votes for every line through it; a real line becomes a bright peak in (rho, theta) accumulator space, with no grouping step in between.

The figure has three planted lines and a decoy ring. All three lines come back at their exact polar coordinates: (rho=1500, θ=0.00), (rho=400, θ=1.57), (rho=283, θ=2.36). The ring is not fitted.

HoughLine.clipTo(w, h) clips the infinite line to the image rectangle, and linesToFLines(...) gives you the visible segments as FLines. houghCircles(...) is there too, for a known radius range.

Morphology, blobs and distance

erode / dilate / open / close with square, disc or cross structuring elements; connectedComponents(...) for labels plus per-blob area, bounding box, centroid and fill ratio; distanceTransform() for the exact Euclidean distance to the nearest background pixel (the Felzenszwalb–Huttenlocher two-pass algorithm, not a chamfer approximation).

:::static
// Three disjoint discs of known radius, then noise on top of them.
val discs = [Triple(55.0, 60.0, 28.0), Triple(150.0, 55.0, 20.0), Triple(110.0, 135.0, 14.0)]
val clean =
    Img.generate(220, 180) { x, y ->
        if (discs.any { (cx, cy, r) -> hypot(x - cx, y - cy) <= r }) 0.85f else 0.15f
    }
val rng = kotlin.random.Random(3)
val noisy = clean.copy()
repeat((noisy.pixelCount * 0.03).toInt()) {
    noisy.data[rng.nextInt(noisy.pixelCount)] = if (rng.nextBoolean()) 0f else 1f
}

val mask = noisy.thresholdOtsu()
val cleaned = mask.open(2, StructuringElement.DISC).close(2, StructuringElement.DISC)
val cc = cleaned.connectedComponents(8)
val dt = cleaned.distanceTransform()

showImgGrid(
    "noisy source" to noisy,
    "thresholdOtsu()" to mask,
    "open(2,DISC).close(2,DISC)" to cleaned,
    "erode(4)" to cleaned.erode(4, StructuringElement.DISC),
    "dilate(4)" to cleaned.dilate(4, StructuringElement.DISC),
    "distanceTransform (magma)" to dt.normalized().colormap("magma"),
    columns = 3,
    cellSide = 250,
)
println("${cc.count} components after cleaning (planted ${discs.size} discs):")
cc.components
    .sortedByDescending { it.area }
    .forEach {
        println(
            "  area=${it.area} centroid=(${it.centroidX.f1}, ${it.centroidY.f1}) " +
                "bbox=${it.boundingWidth}x${it.boundingHeight} fill=${it.fillRatio.f2}",
        )
    }
println("distance-transform peak ${dt.max().f1} (largest planted radius ${discs.maxOf { it.third }})")

Three discs in, three components out, with areas of 2419, 1226 and 584 against the analytic \( \pi r^2 \) of 2463, 1257 and 616 (the shortfall is the open(2) erosion), and centroids at (54.8, 60.1), (150.1, 55.0) and (110.2, 134.9) against planted centres of (55, 60), (150, 55) and (110, 135).

For sub-pixel-smooth contours rather than a pixel-staircase boundary, isoContourFLines(iso) adapts the image to the ScalarField2D that Field's existing marching-squares extractor (contours) already marches, so a distance transform or a blurred response becomes a smooth curve with no glue code.

Multiscale

gaussianPyramid(levels) halves each octave (blurring first, so it does not alias); laplacianPyramid(levels) splits the image into band-pass detail plus a coarse residual, and reconstruct() inverts it. The round-trip is exact to within float precision, at a maximum absolute error of 5.96e-08 on the figure below. That is small enough to build multiband blending and detail transfer on, rather than only to look at.

(Levels below the base are magnified by integer pixel replication for display; each caption states the real size and the zoom factor.)

Geometry — affine and projective warps

Affine2D uses the same six-number layout as java.awt.geom.AffineTransform and SVG's matrix(a, b, c, d, e, f). Homography is the 3×3 projective version. warpAffine and warpPerspective map source → destination and invert internally, doing the correct backwards resample, so there are no holes.

findHomography(src, dst) is the Direct Linear Transform: exact through four correspondences, least-squares through more. The points are Hartley-normalised first; without that, the conditioning on pixel-scale coordinates is bad enough to visibly distort the fit. warpQuadToRect(quad, w, h) wraps it up as the "flatten the photographed page" operation.

The rectified checkerboard matches the original to a mean absolute error of 0.0036 over the interior.

Motion

opticalFlowLK(a, b, points) is pyramidal Lucas–Kanade. Seed it with goodFeaturesToTrack(...), which selects exactly the points where the 2×2 tracking system is well conditioned, and it returns a FlowField of FlowVectors, each of which reports whether it converged. A track is marked ok = false when the system is ill-conditioned (a smooth region, or a straight edge) or the point leaves the frame.

The pyramid is there for large motions. Solve at the coarsest octave, where the displacement is small enough for the linearisation to hold, then propagate the estimate down. A 9-pixel shift that a single level cannot recover comes back to within half a pixel with four.

:::static
// Two frames related by a transform we choose, so the answer is known in advance.
val a =
    Img.generate(200, 200) { x, y ->
        (0.5 + 0.24 * kotlin.math.sin(x * 0.21) * kotlin.math.cos(y * 0.17) + 0.16 * kotlin.math.sin((x + y) * 0.11))
            .toFloat()
            .coerceIn(0f, 1f)
    }
val planted = Affine2D.rotationAbout(0.05, 100.0, 100.0).compose(Affine2D.translation(5.0, -3.0))
val b = a.warpAffine(planted, border = Border.CLAMP)

// Seed only where the 2x2 tracking system is well conditioned.
val seeds =
    a.goodFeaturesToTrack(maxCount = 200, minDistance = 9)
        .filter { it.x in 24.0..176.0 && it.y in 24.0..176.0 }
        .map { it.toVec2() }
val flow = opticalFlowLK(a, b, seeds, windowRadius = 8, levels = 3)
val fit = findHomographyRansac(
    flow.sourcePoints(),
    flow.targetPoints(),
    threshold = 1.5,
    iterations = 1200,
    seed = 3L,
)

// How far is the recovered model from the one we planted, in pixels, on a probe grid?
var worst = 0.0
for (gx in 40..160 step 40) for (gy in 40..160 step 40) {
    val (fx, fy) = fit.homography.apply(gx.toDouble(), gy.toDouble())
    val (tx, ty) = planted.apply(gx.toDouble(), gy.toDouble())
    worst = maxOf(worst, hypot(fx - tx, fy - ty))
}

showImgGrid(
    "frame A" to a,
    "frame B (rot 0.05 rad + shift 5,-3)" to b,
    "|A - B| (inferno)" to frameDifference(a, b).normalized().colormap("inferno"),
    "B warped back by the fitted model" to b.warpPerspective(fit.homography.inverse(), border = Border.CLAMP),
    columns = 2,
    cellSide = 280,
)
val (mdx, mdy) = flow.medianDisplacement()
println(
    "${flow.tracked().size}/${seeds.size} tracked; median displacement " +
        "(${mdx.f2}, ${mdy.f2}) px",
)
println("RANSAC ${fit.inlierCount} inliers of ${flow.tracked().size} after ${fit.iterations} iterations")
println("recovered model agrees with the planted transform to ${worst.f3} px over the probe grid")
println(fit.homography.toString())

With findHomographyRansac(...) half the correspondences can be garbage and the answer is still right, because a wrong model has no reason to be consistently wrong across many points. It samples four correspondences at a time, scores the consensus, keeps the best, and refits least-squares on the inliers. The seed is explicit, so a docs capture or a test gets the same answer every time.

opticalFlowBlockMatch(...) is the brute-force alternative for displacements large enough to break Lucas–Kanade's linearisation, and frameDifference(...) / updateBackground(...) are the crude, reliable motion cues; pair them with open(...) to drop sensor noise and connectedComponents(...) to count what moved.

FlowField.toFLines(scale, ...) draws the whole field as arrows, shaft plus arrowhead, so direction is visible and not just magnitude:

Colour

luma(), rgbToHsv() / hsvToRgb(), and toLinear() / toSrgb(). Blurring, downsampling and averaging are only correct in linear light, so img.toLinear().gaussianBlur(4.0).toSrgb() is measurably different from blurring the encoded values. The panel below shows both, and the amplified difference between them.

colormap(name, lo, hi) runs a 1-channel image through any of the Colormaps ramps (viridis, magma, inferno, plasma, turbo, coolwarm, grayscale); colormapDiverging(name) picks a symmetric range around zero, the right default for anything signed.

Trying it yourself

The imaging_cv_examples starting point builds the examples on this page as live boxes. Change a sigma, a threshold or the planted motion, re-run the box, and watch the picture update.

Speed, against OpenCV

These filters are pure Kotlin with no native dependency. The repository carries a benchmark that puts them beside OpenCV in the same JVM, on the same pixels: tools/imaging-bench, run with

./gradlew :tools:imaging-bench:jmh

The table below compares the two in milliseconds per call.

Bold means we are faster. A dash in the "before" column means that operation has no preserved scalar original to compare against, not that its speedup is unknown.

Gaussian blur, box blur and the Laplacian run the same operation with the same parameters on both sides, and the benchmark fails if their outputs stop agreeing. The rest are matched in work but not in definition, so their ratios compare slightly different operations: our luma is BT.709 where cvtColor is BT.601; our Sobel and Scharr are scaled to true units-per-pixel, so ours is OpenCV's divided by a constant; our pyrDown takes a configurable sigma where OpenCV's is a fixed 5-tap kernel; and median/bilateral have matched parameters but different internal algorithms. Our API also returns a fresh Img each call, so part of our cost is an allocation OpenCV is not paying.

Gaussian blur, dense convolution and the RGB blur are at or slightly past OpenCV; box blur, the Laplacian, Scharr, Sobel, luma and bilateral are all inside about 2x, close enough that the pure-Kotlin path is the reasonable default. bilateral's interval is wide, so read its row as "about 1.2x". Three rows are meaningfully behind, and you should know which before choosing:

• median is an algorithmic gap: it sorts 25 floats per pixel where OpenCV uses a sorting network.
• resize to half has data-dependent inner loops that do not vectorize.
• pyrDown is doing more arithmetic, not doing it worse: ours applies gaussianKernel1D(sigma), 9 taps at the default, where OpenCV's pyrDown is a fixed 5-tap [1,4,6,4,1]/16. That is close to twice the work, in exchange for a better filter and a sigma that is yours to choose.

The benchmark's conditions, the agreement tolerances, and what the absolute milliseconds depend on are in dev_docs/IMAGING_MEASUREMENTS.md.

Reference

Source: imaging.kt, imaging_filters.kt, imaging_features.kt, imaging_geometry.kt, imaging_field.kt in src/jvmMain/kotlin/afield/media/image/.