ref: f9bb63633ce47c39a98d929a98b8392c98b07e87
dir: /rast.c/
/*
* Recommended reading:
*
* Wavelet Rasterization
* (2011) J. Manson, S.Schaefer
*
* A Simple and Fast Line-Clipping Method as a Scratch Extension for Computer Graphics Education
* (2019) Dimitrios Matthes, Vasileios Drakopoulos
*
* Clipping of Bézier Curves
* (1985) Fuhua Cheng, Chi-Cheng Lin
*
* On the numerical condition of polynomials in Bernstein form
* (1987) R. T. Farouki, V. T. Rajan
*
* Numerical Condition of Polynomials in Different Forms
* (2001) Hong Zhang
*
* Computer Aided Geometric Design
* (2017) Thomas W. Sederberg
*/
#include "otfsys.h"
#include "otf.h"
typedef struct SegQ SegQ;
typedef struct SegC SegC;
typedef struct Spt Spt;
typedef double Sval;
struct Spt {
Sval x;
Sval y;
};
struct SegQ {
union {
struct {
union {
Spt p0;
Sval v0[2];
};
union {
Spt p1;
Sval v1[2];
};
union {
Spt p2;
Sval v2[2];
};
};
Sval v[3*2];
};
};
/* the divisions when calculating coefficients are reduced,
* thus pixel value range is [0.0-QBZR_PIX_SCALE] instead of [0.0-1.0]
*/
#define QBZR_PIX_SCALE 4.0
#define ε 1.0e-10
static int
closest²(int v)
{
v--;
v |= v>>1;
v |= v>>2;
v |= v>>4;
v |= v>>8;
v |= v>>16;
return v + 1;
}
static int
qslv(Sval *qs, Sval a, Sval b, Sval c)
{
Sval p, q, d;
int n;
#define ₀(x) (x > -ε && x < ε)
if(₀(a)){
if(₀(b))
return 0;
qs[0] = -c/b;
return qs[0] > 0 && qs[0] < 1;
}
p = b/(2.0*a);
q = c/a;
d = p*p - q;
if(₀(d)){
qs[0] = -p;
return qs[0] > 0 && qs[0] < 1;
}
#undef ₀
if(d < 0.0)
return 0;
d = sqrt(d);
n = 0;
if((q = -d-p) > 0 && q < 1)
qs[n++] = q;
if((q = d-p) > 0 && q < 1)
qs[n++] = q;
return n;
}
static int
Svalcmp(void *a_, void *b_)
{
Sval a = *(Sval*)a_, b = *(Sval*)b_;
return (a > b) - (a < b);
}
static int
qKL(SegQ *s₀, int jj, Sval px, Sval py, Sval *K, Sval *L)
{
Sval qs[1+4*2+1];
Sval a, b, c, w, v;
SegQ s, z;
int i, j, n, r;
/* transform */
for(i = 0; i < nelem(s.v); i += 2){
s.v[i+0] = s₀->v[i+0]*jj - px;
s.v[i+1] = s₀->v[i+1]*jj - py;
}
/* FIXME would it make things faster to do proper convex hull test here? */
if(s.p0.x < 0 && s.p1.x < 0 && s.p2.x < 0 ||
s.p0.x > 1 && s.p1.x > 1 && s.p2.x > 1 ||
s.p0.y < 0 && s.p1.y < 0 && s.p2.y < 0 ||
s.p0.y > 1 && s.p1.y > 1 && s.p2.y > 1)
return 0;
#define e(t,a) (s.p0.a*(1-t)*(1-t) + 2*s.p1.a*(1-t)*t + s.p2.a*t*t)
#define within(v) ((w = e(v, x)) >= -ε && w <= 1+ε && (w = e(v, y)) >= -ε && w <= 1+ε)
/* clip to quadtree cell */
n = 0;
if(s.p0.x >= 0 && s.p0.x <= 1 && s.p0.y >= 0 && s.p0.y <= 1)
qs[n++] = 0;
for(i = 0; i < 2; i++){
c = s.v0[i];
a = c - 2*s.v1[i] + s.v2[i];
b = 2*(s.v1[i] - c);
n += qslv(qs+n, a, b, c);
n += qslv(qs+n, a, b, c-1);
}
qsort(qs, n, sizeof(Sval), Svalcmp);
if(s.p2.x >= 0 && s.p2.x <= 1 && s.p2.y >= 0 && s.p2.y <= 1)
qs[n++] = 1;
j = 0;
for(i = 0; i < n; i++){
v = qs[i];
if(within(v))
qs[j++] = v;
}
#undef within
/* remove duplicates */
n = 1;
for(i = 1; i < j; i++){
v = qs[i];
if(v != qs[i-1])
qs[n++] = v;
}
/* calculate K & L on subsections */
r = 0;
for(i = 0; i < n-1; i++){
a = qs[i];
b = qs[i+1];
z.p0.y = e(a, y);
z.p2.y = e(b, y);
if(fabs(z.p0.y-1) < ε && fabs(z.p2.y-1) < ε)
continue;
z.p0.x = e(a, x);
z.p2.x = e(b, x);
if(fabs(z.p0.x-1) < ε && fabs(z.p2.x-1) < ε)
continue;
#undef e
c = 1 - a;
v = (b - a)/c;
z.p1.x = (1-v)*z.p0.x + v*(c*s.p1.x + a*s.p2.x);
z.p1.y = (1-v)*z.p0.y + v*(c*s.p1.y + a*s.p2.y);
K[0] += z.p2.y - z.p0.y;
K[1] += z.p0.x - z.p2.x;
L[0] += 3*z.p2.x*z.p2.y
+ 2*z.p2.y*z.p1.x
- 2*z.p2.x*z.p1.y
+ z.p2.y*z.p0.x
+ 2*z.p1.y*z.p0.x
- z.p0.y*( z.p2.x + 2*z.p1.x + 3*z.p0.x);
L[1] += 2*z.p1.y*z.p0.x
+ z.p2.y*(2*z.p1.x + z.p0.x)
- 2*z.p1.x*z.p0.y
+ 3*z.p0.x*z.p0.y
- z.p2.x*(3*z.p2.y + 2*z.p1.y + z.p0.y);
r = 1;
}
if(r){
L[0] /= 6.0;
L[1] /= 6.0;
}
return r;
}
static int
lKL(SegQ *s₀, int jj, Sval px, Sval py, Sval *K, Sval *L)
{
Sval x₁, x₂, y₁, y₂, α₁, α₂, β₁, β₂;
/* transform */
α₁ = x₁ = s₀->p0.x*jj - px;
β₁ = y₁ = s₀->p0.y*jj - py;
α₂ = x₂ = s₀->p1.x*jj - px;
β₂ = y₂ = s₀->p1.y*jj - py;
if(α₁ < 0){
if(α₂ < 0) return 0;
α₁ = 0;
β₁ = (y₂-y₁)/(x₂-x₁)*(0-x₁)+y₁;
}else if(α₁ > 1){
if(α₂ > 1) return 0;
α₁ = 1;
β₁ = (y₂-y₁)/(x₂-x₁)*(1-x₁)+y₁;
}
if(α₂ < 0){
α₂ = 0;
β₂ = (y₂-y₁)/(x₂-x₁)*(0-x₁)+y₁;
}else if(α₂ > 1){
α₂ = 1;
β₂ = (y₂-y₁)/(x₂-x₁)*(1-x₁)+y₁;
}
if(β₁ < 0){
if(β₂ < 0) return 0;
β₁ = 0;
α₁ = (x₂-x₁)/(y₂-y₁)*(0-y₁)+x₁;
}else if(β₁ > 1){
if(β₂ > 1) return 0;
β₁ = 1;
α₁ = (x₂-x₁)/(y₂-y₁)*(1-y₁)+x₁;
}
if(β₂ < 0){
β₂ = 0;
α₂ = (x₂-x₁)/(y₂-y₁)*(0-y₁)+x₁;
}else if(β₂ > 1){
β₂ = 1;
α₂ = (x₂-x₁)/(y₂-y₁)*(1-y₁)+x₁;
}
K[0] += β₂ - β₁;
K[1] += α₁ - α₂;
L[0] += (β₂ - β₁)*(α₁ + α₂)/2;
L[1] += (α₁ - α₂)*(β₁ + β₂)/2;
return 1;
}
static u64int
qCxy(SegQ *s, int ns, int jj, int px, int py, Sval *c, u64int *m, u64int tm)
{
int (*f)(SegQ*, int, Sval, Sval, Sval*, Sval*);
u64int tx₀, tx₁, z, all;
Sval K[4][2], L[4][2];
int i;
jj *= 2;
px *= 2;
py *= 2;
if(jj > 8){
tx₀ = 0;
tx₁ = 0;
}else{
tx₀ = 1ULL<<(px*jj + py);
tx₁ = 1ULL<<((px+1)*jj + py);
}
all = 0;
for(i = 0; i < ns; i++, s++){
if((m[i] & tm) == 0)
continue;
K[0][0] = K[0][1] = 0;
K[1][0] = K[1][1] = 0;
K[2][0] = K[2][1] = 0;
K[3][0] = K[3][1] = 0;
L[0][0] = L[0][1] = 0;
L[1][0] = L[1][1] = 0;
L[2][0] = L[2][1] = 0;
L[3][0] = L[3][1] = 0;
if(s->p1.x == s->p2.x && s->p1.y == s->p2.y)
f = lKL;
else
f = qKL;
z = 0;
if(tx₀ == 0){
z |= f(s, jj, px+0, py+0, K[0], L[0]);
z |= f(s, jj, px+0, py+1, K[1], L[1]);
z |= f(s, jj, px+1, py+0, K[2], L[2]);
z |= f(s, jj, px+1, py+1, K[3], L[3]);
}else{
z = 0;
if(f(s, jj, px+0, py+0, K[0], L[0]))
z |= tx₀;
if(f(s, jj, px+0, py+1, K[1], L[1]))
z |= tx₀<<1;
if(f(s, jj, px+1, py+0, K[2], L[2]))
z |= tx₁;
if(f(s, jj, px+1, py+1, K[3], L[3]))
z |= tx₁<<1;
m[ns+i] |= z;
all |= z;
}
if(z != 0){
c[0] += L[0][1] + L[2][1] + K[1][1] - L[1][1] + K[3][1] - L[3][1];
c[1] += L[0][0] + L[1][0] + K[2][0] - L[2][0] + K[3][0] - L[3][0];
c[2] += L[0][0] - L[1][0] + K[2][0] - L[2][0] - K[3][0] + L[3][0];
}
}
return all;
}
static int
qbzr(SegQ *seg, int ns, int dw, int dh, Sval *p)
{
int ₂ju, ₂j, w, h, i, j, ju, px, py;
Sval r₂ju, s₀, s, *c, *c₀x, *c₀y, *cs, zx, zy;
u64int *ma, *mb, all, nall;
₂ju = closest²(dh);
if(₂ju < (j = closest²(dw)))
₂ju = j;
r₂ju = 1.0 / ₂ju;
for(ju = 1; (1<<ju) < ₂ju; ju++);
/* scale */
for(i = 0; i < ns; i++){
for(j = 0; j < nelem(seg->v); j++)
seg[i].v[j] *= r₂ju;
}
/* calculate area */
s₀ = 0;
for(i = 0; i < ns; i++){
#define det(a,b) (a.x*b.y - a.y*b.x)
#define cQ(s) (det(s.p0, s.p1) + det(s.p1, s.p2) + det(s.p0, s.p2)/2.0)/3.0
s₀ -= cQ(seg[i]);
#undef det
#undef cQ
}
s₀ *= QBZR_PIX_SCALE;
/* working space:
* quick is-segment-in-the-cell tests (two 64-bit nums per segment)
* wavelet coefficients (count by geometric series, 3 per "pixel" of every zoom level)
*/
j = 2*ns*sizeof(*ma) + (1-(4<<2*(ju-1)))/(1-4)*3*sizeof(*c);
/* current zoom level bitmasks for quick tests */
if((ma = calloc(1, j)) == nil){
werrstr("no memory");
return -1;
}
/* next zoom level bitmasks */
mb = ma + ns;
/* first (1x1) zoom level - all segments are there */
for(i = 0; i < ns; i++)
ma[i] = 1;
/* coefficients */
c = cs = (Sval*)(ma + 2*ns);
/* bitmasks used per-segment are united (bitwise OR) to skip entire empty cells */
all = 1;
nall = 0;
/* first few loops build the test bitmasks: 1x1 -> 4x4 -> 8x8 */
for(j = 0, ₂j = 1; j < ju && j <= 3; j++, ₂j <<= 1){
for(px = 0; px < ₂j; px++){
for(py = 0; py < ₂j; py++, c += 3){
u64int tm = 1ULL<<(px*₂j + py);
if(all & tm)
nall |= qCxy(seg, ns, ₂j, px, py, c, ma, tm);
}
}
if(j != 3){
for(i = 0; i < ns; i++){
ma[i] = mb[i];
mb[i] = 0;
}
all = nall;
nall = 0;
}
}
/* no more bitmasks, just calculate coefficients */
for(; j < ju; j++, ₂j <<= 1){
for(px = 0; px < ₂j; px++){
/* bitmasks are expanded only up to 8x8 (64 bits)
* so use that for more coarse testing
*/
u64int tm₀ = 1ULL<<((px>>(j-3))*8);
for(py = 0; py < ₂j; py++, c += 3){
u64int tm = tm₀ << (py>>(j-3));
if(all & tm)
qCxy(seg, ns, ₂j, px, py, c, ma, tm);
}
}
}
/* reverse Y axis (opentype coord vs plan 9) */
for(h = dh-1; h >= 0; h--){
for(w = 0; w < dw; w++){
c = cs;
s = s₀;
for(₂j = 1; ₂j < ₂ju; ₂j *= 2){
zx = r₂ju * ₂j * w;
c₀x = c;
px = (int)zx;
zx -= px;
c += px*3*₂j;
for(; zx >= 0; px++, zx--, c = c₀y + 3*₂j){
c₀y = c;
zy = r₂ju * ₂j * h;
py = (int)zy;
zy -= py;
c += 3*py;
for(; zy >= 0; py++, zy--, c += 3){
switch((zy < 0.5)<<1 | (zx < 0.5)){
case 0: s += +c[0] + c[1] - c[2]; break;
case 1: s += +c[0] - c[1] + c[2]; break;
case 2: s += -c[0] + c[1] + c[2]; break;
default: s += -c[0] - c[1] - c[2]; break;
}
}
}
c = c₀x + 3*₂j*₂j;
}
*p++ = s;
}
}
free(ma);
return 0;
}
u8int *
otfdrawglyf(Glyf *g, int h, int *wo)
{
SegQ *s₀, *s;
Sval *fp;
Point *p₀, *p, *prev;
Spt *pts₀, *pts, *e;
int w, i, r, k, npts, npx, ns;
double scale;
/* FIXME those epsilons and ceil is a rather dumb hack for now */
scale = (double)h / (g->yMax - g->yMin) + ε;
w = *wo = ceil((g->xMax - g->xMin) * scale + ε);
npx = w * h;
fp = calloc(1, npx*sizeof(*fp));
ns = 0;
s₀ = nil;
for(k = 0; k < g->numberOfContours; k++){
npts = g->simple->endPtsOfContours[k]+1;
if(k > 0)
npts -= g->simple->endPtsOfContours[k-1]+1;
pts₀ = e = malloc((npts+2)*3*sizeof(*e));
p₀ = p = g->simple->points + (k > 0 ? g->simple->endPtsOfContours[k-1]+1 : 0);
prev = p₀+npts-1;
if(!p->onCurve){
/* fun stuff */
if(prev->onCurve)
*e = (Spt){prev->x, prev->y};
else
*e = (Spt){(Sval)(p->x + prev->x)/2, (Sval)(p->y + prev->y)/2};
e++;
}
for(prev = nil; p < p₀+npts; prev = p, p++, e++){
if(prev != nil && p->onCurve == prev->onCurve){
/* more fun stuff */
if(p->onCurve) /* straight line */
*e = (Spt){p->x, p->y};
else /* a point in the middle */
*e = (Spt){(Sval)(p->x + prev->x)/2, (Sval)(p->y + prev->y)/2};
e++;
}
*e = (Spt){p->x, p->y};
}
if(e[-1].x != pts₀->x || e[-1].y != pts₀->y){
if(p[-1].onCurve) /* close with a straight line */
*e++ = *pts₀;
*e++ = *pts₀;
}
s₀ = realloc(s₀, (ns + (e-pts₀)/2)*sizeof(*s₀));
s = s₀ + ns;
pts = e;
for(e = pts₀; e <= pts-3; e += 2, s++, ns++){
s->v0[0] = (e[0].x - g->xMin) * scale;
s->v0[1] = (e[0].y - g->yMin) * scale;
s->v1[0] = (e[1].x - g->xMin) * scale;
s->v1[1] = (e[1].y - g->yMin) * scale;
s->v2[0] = (e[2].x - g->xMin) * scale;
s->v2[1] = (e[2].y - g->yMin) * scale;
}
free(pts₀);
}
r = qbzr(s₀, ns, w, h, fp);
free(s₀);
if(r != 0)
return nil;
u8int *b = malloc(npx);
for(i = 0; i < npx; i++){
b[i] = fp[i] <= 0 ? 0 : (fp[i] >= QBZR_PIX_SCALE ? 255 : fp[i]/QBZR_PIX_SCALE*255);
b[i] = 255 - b[i];
}
free(fp);
return b;
}