Resources Numerical accuracy Calculation policy
Works out the scale, rotation and position that put a scanned drawing into project coordinates, from a stated scale or from matched points, and converts traced lines with it.
64-bit numbers and closed-form least squares. Unit factors are exact definitions. The similarity fit centers both the pixel and the project coordinates before summing; the affine fit centers the pixel coordinates only. Residuals are reported in project units for every point. Traced arcs become chords of at most 5°.
Open item. Center the project coordinates as well as the pixel coordinates in the affine fit, as the similarity fit already does.
12 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
Tests bin/test-scanstart.mjs
The length of each supported unit in meters.
ft = 0.3048; US survey ft = 1200 ÷ 3937; in = 0.0254; yd = 0.9144; mi = 1609.344
Each factor is the unit's exact definition; the survey foot is written as the fraction 1200 ÷ 3937.
export const UNIT_M = { m: 1, ft: 0.3048, ftUS: 1200 / 3937, mm: 0.001, cm: 0.01, in: 0.0254, km: 1000, mi: 1609.344, yd: 0.9144 };Project units per pixel from two points on the scan and the true distance between them.
units per pixel = distance ÷ √(Δpx² + Δpy²)
One Math.hypot and one division. A zero pixel distance or a distance that is not positive returns null.
export function scaleFromTwo(p1, p2, dist) {
const d = Math.hypot(p2.x - p1.x, p2.y - p1.y);
if (!(d > 0) || !fin(dist) || !(dist > 0)) return null;
return dist / d;
}Project units per pixel from the scan's resolution and a scale such as 1 inch = 50 feet or 1:500.
per ÷ dpi, or ratio × 0.0254 ÷ (m per unit) ÷ dpi
Uses the exact inch and the unit table. The result is only as good as the dpi recorded for the scan.
/** At `dpi`, either 1 inch on paper = `per` project units (an engineer's scale, 1" = 50'), or 1 : `ratio` (paper → real) with the
* project in `unit`. */
export function scaleFromStated({ dpi, per = null, ratio = null, unit = 'ft' } = {}) {
if (!fin(dpi) || !(dpi > 0)) return null;
if (fin(per) && per > 0) return per / dpi;
if (fin(ratio) && ratio > 0 && UNIT_M[unit]) return (ratio * UNIT_M.in / UNIT_M[unit]) / dpi;
return null;
}The transform that puts the sheet's lower-left corner at a point, turned by an angle, at a given scale.
c = s·cos θ, d = s·sin θ; X = x₀ − d·h + c·px + d·py; Y = y₀ + c·h + d·px − c·py
The angle in degrees is multiplied by Math.PI ÷ 180. The pixel row count h enters the offsets so that pixel (0, h) lands on the origin.
/** The sheet (h pixels tall) with its lower-left corner at origin, turned anticlockwise by rotationDeg, upp units per pixel. */
export function gtFromScale({ upp, h, origin = { x: 0, y: 0 }, rotationDeg = 0 }) {
const th = rotationDeg * Math.PI / 180, c = upp * Math.cos(th), d = upp * Math.sin(th);
// X = c·px + d·py + tx, Y = d·px − c·py + ty (y down on the sheet, up in the project); pixel (0, h) → origin
return [origin.x - d * h, c, d, origin.y + c * h, d, -c];
}The project position of a pixel.
X = g₀ + g₁·px + g₂·py; Y = g₃ + g₄·px + g₅·py
Two products added to an offset per coordinate. With large grid coordinates the offset dominates, so the sum carries the offset's absolute precision.
export const apply = (gt, px, py) => [gt[0] + gt[1] * px + gt[2] * py, gt[3] + gt[4] * px + gt[5] * py];The transform from project coordinates back to pixels.
det = g₁·g₅ − g₂·g₄; inverse = adjugate ÷ det
The determinant is tested only against exact zero, so a nearly singular placement is inverted with whatever digits the division leaves.
export function invert(gt) {
const det = gt[1] * gt[5] - gt[2] * gt[4];
if (!det) return null;
const a = gt[5] / det, b = -gt[2] / det, d = -gt[4] / det, e = gt[1] / det;
return [-(a * gt[0] + b * gt[3]), a, b, -(d * gt[0] + e * gt[3]), d, e];
}The project units per pixel a placement stands for, reported with every fit.
scale = √|g₁·g₅ − g₂·g₄|
The square root of the absolute determinant, which is the mean scale for an affine placement that stretches the two directions differently.
export const unitsPerPixel = (gt) => Math.sqrt(Math.abs(gt[1] * gt[5] - gt[2] * gt[4]));The angle a placement turns the sheet by, in degrees.
rotation = atan2(g₄, g₁) × 180 ÷ π
One Math.atan2 of two transform terms; for an affine placement with shear it is the direction of the pixel x axis only.
export const rotationOf = (gt) => Math.atan2(gt[4], gt[1]) * 180 / Math.PI;Solves the normal equations of the affine fit.
Gauss-Jordan elimination with partial pivoting
The largest entry in each column is taken as the pivot. A pivot below 1e-12 in absolute value returns null; the threshold is not scaled to the size of the sums, which grow with the square of the pixel distances.
function solve3(A, b) { // Gaussian elimination with pivoting; A 3×3, b 3 → x | null
const M = A.map((r, i) => [...r, b[i]]);
for (let c = 0; c < 3; c++) {
let p = c; for (let r = c + 1; r < 3; r++) if (Math.abs(M[r][c]) > Math.abs(M[p][c])) p = r;
if (Math.abs(M[p][c]) < 1e-12) return null;
[M[c], M[p]] = [M[p], M[c]];
for (let r = 0; r < 3; r++) if (r !== c) { const f = M[r][c] / M[c][c]; for (let k = c; k < 4; k++) M[r][k] -= f * M[c][k]; }
}
return M.map((r, i) => r[3] / M[i][i]);
}The placement that best matches points picked on the scan to their project coordinates, with each point's residual and the root mean square.
similarity: c = Σ(u·X + v·Y) ÷ Σ(u² + v²), d = Σ(u·Y − v·X) ÷ Σ(u² + v²); rms = √(Σ r² ÷ n)
The similarity fit subtracts the means of both the pixel and the project coordinates before summing, so large grid coordinates do not enter the sums. The affine fit centers the pixel coordinates but multiplies them by the uncentered project coordinates, so its right-hand sums are made of large terms that mostly cancel.
export function fitPoints(pairs, { mode = 'similarity' } = {}) {
const P = (pairs || []).filter((q) => q && [q.px, q.py, q.x, q.y].every(fin));
let gt = null;
if (mode === 'affine') {
if (P.length < 3) return null;
// X = x0 + a·px + b·py; Y = y0 + d·px + e·py — normal equations, centred for conditioning
const n = P.length, mx = P.reduce((s, q) => s + q.px, 0) / n, my = P.reduce((s, q) => s + q.py, 0) / n;
const A = [[0, 0, 0], [0, 0, 0], [0, 0, 0]], bx = [0, 0, 0], by = [0, 0, 0];
for (const q of P) {
const r = [1, q.px - mx, q.py - my];
for (let i = 0; i < 3; i++) { for (let j = 0; j < 3; j++) A[i][j] += r[i] * r[j]; bx[i] += r[i] * q.x; by[i] += r[i] * q.y; }
}
const sx = solve3(A, bx), sy = solve3(A, by);
if (!sx || !sy) return null; // the points lie on one line
gt = [sx[0] - sx[1] * mx - sx[2] * my, sx[1], sx[2], sy[0] - sy[1] * mx - sy[2] * my, sy[1], sy[2]];
} else {
if (P.length < 2) return null;
// u = px, v = −py (y up); X = c·u − d·v + tx, Y = d·u + c·v + ty
const n = P.length;
const mu = P.reduce((s, q) => s + q.px, 0) / n, mv = P.reduce((s, q) => s - q.py, 0) / n;
const mX = P.reduce((s, q) => s + q.x, 0) / n, mY = P.reduce((s, q) => s + q.y, 0) / n;
let sxx = 0, sc = 0, sd = 0;
for (const q of P) {
const u = q.px - mu, v = -q.py - mv, X = q.x - mX, Y = q.y - mY;
sxx += u * u + v * v; sc += u * X + v * Y; sd += u * Y - v * X;
}
if (!(sxx > 0)) return null; // every scan point in one place
const c = sc / sxx, d = sd / sxx;
const tx = mX - (c * mu - d * mv), ty = mY - (d * mu + c * mv);
gt = [tx, c, d, ty, d, -c];
}
const residuals = P.map((q) => { const [X, Y] = apply(gt, q.px, q.py); return Math.hypot(X - q.x, Y - q.y); });
const rms = Math.sqrt(residuals.reduce((s, r) => s + r * r, 0) / residuals.length);
// with exactly the minimum number of points every residual is 0: nothing checks the points (hosts say so); one more point shows
// THAT they disagree but not which one is wrong (any one of them could be); from two more, the furthest off is named
const min = mode === 'affine' ? 3 : 2, checked = P.length > min;
const worst = P.length >= min + 2 ? residuals.reduce((w, r, i) => (r > residuals[w] ? i : w), 0) : null;
return { gt, mode, scale: unitsPerPixel(gt), rotationDeg: rotationOf(gt), residuals, rms, worst, checked };
}Turns an arc found on the scan into points for a line layer.
n = max(2, ⌈|sweep| ÷ 5°⌉); point = center + r·(cos t, sin t)
The arc is stored as chords of at most 5°, so the saved line lies inside the arc by up to r × (1 − cos 2.5°), about 0.1 percent of the radius.
/** Arc (cc-raster trace(): y-down pixel angles, a0 < a1) → points every ≤ 5°. */
function arcPoints(a) {
const span = a.a1 - a.a0, n = Math.max(2, Math.ceil(Math.abs(span) / (Math.PI / 36)));
const out = [];
for (let i = 0; i <= n; i++) { const t = a.a0 + span * i / n; out.push({ x: a.cx + a.r * Math.cos(t), y: a.cy + a.r * Math.sin(t) }); }
return out;
}The lines and arcs traced on the scan as one line layer in project coordinates.
each end point through the placement; kept when its length > the minimum
Every point goes through the placement once. A line is dropped when its length in project units is at or below the minimum length, 0 by default. The pixel coordinates come from the tracing step and carry its resolution.
export function traceLayers(result, gt, { name = 'SCAN-TRACE', source = null, minLen = 0 } = {}) {
const features = [];
const pt = (p) => apply(gt, p.x, p.y);
for (const l of result?.lines || []) {
const a = pt({ x: l.x1, y: l.y1 }), b = pt({ x: l.x2, y: l.y2 });
if (![...a, ...b].every(fin) || Math.hypot(b[0] - a[0], b[1] - a[1]) <= minLen) continue;
features.push({ g: { t: 'ln', c: [[a, b]] }, p: { KIND: 'line' } });
}
for (const r of result?.arcs || []) {
const c = arcPoints(r).map(pt);
if (c.length < 2 || !c.flat().every(fin)) continue;
features.push({ g: { t: 'ln', c: [c] }, p: { KIND: 'arc' } });
}
if (!features.length) return [];
return [{ name, geom: 'ln', hasZ: false, fields: [{ name: 'KIND', type: 'text' }], features, source }];
}Every value is a 64-bit number from input to result. Formulas are closed form or a direct solve, and every comparison uses a stated tolerance.