Reads a scanned image, turns it into black and white, finds how far the sheet is rotated and turns it straight.
How the precision is kept
Pixel values are 8-bit by the nature of the image; the arithmetic that decides them is 64-bit (sums in Float64Array, 64-bit cosine table). The skew angle is found by trying angles on a fixed grid and is reported to 0.01°.
Open item. Round the converted red, green and blue values before they are stored; they are currently cut to whole numbers by the byte array.
Calculations under this policy
6 calculations in 2 source files. 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.
The grey level that best separates ink from paper for the whole sheet.
maximize w_B × w_F × (m_B − m_F)² over the 256 levels
The histogram is a Float64Array of whole-number counts, which are exact up to 2⁵³. The level returned is an integer from 0 to 255.
raster/image.js
lines 39 to 56
CadCraft 0.8.2
export function otsu(grey) {
const hist = new Float64Array(256), d = grey.data;
for (let i = 0; i < d.length; i++) hist[d[i]]++;
const total = d.length;
let sum = 0;
for (let t = 0; t < 256; t++) sum += t * hist[t];
let wB = 0, sumB = 0, best = 0, level = 127;
for (let t = 0; t < 256; t++) {
wB += hist[t];
if (!wB) continue;
const wF = total - wB;
if (!wF) break;
sumB += t * hist[t];
const mB = sumB / wB, mF = (sum - sumB) / wF, between = wB * wF * (mB - mF) * (mB - mF);
if (between > best) { best = between; level = t; }
}
return level; // ink: value <= level
}
Local threshold (Sauvola)
64-bit
Decides ink or paper for each pixel from the mean and spread of the pixels around it.
ink when v ≤ mean × (1 + k × (sd ÷ 128 − 1)); sd = √(Σv² ÷ n − mean²)
Window sums come from two integral images in Float64Array, which hold exact whole numbers. The variance is a difference of two nearly equal terms on flat paper and is floored at 0 before the square root.
raster/image.js
lines 86 to 114
CadCraft 0.8.2
export function threshold(grey, { mode = 'auto', level = 128, window = 0, k = 0.2 } = {}, ws = DEFAULT_WS) {
const { w, h, data } = grey;
checkSize(w, h);
const out = new Uint8Array(w * h);
if (mode === 'level') { for (let i = 0; i < out.length; i++) out[i] = data[i] <= level ? 1 : 0; return { w, h, data: out }; }
if (mode === 'auto') { const f = flatten(grey, { window }, ws), t = otsu(f), o = new Uint8Array(w * h); for (let i = 0; i < o.length; i++) o[i] = f.data[i] <= t ? 1 : 0; return { w, h, data: o }; }
if (mode === 'otsu') { const t = otsu(grey); for (let i = 0; i < out.length; i++) out[i] = data[i] <= t ? 1 : 0; return { w, h, data: out }; }
// Sauvola: local mean and deviation from integral images; the window is odd, about 1/16 of the short side by default
const r = Math.max(7, (window > 0 ? window : Math.round(Math.min(w, h) / 16)) >> 1);
const W1 = w + 1;
const I = ws.f64('thrI', W1 * (h + 1)), I2 = ws.f64('thrI2', W1 * (h + 1));
I.fill(0, 0, W1); I2.fill(0, 0, W1);
for (let y = 0; y < h; y++) {
let s = 0, s2 = 0;
I[(y + 1) * W1] = 0; I2[(y + 1) * W1] = 0;
for (let x = 0; x < w; x++) { const v = data[y * w + x]; s += v; s2 += v * v; I[(y + 1) * W1 + x + 1] = I[y * W1 + x + 1] + s; I2[(y + 1) * W1 + x + 1] = I2[y * W1 + x + 1] + s2; }
}
for (let y = 0; y < h; y++) {
const y0 = Math.max(0, y - r), y1 = Math.min(h, y + r + 1);
for (let x = 0; x < w; x++) {
const x0 = Math.max(0, x - r), x1 = Math.min(w, x + r + 1), n = (x1 - x0) * (y1 - y0);
const sum = I[y1 * W1 + x1] - I[y0 * W1 + x1] - I[y1 * W1 + x0] + I[y0 * W1 + x0];
const sum2 = I2[y1 * W1 + x1] - I2[y0 * W1 + x1] - I2[y1 * W1 + x0] + I2[y0 * W1 + x0];
const mean = sum / n, sd = Math.sqrt(Math.max(0, sum2 / n - mean * mean));
out[y * w + x] = data[y * w + x] <= mean * (1 + k * (sd / 128 - 1)) ? 1 : 0;
}
}
return { w, h, data: out };
}
Sheet skew angle
64-bit
The angle by which the long lines of a scanned sheet are rotated.
score(angle) = Σ rows² + Σ columns² of the ink projected at that angle; best angle on a 0.25° then 0.02° grid
The projected coordinate is truncated to a whole-pixel bin with a bitwise OR ((y + diag) | 0); that is the binning itself and the offset keeps the value non-negative. The answer is the best grid angle, so its resolution is 0.02°, and it is then rounded to 0.01° with Math.round before it is returned. The fine pass adds 0.02 repeatedly, which is not exact in binary.
raster/image.js
lines 163 to 199
CadCraft 0.8.2
/**
* The sheet's skew from its long lines: for each trial angle the ink points are projected across the lines' direction and the
* rows and columns of that projection are counted; long straight lines make a few very full rows/columns, so the angle where the
* sum of squares peaks is the skew. A coarse pass (0.25°) over ±maxDeg, then 0.02° around the best. Points are sampled (at most
* 150 000) from the ink. `confidence` = peak over the median of the trial scores (≈ 1 for a sheet with no lines).
*/
export function findSkew(bm, { maxDeg = 10 } = {}, ws = DEFAULT_WS) {
const { w, h, data } = bm;
let ink = 0;
for (let i = 0; i < data.length; i++) ink += data[i];
if (ink < 50) return { deg: 0, confidence: 0, points: ink };
const step = Math.max(1, Math.ceil(ink / 150000)), px = ws.f64('skx', Math.ceil(ink / step) + 1), py = ws.f64('sky', Math.ceil(ink / step) + 1);
let n = 0, seen = 0;
const cx = w / 2, cy = h / 2;
for (let y = 0; y < h; y++) for (let x = 0; x < w; x++) if (data[y * w + x] && (seen++ % step === 0)) { px[n] = x - cx; py[n] = cy - y; n++; } // y up
const diag = Math.ceil(Math.hypot(w, h)) + 4, hist = ws.f64('skh', 2 * diag + 2), hist2 = ws.f64('skh2', 2 * diag + 2);
const score = (deg) => {
const a = deg * Math.PI / 180, c = Math.cos(a), s = Math.sin(a);
hist.fill(0, 0, 2 * diag + 2); hist2.fill(0, 0, 2 * diag + 2);
for (let i = 0; i < n; i++) {
// rotate the points by −deg, so lines at +deg become level
const xr = px[i] * c + py[i] * s, yr = -px[i] * s + py[i] * c;
hist[(yr + diag) | 0]++; hist2[(xr + diag) | 0]++;
}
let sc = 0;
for (let i = 0; i < 2 * diag + 2; i++) sc += hist[i] * hist[i] + hist2[i] * hist2[i];
return sc;
};
const trial = [];
let best = 0, bestS = -1;
for (let d = -maxDeg; d <= maxDeg + 1e-9; d += 0.25) { const s = score(d); trial.push(s); if (s > bestS) { bestS = s; best = d; } }
const fine0 = best;
for (let d = fine0 - 0.25; d <= fine0 + 0.25 + 1e-9; d += 0.02) { const s = score(d); if (s > bestS) { bestS = s; best = d; } }
trial.sort((p, q) => p - q);
const med = trial[trial.length >> 1] || 1;
return { deg: Math.round(best * 100) / 100, confidence: bestS / med, points: n };
}
Rotating the image
64-bit
Turns a grey image by an angle about its middle.
source = c + R(angle) × (pixel − c); value = bilinear mix of the 4 nearest source pixels
The source position and the bilinear weights are 64-bit; the result is rounded to a whole grey level with Math.round because pixels are 8-bit. A black-and-white image is rotated as grey and cut again at level 127.
raster/image.js
lines 201 to 215
CadCraft 0.8.2
export function rotateGrey(grey, deg, { fill = 255 } = {}) {
const { w, h, data } = grey;
checkSize(w, h);
const out = new Uint8Array(w * h), a = deg * Math.PI / 180, c = Math.cos(a), s = Math.sin(a), cx = (w - 1) / 2, cy = (h - 1) / 2;
for (let y = 0; y < h; y++) for (let x = 0; x < w; x++) {
// the source point that lands on (x, y) when the picture turns anticlockwise on screen (y down)
const dx = x - cx, dy = y - cy;
const sx = cx + dx * c - dy * s, sy = cy + dx * s + dy * c;
if (sx < -0.5 || sy < -0.5 || sx > w - 0.5 || sy > h - 0.5) { out[y * w + x] = fill; continue; }
const x0 = Math.floor(sx), y0 = Math.floor(sy), fx = sx - x0, fy = sy - y0;
const g = (xx, yy) => data[Math.min(h - 1, Math.max(0, yy)) * w + Math.min(w - 1, Math.max(0, xx))];
out[y * w + x] = Math.round((g(x0, y0) * (1 - fx) + g(x0 + 1, y0) * fx) * (1 - fy) + (g(x0, y0 + 1) * (1 - fx) + g(x0 + 1, y0 + 1) * fx) * fy);
}
return { w, h, data: out };
}
Turns the stored frequency coefficients of an 8 × 8 block back into pixel values.
pixel = round(Σ_u Σ_v C(x,u)·C(y,v)·coef(u,v) + 128), clamped to 0–255
A direct (not fast) transform in two passes with a 64-entry Float64Array cosine table, so no fixed-point approximation enters. Each pixel is rounded with Math.round and clamped to 0 to 255.
raster/jpeg.js
lines 138 to 142
CadCraft 0.8.2
// inverse DCT, rows then columns
for (let y = 0; y < 8; y++) for (let x = 0; x < 8; x++) { let s2 = 0; for (let u = 0; u < 8; u++) s2 += COS[x * 8 + u] * coef[y * 8 + u]; tmp[y * 8 + x] = s2; }
for (let x = 0; x < 8; x++) for (let y = 0; y < 8; y++) { let s2 = 0; for (let v = 0; v < 8; v++) s2 += COS[y * 8 + v] * tmp[v * 8 + x]; blk[y * 8 + x] = s2; }
const o = by * 8 * c.pw + bx * 8;
for (let y = 0; y < 8; y++) for (let x = 0; x < 8; x++) { const v = Math.round(blk[y * 8 + x] + 128); c.plane[o + y * c.pw + x] = v < 0 ? 0 : v > 255 ? 255 : v; }
JPEG color conversion
64-bit
Converts the decoded luminance and color-difference values to red, green and blue.
R = Y + 1.402·Cr; G = Y − 0.344136·Cb − 0.714136·Cr; B = Y + 1.772·Cb
The three values are clamped to 0 to 255 and assigned to a Uint8Array without rounding, so the fraction is cut off and not rounded to the nearest level. The coefficients are the standard's six-digit values.
raster/jpeg.js
lines 172 to 173
CadCraft 0.8.2
const r = Y + 1.402 * cr, g = Y - 0.344136 * cb - 0.714136 * cr, bl = Y + 1.772 * cb;
data[o] = r < 0 ? 0 : r > 255 ? 255 : r; data[o + 1] = g < 0 ? 0 : g > 255 ? 255 : g; data[o + 2] = bl < 0 ? 0 : bl > 255 ? 255 : bl;
In detail
Method
Otsu's threshold, Sauvola's local threshold from integral images, a projection-profile search for skew, bilinear resampling for rotation, and a baseline JPEG decoder the file attributes to ITU-T T.81.
Number types
8-bit pixels in Uint8Array; integral images and histograms in Float64Array; decoder coefficients in Int32Array and Float64Array.
Tolerances
Skew search: 0.25° steps over ±10°, then 0.02° steps over ±0.25°; result rounded to 0.01°. Skew sampling: at most about 150 000 ink points.
Limits
Rotation resamples pixels, so each turn blurs the image slightly. Chroma is upsampled by nearest neighbor.
Approach: 64-bit, stated limit
The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.