Fitting

This page covers two helpers. One turns a few numbers into a whole picture, and the other works back to a few numbers from a measure of how wrong they are.

RadialBasis goes outward. You know a value at a handful of scattered places, and it gives you a value everywhere. Fit.minimize goes inward. You have a handful of parameters, and a way to measure how wrong a setting of them is. You want the setting that is least wrong.

Contents

• A field through scattered values
• Choosing a kernel
• Fields of vectors, and warps
• The slope of a field
• Smoothing
• Finding a handful of numbers
• What each one costs

A field through scattered values

You can interpolate a grid by looking at the cells around a spot. Scattered points have no such neighbors, so the answer has to come from all of them at once. A radial basis function does this by giving every known point a bump centered on it. It then weights the bumps so that their sum passes exactly through every value you gave:

let field = RadialBasis(points: [Vector2(120, 140), Vector2(700, 300), Vector2(400, 860)],
                        values: [0.1, 0.9, 0.4])!

for y in stride(from: 0.0, to: height, by: 10) {
    for x in stride(from: 0.0, to: width, by: 10) {
        fill(Color(white: field.value(at: Vector2(x + 5, y + 5))))
        drawRect(x, y, 11, 11)
    }
}

The values can be Doubles, Vector2s, Vector3s, or Colors, and the points can be Vector2 or Vector3. All the channels of a value share one solve, so a field of colors costs the same as a field of numbers.

A type of your own joins them by conforming: a RadialBasisPoint says how many axes it has (axisCount), hands one over by index (component(_:)), and measures a distance(to:) another; a RadialBasisValue says how many numbers it carries (channelCount), hands them over as channels, and builds itself back from them with init(channels:).

The initializer returns nil when the points cannot define a field. That happens when there are no points, or when the two lists have different lengths. It also happens when two points in the same place carry different values, or when all the points sit on one straight line.

Choosing a kernel

The kernel is the shape of each bump. Six are available, and the default is a good first choice for almost everything:

The first three grow without limit away from their point. That sounds wrong, but it is not, because only the sum matters. A flat plane is fitted alongside the bumps, and it cancels the growth. So those three extrapolate sensibly past the edge of the data. The last three fade instead, so the field settles far from every point rather than running off. Each of those three takes a scale in the same units as your points, and a good value is about the spacing between them.

Fields of vectors, and warps

A field of Vector2s is a warp. You pin a few places to where they should move, and everything between them follows smoothly:

let warp = RadialBasis(points: pins, values: targets)!
let bent = straightLine.map { warp.value(at: $0) }

When you read every point of a shape through the warp, the shape bends. When you read a grid through it, you get the classic rubber-sheet picture. Each pinned place lands exactly on its target, because the field interpolates the values you gave.

The slope of a field

A fitted field can tell you how fast it changes, not only what it is. For a field of numbers, gradient(at:) returns the direction the field climbs fastest, as a vector whose length is how steep the climb is. The perpendicular to it runs along the field's contour, so a particle can climb a fitted surface or circle on it:

let height = RadialBasis(points: [Vector2(200, 300), Vector2(800, 260), Vector2(500, 800)],
                         values: [0.9, 0.2, 0.6])!
var walker = Vector2(540, 540)
let uphill = height.gradient(at: walker)
walker += uphill.perpendicular.normalized * 2    // one step along the contour

derivative(at:along:) works on every field, whatever its values are. It says how the value changes as the point moves along a direction. A unit direction gives the rate per unit of distance, and a longer one scales the answer with its length. For a warp, the derivatives along the two axes are the warp's Jacobian. It says how a small square at that point is stretched and turned:

let warp = RadialBasis(points: [Vector2(300, 300), Vector2(780, 320), Vector2(540, 800)],
                       values: [Vector2(340, 280), Vector2(760, 360), Vector2(540, 760)])!
let p = Vector2(540, 540)
let dx = warp.derivative(at: p, along: Vector2(1, 0))
let dy = warp.derivative(at: p, along: Vector2(0, 1))
let stretch = dx.x * dy.y - dx.y * dy.x          // area grows by this factor here

Both are exact, because each kernel has a closed-form derivative, and a read costs about the same as reading the value. At a known point, the linear kernel's cone has no single slope, and the answer there leaves that point's own bump out. For a function you wrote yourself rather than fitted, the numerical helpers measure the same slope by sampling.

Smoothing

At smoothing: 0 the field passes through every value exactly. That is right for clean data, but it makes noisy data ring. When you raise the value, the field is allowed to miss the points, and the higher the value, the more it may miss. In exchange, there are fewer wobbles between the points:

RadialBasis(points: samples, values: readings, smoothing: 0.05)

The value is measured against the scale of the fit, not in raw units. So the same number means the same thing whether your points are spread over a thousand pixels or over one. Around 0.01 is light smoothing, and 1 is heavy.

Finding a handful of numbers

Fit.minimize goes the other direction. You give it a starting setting and a cost function, and it moves the setting downhill, toward a lower cost:

// The circle that passes closest to a set of marks.
let best = Fit.minimize(from: [width / 2, height / 2, 100]) { p in
    marks.reduce(0.0) { total, mark in
        let off = Vector2(p[0], p[1]).distance(to: mark) - p[2]
        return total + off * off
    }
}
drawCircle(best.values[0], best.values[1], best.values[2])

The result carries the best setting found, its cost, and how many steps it took. It also tells you whether it isSettled on its own rather than running out of steps. A walk that used every step may need more steps, or a larger rate.

Two things are worth knowing. First, each parameter moves by about rate per step, whatever the slope is at that point. So one rate works for a parameter measured in pixels beside a parameter measured in turns. Second, the walk goes downhill from where you start. So when a cost has several separate valleys, the walk settles in whichever valley your starting point sat in. When that matters, start the walk from several places and keep the best answer.

What each one costs

Fitting a field solves a dense system, so its cost grows with about the cube of the number of points. A few hundred points fit in milliseconds. For tens of thousands of points, this is the wrong tool. Reading the field afterward costs one term per point, on every read. So for a large fit read over a whole canvas, the reading is the expensive half. Fit once in setup(), and read in draw(). A slope read (gradient(at:) or derivative(at:along:)) costs about the same as a value read.

Minimizing measures the slope rather than deriving it. It does this by trying each parameter a little to either side of where it stands. So your cost function is called about twice per parameter per step. Keep the cost function cheap. Remember that 300 steps over 3 parameters is already 1,800 calls.

See also Geometry for the value types both helpers work with, Spatial queries for finding which scattered points are near a place, and Flow fields for a field generated from noise rather than fitted through values you chose.