ref: 0d17a04bbcbb6ee26415d356f40e75d891c5f016
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 MAXCOMPONENTS 16
#define ε 1.0e-6
#define min(a,b) ((a)<(b)?(a):(b))
#define max(a,b) ((a)>(b)?(a):(b))
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;
}
static SegQ *
addpoints(SegQ *s, Spt *pts₀, int max, Glyf *g)
{
Point *p₀, *p, *prev;
Spt *e, *pts;
int k, npts;
for(k = 0; k < g->numberOfContours; k++){
npts = g->simple->endPtsOfContours[k]+1;
if(k > 0)
npts -= g->simple->endPtsOfContours[k-1]+1;
if(npts < 2)
continue;
e = pts = pts₀;
p₀ = p = g->simple->points + (k > 0 ? g->simple->endPtsOfContours[k-1]+1 : 0);
prev = p₀+npts-1;
if(!p->onCurve){
/* fun stuff */
if(max-- < 1)
goto noplace;
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(max-- < 1)
goto noplace;
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++;
}
if(max-- < 1)
goto noplace;
*e = (Spt){p->x, p->y};
}
if(e > pts && (e[-1].x != pts->x || e[-1].y != pts->y)){
if(p[-1].onCurve){ /* close with a straight line */
if(max-- < 1)
goto noplace;
*e++ = *pts;
}
if(max-- < 1)
goto noplace;
*e++ = *pts;
}
for(; pts <= e-3; pts += 2, s++){
s->v0[0] = pts[0].x;
s->v0[1] = pts[0].y;
s->v1[0] = pts[1].x;
s->v1[1] = pts[1].y;
s->v2[0] = pts[2].x;
s->v2[1] = pts[2].y;
}
}
return s;
noplace:
werrstr("points don't fit");
return nil;
}
static int
glyfmaxnpts(Glyf *g)
{
int k, npts, d;
npts = 0;
for(k = 0; k < g->numberOfContours; k++){
d = g->simple->endPtsOfContours[k]+1;
if(k > 0)
d -= g->simple->endPtsOfContours[k-1]+1;
if(d > 1)
npts += 1+2*d+2; /* two extra to close off the contour */
}
return npts < 2 ? 0 : npts;
}
int
otfdrawglyf(Otf *o, Glyf *g, int ppem, int gap, GlyfImage *im)
{
int i, j, maxptstotal, maxpts, ngs, w, h, r, npx, baseline;
Glyf *gs[MAXCOMPONENTS];
ComponentGlyph *cg;
Sval scale, offY;
SegQ *s₀, *s, *p;
Spt *pts;
Sval *fp;
u8int *b;
r = -1;
ngs = 0;
pts = nil;
fp = nil;
if(g->simple != nil){
gs[ngs++] = g;
}else{
for(cg = g->component; cg != nil && ngs < nelem(gs); cg = cg->next, ngs++){
if((gs[ngs] = otfglyf(o, cg->glyphIndex)) == nil)
goto done;
}
}
for(maxptstotal = maxpts = i = 0; i < ngs; i++){
j = glyfmaxnpts(gs[i]);
maxptstotal += j;
if(maxpts < j)
maxpts = j;
}
/*
FIXME metrics (lsb+rsb)
for(cg = g->component, i = 0; cg != nil && i < ngs; cg = cg->next, i++){
if(cg->flags & CGLYPH_FL_METRICS){
...
break;
}
}
*/
double xMin = 99999, yMin = 99999, xMax = -99999, yMax = -99999;
scale = (Sval)ppem / otfupem(o);
pts = calloc(1, maxpts*sizeof(*pts) + maxptstotal/2*sizeof(*s));
p = s = s₀ = (SegQ*)(pts + maxpts);
cg = g->component;
for(i = 0; i < ngs; i++){
double dx = 0, dy = 0, gscaleX = scale, gscaleY = scale;
if(cg != nil){
if(cg->flags & CGLYPH_FL_SIGNED_XY){
dx = cg->dx;
dy = cg->dy;
}else{
dx = 0;
dy = 0;
}
if(cg->flags & CGLYPH_FL_SCALE_XY){
gscaleX *= cg->scaleX;
gscaleY *= cg->scaleY;
}
if(cg->flags & CGLYPH_FL_SCALED_COMPONENT_OFFSET){
dx *= gscaleX;
dy *= gscaleY;
}else{
dx *= scale;
dy *= scale;
}
/* FIXME rounding
if(cg->flags & CGLYPH_FL_ROUND_TO_GRID_XY){
...
}
*/
cg = cg->next;
}
if((s = addpoints(s, pts, maxpts, gs[i])) == nil)
goto done;
for(; p < s; p++){
p->p0.x = p->p0.x*gscaleX + dx;
xMin = min(xMin, p->p0.x);
xMax = max(xMax, p->p0.x);
p->p0.y = p->p0.y*gscaleY + dy;
yMin = min(yMin, p->p0.y);
yMax = max(yMax, p->p0.y);
p->p1.x = p->p1.x*gscaleX + dx;
xMin = min(xMin, p->p1.x);
xMax = max(xMax, p->p1.x);
p->p1.y = p->p1.y*gscaleY + dy;
yMin = min(yMin, p->p1.y);
yMax = max(yMax, p->p1.y);
p->p2.x = p->p2.x*gscaleX + dx;
xMin = min(xMin, p->p2.x);
xMax = max(xMax, p->p2.x);
p->p2.y = p->p2.y*gscaleY + dy;
yMin = min(yMin, p->p2.y);
yMax = max(yMax, p->p2.y);
}
}
if(s == s₀){
w = h = 1;
baseline = 0;
}else{
xMin -= gap;
yMin -= gap;
xMax += gap;
yMax += gap;
/* height+baseline is where it is in the image */
baseline = floor(yMin);
offY = -baseline;
yMax += offY;
w = ceil(xMax - xMin);
h = ceil(yMax);
for(p = s₀; p != s; p++){
p->p0.x -= xMin;
p->p1.x -= xMin;
p->p2.x -= xMin;
p->p0.y += offY;
p->p1.y += offY;
p->p2.y += offY;
}
}
npx = w*h;
fp = malloc(npx*sizeof(*fp));
if(qbzr(s₀, s-s₀, w, h, fp) == 0){
b = (u8int*)fp;
for(i = 0; i < npx; i++)
b[i] = 255 - (fp[i] <= ε ? 0 : (fp[i] >= QBZR_PIX_SCALE ? 255 : fp[i]/QBZR_PIX_SCALE*255));
if((b = realloc(fp, npx)) != nil){
im->b = b;
im->w = w;
im->h = h;
im->baseline = baseline;
fp = nil;
r = 0;
}
}
done:
for(i = 0; i < ngs; i++){
if(gs[i] != g)
free(gs[i]);
}
free(pts);
free(fp);
return r;
}