283 lines
		
	
	
		
			9.8 KiB
		
	
	
	
		
			C
		
	
	
			
		
		
	
	
			283 lines
		
	
	
		
			9.8 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      15
 | |
| #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 + 1;
 | |
|   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, short tmp0[restrict dyn][sxn],
 | |
|                          short tmp1[restrict dyn][sxn],
 | |
|                          short 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(short) * dyn * sxn)),
 | |
|                    gc(xmemalign(64, sizeof(short) * dyn * sxn)),
 | |
|                    gc(xmemalign(64, sizeof(short) * 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;
 | |
| }
 |