ref: 2a0ed585ce506b55d90e1c10142c2150ffa892e2
parent: 0a6bcc7b83d9b5fac12229cac57207b9dc91ba2f
author: rodri <rgl@antares-labs.eu>
date: Wed Aug 13 19:03:23 EDT 2025
libmemdraw: use fixed-point in memaffinewarp() The previous version relied on libgeometry's invm(2) and xform(2) to do its warping, which use double-precision floating-point numbers. They have been imported and its operations replaced by fixed-point ones, so that it can be included in the kernel for future use in draw(3) (most of the architectures lack support for floating-point at the moment). A new data structure, Warp, is now available to store the affine matrix; along with a helper routine, mkwarp(), to translate libgeometry Matrices into Warps that have the correct format necessary to do the transformations.
--- a/sys/include/memdraw.h
+++ b/sys/include/memdraw.h
@@ -7,6 +7,7 @@
typedef struct Memlayer Memlayer;
typedef struct Memcmap Memcmap;
typedef struct Memdrawparam Memdrawparam;
+typedef long Warp[3][3];
#pragma incomplete Memlayer
@@ -150,7 +151,8 @@
extern void memarc(Memimage*, Point, int, int, int, Memimage*, Point, int, int, int);
extern Rectangle memlinebbox(Point, Point, int, int, int);
extern int memlineendsize(int);
-extern int memaffinewarp(Memimage*, Rectangle, Memimage*, Point, uchar*);
+extern int memaffinewarp(Memimage*, Rectangle, Memimage*, Point, Warp);
+extern void mkwarp(Warp, double[3][3]);
extern void _memmkcmap(void);
extern int memimageinit(void);
--- a/sys/man/2/memdraw
+++ b/sys/man/2/memdraw
@@ -25,6 +25,7 @@
memimageline,
memimagedraw,
memaffinewarp,
+mkwarp,
drawclip,
drawclipnorepl,
memlinebbox,
@@ -89,6 +90,8 @@
\fI...\fP
} Memdrawparam;
+typedef long Warp[3][3];
+
.ta \w'\fLMemsubfont* 'u
int drawdebug;
.ft
@@ -135,7 +138,8 @@
void memimagedraw(Memimage *dst, Rectangle r, Memimage *src,
Point sp, Memimage *mask, Point mp, Drawop op)
int memaffinewarp(Memimage *dst, Rectangle r, Memimage *src,
- Point sp, uchar *m)
+ Point sp, Warp w)
+void mkwarp(Warp w, double m[3][3])
.PP
.ft L
.nf
@@ -383,17 +387,24 @@
.PP
.I Memaffinewarp
applies an affine transformation defined by
-.B m
+.B w
to the
.BR src ,
-where
-.B m
-is a 2×3 row-major matrix of
+and stores the result in the
+.BR dst .
+.I Mkwarp
+converts a 3×3 row-major matrix of
.BR double s
(usually a
.IR geometry (2)
-Matrix), and stores the result in the
-.BR dst .
+.BR Matrix )
+into the fixed-point representation required by
+.I memaffinewarp
+in
+.BR w ,
+ignoring the last row which is always set to
+.B "[0 0 1]"
+for affine transformations.
Both
.BR Memimage s
must have the same channel description and use 8-bit channels.
--- a/sys/src/libmemdraw/warp.c
+++ b/sys/src/libmemdraw/warp.c
@@ -2,11 +2,27 @@
#include <libc.h>
#include <draw.h>
#include <memdraw.h>
-#include <geometry.h>
-#define lerp(a,b,t) ((a) + ((((b) - (a))*(t))>>16))
+/* 25.7 fixed-point number operations */
+#define FMASK ((1<<7) - 1)
+#define flt2fix(n) ((long)((n)*(1<<7) + ((n) < 0? -0.5: 0.5)))
+#define int2fix(n) ((vlong)(n)<<7)
+#define fix2int(n) ((n)>>7)
+#define fixmul(a,b) ((vlong)(a)*(vlong)(b) >> 7)
+#define fixdiv(a,b) (((vlong)(a) << 7)/(vlong)(b))
+#define fixfrac(n) ((n)&FMASK)
+
+#define lerp(a,b,t) ((a) + ((((b) - (a))*(t))>>7))
+
+typedef struct Foint Foint;
typedef struct Sampler Sampler;
+
+struct Foint
+{+ long x, y, w;
+};
+
struct Sampler
{Memimage *i;
@@ -16,6 +32,73 @@
int cmask;
};
+static void
+fix_smulm(Warp m, long s)
+{+ m[0][0] = fixmul(m[0][0], s); m[0][1] = fixmul(m[0][1], s); m[0][2] = fixmul(m[0][2], s);
+ m[1][0] = fixmul(m[1][0], s); m[1][1] = fixmul(m[1][1], s); m[1][2] = fixmul(m[1][2], s);
+ m[2][0] = fixmul(m[2][0], s); m[2][1] = fixmul(m[2][1], s); m[2][2] = fixmul(m[2][2], s);
+}
+
+static long
+fix_detm(Warp m)
+{+ return fixmul(m[0][0], (fixmul(m[1][1], m[2][2]) - fixmul(m[1][2], m[2][1])))+
+ fixmul(m[0][1], (fixmul(m[1][2], m[2][0]) - fixmul(m[1][0], m[2][2])))+
+ fixmul(m[0][2], (fixmul(m[1][0], m[2][1]) - fixmul(m[1][1], m[2][0])));
+}
+
+static void
+fix_adjm(Warp m)
+{+ Warp tmp;
+
+ tmp[0][0] = fixmul(m[1][1], m[2][2]) - fixmul(m[1][2], m[2][1]);
+ tmp[0][1] = -fixmul(m[0][1], m[2][2]) + fixmul(m[0][2], m[2][1]);
+ tmp[0][2] = fixmul(m[0][1], m[1][2]) - fixmul(m[0][2], m[1][1]);
+
+ tmp[1][0] = -fixmul(m[1][0], m[2][2]) + fixmul(m[1][2], m[2][0]);
+ tmp[1][1] = fixmul(m[0][0], m[2][2]) - fixmul(m[0][2], m[2][0]);
+ tmp[1][2] = -fixmul(m[0][0], m[1][2]) + fixmul(m[0][2], m[1][0]);
+
+ tmp[2][0] = fixmul(m[1][0], m[2][1]) - fixmul(m[1][1], m[2][0]);
+ tmp[2][1] = -fixmul(m[0][0], m[2][1]) + fixmul(m[0][1], m[2][0]);
+ tmp[2][2] = fixmul(m[0][0], m[1][1]) - fixmul(m[0][1], m[1][0]);
+
+ memmove(m, tmp, sizeof tmp);
+}
+
+/* Cramer's */
+static void
+fix_invm(Warp m)
+{+ long det;
+
+ det = fix_detm(m);
+ if(det == 0)
+ return; /* singular matrices are not invertible */
+ fix_adjm(m);
+ fix_smulm(m, fixdiv(1<<7, det));
+}
+
+static Foint
+fix_xform(Foint p, Warp m)
+{+ return (Foint){+ fixmul(p.x, m[0][0]) + fixmul(p.y, m[0][1]) + fixmul(p.w, m[0][2]),
+ fixmul(p.x, m[1][0]) + fixmul(p.y, m[1][1]) + fixmul(p.w, m[1][2]),
+ fixmul(p.x, m[2][0]) + fixmul(p.y, m[2][1]) + fixmul(p.w, m[2][2])
+ };
+}
+
+void
+mkwarp(Warp w, double m[3][3])
+{+ w[0][0] = flt2fix(m[0][0]); w[0][1] = flt2fix(m[0][1]); w[0][2] = flt2fix(m[0][2]);
+ w[1][0] = flt2fix(m[1][0]); w[1][1] = flt2fix(m[1][1]); w[1][2] = flt2fix(m[1][2]);
+ w[2][0] = 0; w[2][1] = 0; w[2][2] = 1<<7;
+}
+
static ulong
edgehandler(Sampler*, Point*)
{@@ -27,7 +110,8 @@
{ulong c;
- if(ptinrect(p, s->i->r)){+ if(p.x >= s->i->r.min.x && p.x < s->i->r.max.x
+ && p.y >= s->i->r.min.y && p.y < s->i->r.max.y){c = *(ulong*)(s->a + p.y*s->bpl + p.x*s->bpp);
c &= s->cmask;
}else
@@ -36,16 +120,16 @@
}
int
-memaffinewarp(Memimage *d, Rectangle r, Memimage *s, Point sp0, uchar *ma)
+memaffinewarp(Memimage *d, Rectangle r, Memimage *s, Point sp0, Warp ma)
{Sampler samp;
- Matrix m;
+ Warp m;
Point sp, dp;
- Point2 p2;
- double Δx, Δy;
- int Δx2, Δy2;
+ Foint p2;
+ long Δx, Δy;
ulong c00, c01, c10, c11;
uchar c0, c1, c, *a;
+ int bpp, bploff;
if(s->depth < 8 || (s->depth & 3) != 0){ werrstr("unsupported image format");@@ -64,25 +148,23 @@
samp.a = s->data->bdata + s->zero;
samp.bpp = s->depth >> 3;
samp.bpl = sizeof(ulong)*s->width;
- samp.cmask = (1ULL << 8*s->nchan) - 1;
+ samp.cmask = (1ULL << s->depth) - 1;
+ a = d->data->bdata + d->zero;
+ bpp = d->depth >> 3;
+ bploff = sizeof(ulong)*d->width - Dx(d->r)*bpp;
- memmove(m, ma, 2*3*8);
- m[2][0] = 0; m[2][1] = 0; m[2][2] = 1;
- invm(m);
+ memmove(m, ma, sizeof(Warp));
+ fix_invm(m);
- for(dp.y = r.min.y; dp.y < r.max.y; dp.y++)
- for(dp.x = r.min.x; dp.x < r.max.x; dp.x++){- a = byteaddr(d, dp);
+ for(dp.y = r.min.y; dp.y < r.max.y; dp.y++, a += bploff)
+ for(dp.x = r.min.x; dp.x < r.max.x; dp.x++, a += bpp){+ p2 = fix_xform((Foint){int2fix(dp.x - r.min.x), int2fix(dp.y - r.min.y), 1<<7}, m);- p2 = xform((Point2){dp.x - r.min.x, dp.y - r.min.y, 1}, m);- sp.x = p2.x;
- sp.y = p2.y;
+ Δx = fixfrac(p2.x);
+ Δy = fixfrac(p2.y);
- Δx = p2.x - sp.x;
- Δx2 = Δx * (1<<16);
- Δy = p2.y - sp.y;
- Δy2 = Δy * (1<<16);
-
+ sp.x = fix2int(p2.x);
+ sp.y = fix2int(p2.y);
sp.x += sp0.x;
sp.y += sp0.y;
c00 = sample(&samp, sp);
@@ -95,24 +177,24 @@
switch(s->nchan){case 4:
- c0 = c00 >> 8*3 & 0xFF; c0 = lerp(c0, c01 >> 8*3 & 0xFF, Δx2);
- c1 = c10 >> 8*3 & 0xFF; c1 = lerp(c1, c11 >> 8*3 & 0xFF, Δx2);
- c = lerp(c0, c1, Δy2);
+ c0 = c00 >> 8*3 & 0xFF; c0 = lerp(c0, c01 >> 8*3 & 0xFF, Δx);
+ c1 = c10 >> 8*3 & 0xFF; c1 = lerp(c1, c11 >> 8*3 & 0xFF, Δx);
+ c = lerp(c0, c1, Δy);
*(a + 3) = c;
case 3:
- c0 = c00 >> 8*2 & 0xFF; c0 = lerp(c0, c01 >> 8*2 & 0xFF, Δx2);
- c1 = c10 >> 8*2 & 0xFF; c1 = lerp(c1, c11 >> 8*2 & 0xFF, Δx2);
- c = lerp(c0, c1, Δy2);
+ c0 = c00 >> 8*2 & 0xFF; c0 = lerp(c0, c01 >> 8*2 & 0xFF, Δx);
+ c1 = c10 >> 8*2 & 0xFF; c1 = lerp(c1, c11 >> 8*2 & 0xFF, Δx);
+ c = lerp(c0, c1, Δy);
*(a + 2) = c;
case 2:
- c0 = c00 >> 8*1 & 0xFF; c0 = lerp(c0, c01 >> 8*1 & 0xFF, Δx2);
- c1 = c10 >> 8*1 & 0xFF; c1 = lerp(c1, c11 >> 8*1 & 0xFF, Δx2);
- c = lerp(c0, c1, Δy2);
+ c0 = c00 >> 8*1 & 0xFF; c0 = lerp(c0, c01 >> 8*1 & 0xFF, Δx);
+ c1 = c10 >> 8*1 & 0xFF; c1 = lerp(c1, c11 >> 8*1 & 0xFF, Δx);
+ c = lerp(c0, c1, Δy);
*(a + 1) = c;
case 1:
- c0 = c00 >> 8*0 & 0xFF; c0 = lerp(c0, c01 >> 8*0 & 0xFF, Δx2);
- c1 = c10 >> 8*0 & 0xFF; c1 = lerp(c1, c11 >> 8*0 & 0xFF, Δx2);
- c = lerp(c0, c1, Δy2);
+ c0 = c00 >> 8*0 & 0xFF; c0 = lerp(c0, c01 >> 8*0 & 0xFF, Δx);
+ c1 = c10 >> 8*0 & 0xFF; c1 = lerp(c1, c11 >> 8*0 & 0xFF, Δx);
+ c = lerp(c0, c1, Δy);
*(a + 0) = c;
}
}
--- a/sys/src/libmemdraw/warptest.c
+++ b/sys/src/libmemdraw/warptest.c
@@ -79,6 +79,7 @@
main(int argc, char *argv[])
{Memimage *dst, *src;
+ Warp w;
Rectangle dr;
double sx, sy, tx, ty, θ, c, s;
@@ -130,9 +131,11 @@
dr = src->r;
dst = allocmemimage(dr, src->chan);
+ mkwarp(w, T);
profbegin();
- memaffinewarp(dst, dst->r, src, src->r.min, (uchar*)T);
+ if(memaffinewarp(dst, dst->r, src, src->r.min, w) < 0)
+ sysfatal("memaffinewarp: %r");profend();
writememimage(1, dst);
--
⑨