ref: 2cfe3c2c1cb0de08ef9a4dca2d22dac286736b09
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"
#include "otfpriv.h"
#pragma clang float_control(precise, on)
typedef struct SegQ SegQ;
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-9
#define is₀(v) ((v) > -ε && (v) < ε)
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;
if(unlikely(is₀(a))){
if(is₀(b))
return 0;
qs[0] = -c/b;
return qs[0] > 0 && qs[0] < 1;
}
p = b/(2.0*a);
d = p*p - c/a;
if(d < ε)
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 */
s.p0.x = s₀->p0.x*jj - px;
s.p0.y = s₀->p0.y*jj - py;
s.p1.x = s₀->p1.x*jj - px;
s.p1.y = s₀->p1.y*jj - py;
s.p2.x = s₀->p2.x*jj - px;
s.p2.y = s₀->p2.y*jj - py;
#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;
c = s.p0.x;
a = c - 2*s.p1.x + s.p2.x;
b = 2*(s.p1.x - c);
n += qslv(qs+n, a, b, c);
n += qslv(qs+n, a, b, c-1);
c = s.p0.y;
a = c - 2*s.p1.y + s.p2.y;
b = 2*(s.p1.y - 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(is₀(z.p0.y) && is₀(z.p2.y))
continue;
z.p0.x = e(a, x);
z.p2.x = e(b, x);
if(is₀(z.p0.x) && is₀(z.p2.x))
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
Cxy(SegQ *s, int ns, int jj, int px, int py, Sval *c, u64int *m, u64int tm)
{
int (*f)(SegQ*, int, Sval, Sval, Sval*, Sval*);
Sval K[4][2], L[4][2], q[6], j;
u64int tx₀, tx₁, z, all;
u8int w;
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);
}
j = 1.0/(Sval)jj;
q[0] = j*px;
q[1] = j*py;
q[2] = q[0]+j;
q[3] = q[1]+j;
q[4] = q[2]+j;
q[5] = q[3]+j;
all = 0;
for(i = 0; i < ns; i++, s++){
if((m[i] & tm) == 0 ||
s->p0.x <= q[0] && s->p1.x <= q[0] && s->p2.x <= q[0] ||
s->p0.x >= q[4] && s->p1.x >= q[4] && s->p2.x >= q[4] ||
s->p0.y <= q[1] && s->p1.y <= q[1] && s->p2.y <= q[1] ||
s->p0.y >= q[5] && s->p1.y >= q[5] && s->p2.y >= q[5])
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;
z = 0;
f = s->p1.x == s->p2.x && s->p1.y == s->p2.y ? lKL : qKL;
w =
(s->p0.x <= q[2] || s->p1.x <= q[2] || s->p2.x <= q[2])<<0 |
(s->p0.x >= q[2] || s->p1.x >= q[2] || s->p2.x >= q[2])<<1 |
(s->p0.y <= q[3] || s->p1.y <= q[3] || s->p2.y <= q[3])<<2 |
(s->p0.y >= q[3] || s->p1.y >= q[3] || s->p2.y >= q[3])<<3;
if(tx₀ == 0){
if((w & 5) == 5) z |= f(s, jj, px+0, py+0, K[0], L[0]);
if((w & 6) == 6) z |= f(s, jj, px+1, py+0, K[2], L[2]);
if((w & 9) == 9) z |= f(s, jj, px+0, py+1, K[1], L[1]);
if((w & 10) == 10) z |= f(s, jj, px+1, py+1, K[3], L[3]);
}else{
if((w & 5) == 5 && f(s, jj, px+0, py+0, K[0], L[0])) z |= tx₀;
if((w & 6) == 6 && f(s, jj, px+1, py+0, K[2], L[2])) z |= tx₁;
if((w & 9) == 9 && f(s, jj, px+0, py+1, K[1], L[1])) z |= tx₀<<1;
if((w & 10) == 10 && 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[1][1] + L[2][1] - L[3][1] + K[1][1] + K[3][1];
c[1] += L[0][0] + L[1][0] - L[2][0] - L[3][0] + K[2][0] + K[3][0];
c[2] += L[0][0] - L[1][0] - L[2][0] + L[3][0] + K[2][0] - K[3][0];
}
}
return all;
}
static int
qbzr(SegQ *seg, int ns, int dw, int dh, u8int *p)
{
int ju², j², w, h, i, j, ju, px, py;
Sval rju², 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;
rju² = 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] *= rju²;
}
/* calculate area */
s₀ = 0;
for(i = 0; i < ns; i++){
#define det(a,b) (a.x*b.y - a.y*b.x)
#define cL(s) det(s.p0, s.p2)/2.0
#define cQ(s) (2*det(s.p0, s.p1) + 2*det(s.p1, s.p2) + det(s.p0, s.p2))/6.0
if(seg[i].p1.x == seg[i].p2.x && seg[i].p1.y == seg[i].p2.y)
s₀ -= cL(seg[i]);
else
s₀ -= cQ(seg[i]);
#undef det
#undef cL
#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 |= Cxy(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)
Cxy(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 = rju² * 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 = rju² * 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++ = 255 -
(s < QBZR_PIX_SCALE/255 ? 0 : (s >= QBZR_PIX_SCALE ? 255 : s/QBZR_PIX_SCALE*255));
}
}
free(ma);
return 0;
}
static SegQ *
addpoints(SegQ *s, Spt *pts₀, Glyf *g, Sval dx, Sval dy, Sval sx, Sval sy, Sval *bb)
{
int k, npts, xMin, yMin, xMax, yMax;
Point *p₀, *p, *prev;
Spt *e, *pts;
Sval v;
xMin = yMin = 99999;
xMax = yMax = -99999;
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(prev->onCurve)
*e = (Spt){dx+sx*prev->x, dy+sy*prev->y};
else
*e = (Spt){dx+sx/2*(p->x+prev->x), dy+sy/2*(p->y+prev->y)};
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){dx+sx*p->x, dy+sy*p->y};
else /* a point in the middle */
*e = (Spt){dx+sx/2*(p->x+prev->x), dy+sy/2*(p->y+prev->y)};
e++;
}
*e = (Spt){dx+sx*p->x, dy+sy*p->y};
if(bb != nil){
if(xMin > p->x) xMin = p->x;
if(yMin > p->y) yMin = p->y;
if(xMax < p->x) xMax = p->x;
if(yMax < p->y) yMax = p->y;
}
}
/* close with a straight line */
if(p[-1].onCurve)
*e++ = *pts;
*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;
}
}
if(bb != nil){
if(bb[0] > (v = dx+sx*xMin)) bb[0] = v;
if(bb[1] > (v = dy+sy*yMin)) bb[1] = v;
if(bb[2] < (v = dx+sx*xMax)) bb[2] = v;
if(bb[3] < (v = dy+sy*yMax)) bb[3] = v;
}
return s;
}
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;
}
GlyfImage *
otfdrawglyf(Otf *o, Glyf *g, double ppem)
{
int i, j, maxptstotal, maxpts, ngs, w, h, npx, baseline;
Glyf *gs[MAXCOMPONENTS];
ComponentGlyph *cg;
Sval scale, offY;
SegQ *s₀, *s, *p;
GlyfImage *im;
Sval bb[4];
Spt *pts;
ngs = 0;
pts = nil;
im = 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;
}
}
*/
pts = calloc(1, maxpts*sizeof(*pts) + maxptstotal/2*sizeof(*s));
scale = ppem / o->td.head->unitsPerEm;
bb[0] = g->xMin*scale;
bb[1] = g->yMin*scale;
bb[2] = g->xMax*scale;
bb[3] = g->yMax*scale;
s = s₀ = (SegQ*)(pts + maxpts);
cg = g->component;
for(i = 0; i < ngs; i++){
Sval 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, gs[i], dx, dy, gscaleX, gscaleY, ngs > 1 ? bb : nil)) == nil)
goto done;
}
if(s == s₀){
w = h = 1;
baseline = 0;
}else{
/* height+baseline is where it is in the image */
baseline = floor(bb[1]);
offY = -baseline;
if(offY+bb[1] < ε){
/*
* very rarely a glyph will have one of its segments
* go over the edge because of imprecise arithmetic.
* it's much easier to solve this by always having a small gap.
*/
offY += 4*ε;
bb[3] += 8*ε;
}
bb[3] += offY;
w = ceil(bb[2] - bb[0]);
h = ceil(bb[3]);
for(p = s₀; p != s; p++){
p->p0.x -= bb[0];
p->p1.x -= bb[0];
p->p2.x -= bb[0];
p->p0.y += offY;
p->p1.y += offY;
p->p2.y += offY;
}
}
npx = w*h;
im = malloc(sizeof(*im) + npx);
im->b[0] = 0;
if(qbzr(s₀, s-s₀, w, h, im->b) == 0){
im->w = w;
im->h = h;
im->baseline = baseline;
}else{
free(im);
im = nil;
}
done:
for(i = 0; i < ngs; i++){
if(gs[i] != g)
free(gs[i]);
}
free(pts);
return im;
}