进展！
自己写的网页终于功能全部完成啦！
现在已经把成品网页放到博客上来了，上面几个链接最右边那个 PageDewarp 就是，戳进去就行了。
放几张成品：
效果还是相当不错的。如果有二值化就更好了。
卷曲复原
在确定了曲面的各个参数之后，整个纸面的曲面就有了，剩下的就是把卷曲复原，输出结果。
卷曲复原的过程是一个图像上坐标映射的过程。我们遍历结果图像上的每一个坐标，唯一地确立其在输入图像上的源坐标，把源坐标的像素复制进来（或许再使用一些采样的技巧），拼成最终的结果，就成了。
更具体一点，如果输出图像是 \(w\times h\)，我们就需要把页面曲面水平分成 \(w\) 等分，垂直分成 \(h\) 等分，形成一个网格。结果图像上的 \((i, j)\) 像素就是网格中 \((i,j)\) 对应点的三维坐标在原始图像上投影点对应的像素。
垂直等分是简单的，因为在最早的假设里面我就假设页面的垂直边沿是直线。
但是水平边沿是一条曲线，或者更精确的说，可以被近似看为一条分段三次样条曲线。怎么等分这样一个曲线呢？
这种事情用数学解决是最理想的，只要能够解析地计算出任意区间曲线的长度，配合二分这个问题就解决了。但是三次函数曲线长度的积分是非初等的且无比复杂的，所以数学方法行不通。
（Well，好像 Mathematica 告诉我算还是可以算的……）
于是只能使用不那么优雅的数值方法解决了。飞快地写了一个 Simpson Integration，效果看起来还不错，嘛反正精度不够缩小步长就可以了，和后面的计算比这个又不是瓶颈。
// spline: 两端为 (0, 0) 和 (1, 0) 的分段样条曲线
// res: 等分的数量
const step = 0.0001;
let len = 0, t = [];
for (let i = 0; i <= this.nodeCount; i++) {
let [a, b, c, d] = spline[i]; // a + bx + cx^2 + dx^3
let xLeft = i === 0 ? 0 : this.splineNodes[i - 1][0];
let xRight = i === this.nodeCount ? 1 : this.splineNodes[i][0];
let f = x => Math.hypot(1, ((3 * d * x + 2 * c) * x + b));
for (let x = xLeft; x < xRight; x += step) {
let h = Math.min(xRight - x, step);
len += (f(x) + 4 * f(x + h / 2) + f(x + h)) * h / 6;
t.push([x, len, ((d * x + c) * x + b) * x + a]);
}
}
let ret = [], p = 0;
for (let i = 1; i <= res; i++) {
let l = i / res * len;
while (p + 1 < t.length && t[p + 1][1] <= l) p++;
let lt = t[p], rt = p + 1 === t.length ? [1, len, 0] : t[p + 1];
ret.push(l - lt[1] <= rt[1] - l ? [lt[0], lt[2]] : [rt[0], rt[2]]);
}
// ret 就是曲线上的点了anyway，反正积的函数也没有那么畸形，也没那么着急优化，代码写得乱一点也就凑活着过去了（疯狂为自己的懒找借口）
求出曲线一连串等分点之后就可以做复原了：
// vadd, multiply: 向量和，数乘
// criticalPoints: 曲线等分点
for (let x = 0; x < resultWidth; x++) { // 遍历结果的每一列像素
const top = vadd(this.upperLeft3D,
multiply(this.baseVec, criticalPoints[x][0]), multiply(this.normalVec, criticalPoints[x][1]));
for (let y = 0; y < resultHeight; y++) {
let coord = multiply(this.toCanvasCoord(this.project(
vadd(top, multiply(vertical, y / resultHeight)))), width / this.canvasWidth); // (x, y) 对应的图像的原坐标
let di = (y * resultWidth + x) * 4;
let si = (Math.round(coord[1]) * width + Math.round(coord[0])) * 4;
// R, G, B, alpha 四通道
imageDataResult.data[di] = imageDataOriginal.data[si];
imageDataResult.data[di + 1] = imageDataOriginal.data[si + 1];
imageDataResult.data[di + 2] = imageDataOriginal.data[si + 2];
imageDataResult.data[di + 3] = imageDataOriginal.data[si + 3];
}
}大概是 OI 写惯了觉得 \(10^3\times 10^3\) 级别的运算没啥，直觉告诉我这一段代码跑得会很快。然而现实是残酷的，这段代码要跑三十秒左右，也不知道是哪里出了问题，V8 的 JIT 都拯救不了。
目前只能先用这么慢的代码将就着了，回头有空用 WebAssembly 试图优化一下吧。接下来应该都是工程问题了。