cosmopolitan/dsp/scale/gyarados.c

282 lines
9.7 KiB
C

/*-*- mode:c;indent-tabs-mode:nil;c-basic-offset:2;tab-width:8;coding:utf-8 -*-│
│vi: set net ft=c ts=2 sts=2 sw=2 fenc=utf-8 :vi│
╞══════════════════════════════════════════════════════════════════════════════╡
│ Copyright 2020 Justine Alexandra Roberts Tunney │
│ │
│ This program is free software; you can redistribute it and/or modify │
│ it under the terms of the GNU General Public License as published by │
│ the Free Software Foundation; version 2 of the License. │
│ │
│ This program is distributed in the hope that it will be useful, but │
│ WITHOUT ANY WARRANTY; without even the implied warranty of │
│ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU │
│ General Public License for more details. │
│ │
│ You should have received a copy of the GNU General Public License │
│ along with this program; if not, write to the Free Software │
│ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA │
│ 02110-1301 USA │
╚─────────────────────────────────────────────────────────────────────────────*/
#include "dsp/core/c161.h"
#include "dsp/core/core.h"
#include "dsp/core/ituround.h"
#include "dsp/core/q.h"
#include "dsp/core/twixt8.h"
#include "libc/limits.h"
#include "libc/log/check.h"
#include "libc/log/log.h"
#include "libc/macros.h"
#include "libc/math.h"
#include "libc/mem/mem.h"
#include "libc/nexgen32e/bsr.h"
#include "libc/runtime/gc.h"
#include "libc/str/str.h"
#include "libc/testlib/testlib.h"
#include "libc/x/x.h"
#include "third_party/dtoa/dtoa.h"
#include "tool/viz/lib/knobs.h"
/**
* @fileoverview Gyarados resizes graphics.
* @note H/T John Costella, Facebook, Photoshop, Carl Friedrich Gauss
* @note Eric Brasseur has an interesting blog post on tip of iceberg
* @see Magikarp
*/
#define M 14
#define SQR(X) ((X) * (X))
struct SamplingSolution {
int n, s;
void *weights;
void *indices;
};
static double ComputeWeight(double x) {
if (-1.5 < x && x < 1.5) {
if (-.5 < x && x < .5) {
return .75 - SQR(x);
} else if (x < 0) {
return .5 * SQR(x + 1.5);
} else {
return .5 * SQR(x - 1.5);
}
} else {
return 0;
}
}
static struct SamplingSolution *NewSamplingSolution(long n, long s) {
struct SamplingSolution *ss;
ss = xcalloc(1, sizeof(struct SamplingSolution));
ss->n = n;
ss->s = s;
ss->weights = xcalloc(n * s, sizeof(short));
ss->indices = xcalloc(n * s, sizeof(short));
return ss;
}
static bool IsNormalized(int n, double A[n]) {
int i;
double x;
for (x = i = 0; i < n; ++i) x += A[i];
return fabs(x - 1) < 1e-4;
}
void FreeSamplingSolution(struct SamplingSolution *ss) {
long i;
if (ss) {
free(ss->indices);
free(ss->weights);
free(ss);
}
}
struct SamplingSolution *ComputeSamplingSolution(long dn, long sn, double dar,
double off, double par) {
double *fweights;
double sum, hw, w, x, f;
short *weights, *indices;
struct SamplingSolution *res;
long j, i, k, n, min, max, s, N[6];
if (!dar) dar = sn, dar /= dn;
if (!off) off = (dar - 1) / 2;
f = dar < 1 ? 1 / dar : dar;
s = 3 * f + 4;
fweights = gc(xcalloc(s, sizeof(double)));
res = NewSamplingSolution(dn, s);
weights = res->weights;
indices = res->indices;
for (i = 0; i < dn; ++i) {
x = off + i * dar;
hw = 1.5 * f;
min = ceil(x - hw);
max = floor(x + hw);
n = max - min + 1;
CHECK_LE(n, s);
for (k = 0, j = min; j <= max; ++j) {
fweights[k++] = ComputeWeight((j - x) / (f / par));
}
for (sum = k = 0; k < n; ++k) sum += fweights[k];
for (j = 0; j < n; ++j) fweights[j] *= 1 / sum;
DCHECK(IsNormalized(n, fweights));
for (j = 0; j < n; ++j) {
indices[i * s + j] = MIN(sn - 1, MAX(0, min + j));
}
for (j = 0; j < n; j += 6) {
GetIntegerCoefficients(N, fweights + j, M, 0, 255);
for (k = 0; k < MIN(6, n - j); ++k) {
weights[i * s + j + k] = N[k];
}
}
}
return res;
}
static void *ZeroMatrix(long yw, long xw, int p[yw][xw], long yn, long xn) {
long y;
for (y = 0; y < yn; ++y) {
memset(p[y], 0, xn);
}
return p;
}
static int Sharpen(int ax, int bx, int cx) {
return (-1 * ax + 6 * bx + -1 * cx + 2) / 4;
}
static void GyaradosImpl(long dyw, long dxw, int dst[dyw][dxw], long syw,
long sxw, const int src[syw][sxw], long dyn, long dxn,
long syn, long sxn, int tmp0[restrict dyn][sxn],
int tmp1[restrict dyn][sxn],
int tmp2[restrict dyn][dxn], long yfn, long xfn,
const short fyi[dyn][yfn], const short fyw[dyn][yfn],
const short fxi[dxn][xfn], const short fxw[dxn][xfn],
bool sharpen) {
long i, j;
int eax, dy, dx, sy, sx;
for (sx = 0; sx < sxn; ++sx) {
for (dy = 0; dy < dyn; ++dy) {
for (eax = i = 0; i < yfn; ++i) {
eax += fyw[dy][i] * src[fyi[dy][i]][sx];
}
tmp0[dy][sx] = QRS(M, eax);
}
}
for (dy = 0; dy < dyn; ++dy) {
for (sx = 0; sx < sxn; ++sx) {
tmp1[dy][sx] = sharpen ? Sharpen(tmp0[MIN(dyn - 1, MAX(0, dy - 1))][sx],
tmp0[dy][sx],
tmp0[MIN(dyn - 1, MAX(0, dy + 1))][sx])
: tmp0[dy][sx];
}
}
for (dx = 0; dx < dxn; ++dx) {
for (dy = 0; dy < dyn; ++dy) {
for (eax = i = 0; i < xfn; ++i) {
eax += fxw[dx][i] * tmp1[dy][fxi[dx][i]];
}
tmp2[dy][dx] = QRS(M, eax);
}
}
for (dx = 0; dx < dxn; ++dx) {
for (dy = 0; dy < dyn; ++dy) {
dst[dy][dx] = sharpen ? Sharpen(tmp2[dy][MIN(dxn - 1, MAX(0, dx - 1))],
tmp2[dy][dx],
tmp2[dy][MIN(dxn - 1, MAX(0, dx + 1))])
: tmp2[dy][dx];
}
}
}
/**
* Scales image.
*
* @note gyarados is magikarp in its infinite form
* @see Magikarp2xY(), Magikarp2xX()
*/
void *Gyarados(long dyw, long dxw, int dst[dyw][dxw], long syw, long sxw,
const int src[syw][sxw], long dyn, long dxn, long syn, long sxn,
struct SamplingSolution *cy, struct SamplingSolution *cx,
bool sharpen) {
if (dyn > 0 && dxn > 0) {
if (syn > 0 && sxn > 0) {
CHECK_LE(syn, syw);
CHECK_LE(sxn, sxw);
CHECK_LE(dyn, dyw);
CHECK_LE(dxn, dxw);
CHECK_LT(bsrl(syn) + bsrl(sxn), 32);
CHECK_LT(bsrl(dyn) + bsrl(dxn), 32);
CHECK_LE(dyw, 0x7fff);
CHECK_LE(dxw, 0x7fff);
CHECK_LE(syw, 0x7fff);
CHECK_LE(sxw, 0x7fff);
CHECK_LE(dyn, 0x7fff);
CHECK_LE(dxn, 0x7fff);
CHECK_LE(syn, 0x7fff);
CHECK_LE(sxn, 0x7fff);
GyaradosImpl(dyw, dxw, dst, syw, sxw, src, dyn, dxn, syn, sxn,
gc(xmemalign(64, sizeof(int) * dyn * sxn)),
gc(xmemalign(64, sizeof(int) * dyn * sxn)),
gc(xmemalign(64, sizeof(int) * dyn * dxn)), cy->s, cx->s,
cy->indices, cy->weights, cx->indices, cx->weights, sharpen);
} else {
ZeroMatrix(dyw, dxw, dst, dyn, dxn);
}
}
return dst;
}
void *GyaradosUint8(long dyw, long dxw, unsigned char dst[dyw][dxw], long syw,
long sxw, const unsigned char src[syw][sxw], long dyn,
long dxn, long syn, long sxn, long lo, long hi,
struct SamplingSolution *cy, struct SamplingSolution *cx,
bool sharpen) {
static bool once;
static int Tin[256];
static unsigned char Tout[32768];
long i, y, x;
int(*tmp)[MAX(dyn, syn)][MAX(dxn, sxn)];
if (!once) {
for (i = 0; i < ARRAYLEN(Tin); ++i) {
Tin[i] = F2Q(15, rgb2linpc(i / 255., 2.4));
}
for (i = 0; i < ARRAYLEN(Tout); ++i) {
Tout[i] = MIN(255, MAX(0, round(rgb2stdpc(Q2F(15, i), 2.4) * 255.)));
}
once = true;
}
tmp = xmemalign(64, sizeof(int) * MAX(dyn, syn) * MAX(dxn, sxn));
for (y = 0; y < syn; ++y) {
for (x = 0; x < sxn; ++x) {
(*tmp)[y][x] = Tin[src[y][x]];
}
}
Gyarados(MAX(dyn, syn), MAX(dxn, sxn), *tmp, MAX(dyn, syn), MAX(dxn, sxn),
*tmp, dyn, dxn, syn, sxn, cy, cx, sharpen);
for (y = 0; y < dyn; ++y) {
for (x = 0; x < dxn; ++x) {
dst[y][x] = Tout[MIN(32767, MAX(0, (*tmp)[y][x]))];
}
}
free(tmp);
return dst;
}
void *EzGyarados(long dcw, long dyw, long dxw, unsigned char dst[dcw][dyw][dxw],
long scw, long syw, long sxw,
const unsigned char src[scw][syw][sxw], long c0, long cn,
long dyn, long dxn, long syn, long sxn, double ry, double rx,
double oy, double ox) {
long c;
struct SamplingSolution *cy, *cx;
cy = ComputeSamplingSolution(dyn, syn, ry, oy, 1);
cx = ComputeSamplingSolution(dxn, sxn, rx, ox, 1);
for (c = c0; c < cn; ++c) {
GyaradosUint8(dyw, dxw, dst[c], syw, sxw, src[c], dyn, dxn, syn, sxn, 0,
255, cy, cx, true);
}
FreeSamplingSolution(cx);
FreeSamplingSolution(cy);
return dst;
}