summaryrefslogtreecommitdiff
path: root/ccast/filt.c
diff options
context:
space:
mode:
Diffstat (limited to 'ccast/filt.c')
-rw-r--r--ccast/filt.c531
1 files changed, 531 insertions, 0 deletions
diff --git a/ccast/filt.c b/ccast/filt.c
new file mode 100644
index 0000000..ab075ae
--- /dev/null
+++ b/ccast/filt.c
@@ -0,0 +1,531 @@
+
+/*
+ * Argyll Color Correction System
+ * ChromCast up filter test code.
+ *
+ * Author: Graeme W. Gill
+ * Date: 28/8/2014
+ *
+ * Copyright 2014 Graeme W. Gill
+ * All rights reserved.
+ *
+ * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
+ * see the License2.txt file for licencing details.
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <sys/types.h>
+#include <time.h>
+#include "copyright.h"
+#include "aconfig.h"
+#ifndef SALONEINSTLIB
+#include "numlib.h"
+#else
+#include "numsup.h"
+#endif
+#include "yajl.h"
+#include "conv.h"
+#include "base64.h"
+#include "ccmdns.h"
+#include "ccpacket.h"
+#include "ccmes.h"
+#include "yajl.h"
+
+#define DO_WEIGHTING
+#define SUM_CONSTRAINT
+
+//#define VERT 1 /* 1 for vertical */
+
+#ifndef DBL_PI
+# define DBL_PI 3.1415926535897932384626433832795
+#endif
+
+double lanczos3(double wi, double x) {
+ double y;
+
+ x = fabs(1.0 * x/wi);
+ if (x >= 3.0)
+ return 0.0;
+ if (x < 1e-6)
+ return 1.0;
+ y = sin(DBL_PI * x)/(DBL_PI * x) * sin(DBL_PI * x/3.0)/(DBL_PI * x/3.0);
+
+ return y;
+}
+
+double lanczos2(double wi, double x) {
+ double y;
+
+ x = fabs(1.0 * x/wi);
+ if (x >= 2.0)
+ return 0.0;
+ if (x < 1e-6)
+ return 1.0;
+ y = sin(DBL_PI * x)/(DBL_PI * x) * sin(DBL_PI * x/2.0)/(DBL_PI * x/2.0);
+
+ return y;
+}
+
+double in[2][72] = {
+{ // Horizontal
+255, 0, 0, 0, 0, 0, 0, 0, 0,
+255, 0, 0, 0, 0, 0, 0, 0, 0,
+254, 0, 0, 0, 0, 0, 0, 0, 0,
+254, 0, 0, 0, 0, 0, 0, 0, 0,
+253, 0, 0, 0, 0, 0, 0, 0, 0,
+253, 0, 0, 0, 0, 0, 0, 0, 0,
+252, 0, 0, 0, 0, 0, 0, 0, 0,
+252, 0, 0, 0, 0, 0, 0, 0, 0
+}, { // Vertical slice input target
+255, 0, 0, 0, 0, 0, 0, 0, 0,
+255, 0, 0, 0, 0, 0, 0, 0, 0,
+254, 0, 0, 0, 0, 0, 0, 0, 0,
+254, 0, 0, 0, 0, 0, 0, 0, 0,
+253, 0, 0, 0, 0, 0, 0, 0, 0,
+253, 0, 0, 0, 0, 0, 0, 0, 0,
+252, 0, 0, 0, 0, 0, 0, 0, 0,
+252, 0, 0, 0, 0, 0, 0, 0, 0
+} };
+
+#ifdef NEVER
+double out[2][] = {
+{ // Horizontal
+170, 110, 23, 5, 15, 17, 16, 16, 16, 17, 10, 6, 57,
+151, 153, 61, 7, 10, 16, 16, 16, 16, 17, 15, 5, 22,
+107, 169, 109, 23, 5, 15, 17, 16, 16, 16, 17, 10, 6, 57,
+150, 152, 61, 7, 10, 16, 16, 16, 16, 17, 15, 5, 22,
+107, 168, 109, 23, 5, 15, 17, 16, 16, 16, 17, 10, 6, 56,
+149, 151, 61, 7, 10, 16, 16, 16, 16, 17, 15, 6, 22,
+106, 167, 108, 23, 5, 15, 17, 16, 16, 16, 17, 10, 6, 56,
+148, 150
+},
+{ // Vertical slice target output
+235, 100, 0, 16, 16, 16, 16, 16, 16, 16, 20, 1, 13,
+191, 194, 19, 0, 20, 16, 16, 16, 16, 16, 16, 16, 0, 95,
+234, 99, 0, 16, 16, 16, 16, 16, 16, 16, 20, 1, 13,
+190, 194, 19, 0, 20, 16, 16, 16, 16, 16, 16, 16, 0, 94,
+233, 99, 0, 16, 16, 16, 16, 16, 16, 16, 20, 1, 13,
+189, 193, 19, 0, 20, 16, 16, 16, 16, 16, 16, 16, 0, 94,
+232, 99, 0, 16, 16, 16, 16, 16, 16, 16, 20, 1, 13,
+188, 192
+}
+};
+#else
+//#define N1 -9
+//#define N2 -8
+#define N1 0
+#define N2 0
+
+//#define N3 1
+#define N4 0
+
+double out[2][96] = {
+{ // Horozontal
+170, 110, 23, 5, 15, 17, 16, 16, 16, 17, 10, 6, 57,
+151, 153, 61, 7, 10, 16, 16, 16, 16, 17, 15, 5, 22,
+107, 169, 109, 23, 5, 15, 17, 16, 16, 16, 17, 10, 6, 57,
+150, 152, 61, 7, 10, 16, 16, 16, 16, 17, 15, 5, 22,
+107, 168, 109, 23, 5, 15, 17, 16, 16, 16, 17, 10, 6, 56,
+149, 151, 61, 7, 10, 16, 16, 16, 16, 17, 15, 6, 22,
+106, 167, 108, 23, 5, 15, 17, 16, 16, 16, 17, 10, 6, 56,
+148, 150
+},
+{ // Vertical slice target output
+235, 100, N2, 16, 16, 16, 16, 16, 16, 16, 20, 1, 13,
+191, 194, 19, N4, 20, 16, 16, 16, 16, 16, 16, 16, N1, 95,
+234, 99, N2, 16, 16, 16, 16, 16, 16, 16, 20, 1, 13,
+190, 194, 19, N4, 20, 16, 16, 16, 16, 16, 16, 16, N1, 94,
+233, 99, N2, 16, 16, 16, 16, 16, 16, 16, 20, 1, 13,
+189, 193, 19, N4, 20, 16, 16, 16, 16, 16, 16, 16, N1, 94,
+232, 99, N2, 16, 16, 16, 16, 16, 16, 16, 20, 1, 13,
+188, 192
+}
+};
+#endif
+
+// Computed the input to output filters.
+// There will be two phases, depending on whether
+// the input pixel has an even or odd address.
+// Although in this case they seem like they
+// almost interleave, they are actually mirror
+// images with offset, so it's easier to keep them
+// separate, rather than trying to figure the offset out.
+// The index is the output pixel around the closest one
+// to the scaled input pixel - ie. floor(in * 1.5 + 0.5)
+#define FWIDTH 6
+#define NWIDTH (2 * FWIDTH + 1)
+
+double filt_v[2][2][NWIDTH];
+double filt_vx[2][2][NWIDTH]; /* max */
+double filt_vn[2][2][NWIDTH]; /* min */
+
+double *filt[2][2] = { { &filt_v[0][0][FWIDTH], &filt_v[0][1][FWIDTH] },
+ { &filt_v[1][0][FWIDTH], &filt_v[1][1][FWIDTH] } };
+
+double *filtx[2][2] = { { &filt_vx[0][0][FWIDTH], &filt_vx[0][1][FWIDTH] },
+ { &filt_vx[1][0][FWIDTH], &filt_vx[1][1][FWIDTH] } };
+
+double *filtn[2][2] = { { &filt_vn[0][0][FWIDTH], &filt_vn[0][1][FWIDTH] },
+ { &filt_vn[1][0][FWIDTH], &filt_vn[1][1][FWIDTH] } };
+
+//int fneg[2] = { -5, -4 }; /* Negative index range (inclusive) */
+//int fpos[2] = { 5, 3 }; /* Positive index range (inclusive) */
+int fneg[2][2] = { { -4, -4 }, /* Negative index range (inclusive) */
+ { -4, -4 } };
+int fpos[2][2] = { { 4, 3 }, /* Positive index range (inclusive) */
+ { 4, 3 } };
+
+/* Weightings [horiz/vert][phase] */
+double filtw_v[2][2][NWIDTH] = {
+ /* -6, -5, -4, -3, -2, -1, 0, +1, +2, +3, +4, +5, +6 */
+// { { 1.0, 1.0, 1.0, 5.0, 3.0, 42.0, 80.0, 43.0, 3.0, 5.0, 1.0, 1.0, 1.0 },
+// { 1.0, 1.0, 3.0, 4.0, 19.0, 62.0, 62.0, 21.0, 4.0, 3.0, 1.0, 1.0, 1.0 } },
+ { { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 },
+ { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 } },
+ /* -6, -5, -4, -3, -2, -1, 0, +1, +2, +3, +4, +5, +6 */
+// { { 1.0, 1.0, 1.0, 3.0, 17.0, 36.0, 100.0, 40.0, 15.0, 1.0, 3.0, 1.0, 1.0 },
+// { 1.0, 1.0, 1.0, 8.0, 2.0, 80.0, 81.0, 2.0, 8.0, 2.0, 1.0, 1.0, 1.0 } }
+ { { 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 4.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0 },
+ { 1.0, 1.0, 1.0, 1.0, 1.0, 3.0, 3.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 } },
+};
+double *filtw[2][2] = { { &filtw_v[0][0][FWIDTH], &filtw_v[0][1][FWIDTH] },
+ { &filtw_v[1][0][FWIDTH], &filtw_v[1][1][FWIDTH] } };
+
+#define MX 2.0
+#define MN -0.5
+
+// Compute the filter shapes
+// vv = 0 for horizontal, 1 for vertical
+static void compute(int vv) {
+ int ph, ii, i, jj, j;
+ double iv, ov;
+ int fcount[2];
+ int niv = sizeof(in[vv])/sizeof(double);
+ int nov = sizeof(out[vv])/sizeof(double);
+
+ // Clear the filters
+ for (ph = 0; ph < 2; ph++) {
+ for (j = -FWIDTH; j <= FWIDTH; j++) {
+ filt[vv][ph][j] = 0.0;
+ filtx[vv][ph][j] = MX;
+ filtn[vv][ph][j] = MN;
+ }
+ fcount[ph] = 0;
+ }
+
+ // Discover an input value
+ for (ii = 0; ii < niv; ii++) {
+ double rgb[3], ycc[3];
+ double prop, propx, propn;
+
+ if (in[ii] == 0)
+ continue;
+
+ iv = in[vv][ii];
+ rgb[0] = rgb[1] = rgb[2] = iv;
+ ccast2YCbCr(NULL, ycc, rgb);
+ iv = ycc[0];
+
+ ph = ii & 1;
+ jj = (int)floor(ii * 1.5 + 0.5);
+
+ for (j = -FWIDTH; j <= FWIDTH; j++) {
+ int k = jj + j;
+ if (k < 0 || k >= nov)
+ continue;
+ ov = out[vv][k];
+
+ prop = ((ov - 16.0)/219.0)/((iv - 16.0)/219.0);
+ propx = ((ov - 16.0)/219.0)/((iv - 0.5 - 16.0)/219.0);
+ propn = ((ov - 16.0)/219.0)/((iv + 0.5 - 16.0)/219.0);
+
+ if (propx < propn) {
+ double tt = propn;
+ propn = propx;
+ propx = tt;
+ }
+
+//printf("~1 phase %d, off %d, iv %f ov %f, prop %f\n",ph,j,iv,ov,prop);
+ filt[vv][ph][j] += prop;
+ if (propx < filtx[vv][ph][j])
+ filtx[vv][ph][j] = propx;
+ if (propn > filtn[vv][ph][j])
+ filtn[vv][ph][j] = propn;
+ }
+ fcount[ph]++;
+ }
+
+ // Compute average values
+ for (ph = 0; ph < 2; ph++) {
+ for (j = -FWIDTH; j <= FWIDTH; j++)
+ filt[vv][ph][j] /= (double)fcount[ph];
+ }
+}
+
+#define FCO2IX(bank, off) (bank ? fpos[vv][0] - fneg[vv][0] + 1 + off - fneg[vv][1] : off - fneg[vv][0])
+
+// Compute the filter shapes using SVD
+// vv = 0 for horizontal, 1 for vertical
+static void compute2(int vv) {
+ int niv = sizeof(in[vv])/sizeof(double); /* Number of input values */
+ int nov = sizeof(out[vv])/sizeof(double); /* Number of output values */
+ int novextra = 0;
+ int nfc = fpos[vv][0] - fneg[vv][0] + 1
+ + fpos[vv][1] - fneg[vv][1] + 1; /* Number of filter coeficients */
+ double **A, *b;
+ int oe; /* Even or odd */
+ int j, i;
+
+#ifdef SUM_CONSTRAINT
+ novextra = 3;
+#endif
+
+ nov += novextra; /* Extra constraint of sum */
+
+ /* We assume nov > nfc */
+ A = dmatrixz(0, nov-1, 0, nfc-1);
+ b = dvectorz(0, nov-1);
+
+ /* For each output value */
+ for (j = 0; j < (nov-novextra); j++) {
+ double ww = 1.0;
+
+#ifdef DO_WEIGHTING
+ /* Figure out the weighting */
+ for (oe = 0; oe < 2; oe++) {
+ int fix; /* Filter index */
+
+ /* For offset range of filter */
+ for (fix = fneg[vv][oe]; fix <= fpos[vv][oe]; fix++) {
+ int ocx, icx, ph;
+
+ ocx = j - fix; /* Output center of filter */
+ if (ocx < 0 || ocx >= nov)
+ continue; /* Filter would never get applied */
+ if (((2 * ocx) % 3) == 2)
+ continue; /* Would never get applied */
+ icx = (int)floor(ocx / 1.5); /* Input center index for this output */
+ if (icx < 0 || icx >= niv)
+ continue; /* Filter would never get applied */
+ ph = icx & 1; /* Phase of filter */
+ if (ph != oe) /* Not a filter that would appear at this ouput */
+ continue;
+ if (in[vv][icx] >= 200.0)
+ ww = filtw[vv][ph][fix];
+ }
+ }
+#endif /* DO_WEIGHTING */
+
+ /* For even and odd filters */
+ for (oe = 0; oe < 2; oe++) {
+ int fix; /* Filter index */
+
+ /* For offset range of filter */
+ for (fix = fneg[vv][oe]; fix <= fpos[vv][oe]; fix++) {
+ double rgb[3], ycc[3];
+ int ocx, icx, ph;
+
+ ocx = j - fix; /* Output center of filter */
+ if (ocx < 0 || ocx >= nov)
+ continue; /* Filter would never get applied */
+ if (((2 * ocx) % 3) == 2)
+ continue; /* Would never get applied */
+ icx = (int)floor(ocx / 1.5); /* Input center index for this output */
+ if (icx < 0 || icx >= niv)
+ continue; /* Filter would never get applied */
+ ph = icx & 1; /* Phase of filter */
+ if (ph != oe) /* Not a filter that would appear at this ouput */
+ continue;
+//printf("j = %d/%d, k = %d/%d, ix = %d/%d\n",j,nov,oe * NWIDTH + FWIDTH + fix,nfc,ix,niv);
+ rgb[0] = rgb[1] = rgb[2] = in[vv][icx];
+ ccast2YCbCr(NULL, ycc, rgb);
+ A[j][FCO2IX(oe, fix)] += ww * ycc[0];
+//printf("A[%d][%d] = %f\n",j,FCO2IX(oe, fix),A[j][FCO2IX(oe, fix)]);
+ }
+ }
+ b[j] = ww * out[vv][j];
+//printf("b[%d] = %f\n",j,b[j]);
+ }
+
+#ifdef SUM_CONSTRAINT
+ /* Add sum constraints */
+ /* For 3 repeating output slots */
+ for (j = nov-novextra; j < nov; j++) {
+ double ww = 10000.0;
+ int jj = j - (nov-novextra);
+
+ b[j] = ww;
+
+ /* For even and odd filters */
+ for (oe = 0; oe < 2; oe++) {
+ int fix; /* Filter index */
+
+ /* For offset range of filter */
+ for (fix = fneg[vv][oe]; fix <= fpos[vv][oe]; fix++) {
+ double rgb[3], ycc[3];
+ int ocx, ocx2, icx, ph;
+
+ ocx = j - fix; /* Output center of filter */
+ ocx2 = 2 * ocx;
+
+ while (ocx2 < 0)
+ ocx2 += 3;
+ while (ocx2 >= 3)
+ ocx2 -= 3;
+
+ if (ocx2 == 2)
+ continue; /* Would never get applied */
+
+ while (ocx < 0)
+ ocx += 3;
+ while (ocx >= 3)
+ ocx -= 3;
+
+ icx = (int)floor(ocx / 1.5); /* Input center index for this output */
+ ph = icx & 1; /* Phase of filter */
+ if (ph != oe) /* Not a filter that would appear at this ouput */
+ continue;
+ A[j][FCO2IX(oe, fix)] = ww;
+printf("A[%d][%d] = %f\n",j,FCO2IX(oe, fix),A[j][FCO2IX(oe, fix)]);
+ }
+ }
+ }
+#endif /* SUM_CONSTRAINT */
+
+ /* Solve the equation A.x = b using SVD */
+ /* (The w[] values are thresholded for best accuracy) */
+ /* Return non-zero if no solution found */
+ if (svdsolve(A, b, nov, nfc))
+ error("svdsolve failed");
+
+ /* Print the filter shape */
+ /* and copy to the filter */
+ printf("SVD computed for %s:\n", vv ? "vertical" : "horizontal");
+ for (oe = 0; oe < 2; oe++) {
+ int fix; /* Filter index */
+ double sum = 0.0;
+
+ printf("Phase %d\n",oe);
+// for (fix = -FWIDTH; fix <= FWIDTH; fix++) {
+ for (fix = fneg[vv][oe]; fix <= fpos[vv][oe]; fix++) {
+ printf(" %d -> %f\n",fix, b[FCO2IX(oe, fix)]);
+ sum += b[FCO2IX(oe, fix)];
+ filt[vv][oe][fix] = b[FCO2IX(oe, fix)];
+ }
+ printf("sum = %f\n",sum);
+ }
+}
+
+void check(int vv) {
+ double *chout;
+ int niv = sizeof(in[vv])/sizeof(double); /* Number of input values */
+ int nov = sizeof(out[vv])/sizeof(double); /* Number of output values */
+
+ int range, i, ii, j;
+ double xc, x, iv, tw, w, y;
+ double cout, terr = 0.0;
+
+printf("~1 nov = %d\n",nov);
+ if ((chout = (double *)malloc(sizeof(double) * nov)) == NULL)
+ error("Malloc failed");
+
+ // Clear the output
+ for (i = 0; i < nov; i++)
+ chout[i] = 0.0;
+
+ // For all the input value
+ for (ii = 0; ii < niv; ii++) {
+ int ph, jj;
+ double rgb[3], ycc[3];
+ double prop;
+
+ iv = in[vv][ii];
+ rgb[0] = rgb[1] = rgb[2] = iv;
+ ccast2YCbCr(NULL, ycc, rgb);
+ iv = ycc[0];
+
+ ph = ii & 1;
+ jj = (int)floor(ii * 1.5 + 0.5);
+
+ for (j = -FWIDTH; j <= FWIDTH; j++) {
+ int k = jj + j;
+ if (k < 0 || k >= nov)
+ continue;
+
+//if ((jj + j) == 4) printf("[%d] += w %f * iv %f\n",jj + j, filt[ph][j], iv);
+ chout[k] += filt[vv][ph][j] * iv;
+ }
+ }
+
+ for (i = 0; i < nov; i++) {
+ double ov, ee;
+ ov = chout[i];
+ ov = floor(ov + 0.5);
+#ifdef NEVER
+ if (ov < 0.0)
+ ov = 0.0;
+ else if (ov > 255.0)
+ ov = 255.0;
+#endif
+ ee = ov - out[vv][i];
+ terr += ee * ee;
+ printf("out %d = %f should be %f err %f\n",i,chout[i],out[vv][i],ee);
+ }
+
+ printf("Total err = %f RMS\n",sqrt(terr));
+}
+
+int
+main(int argc,
+ char *argv[]
+) {
+ int ph, vv, j;
+
+ double err;
+ double cp[2]; /* Initial starting point */
+ double s[2]; /* Size of initial search area */
+
+ printf("Hi there\n");
+
+#ifdef NEVER
+ compute(1);
+
+ // Print filter shape
+ printf("Directly computed:\n");
+ for (ph = 0; ph < 2; ph++) {
+ printf("Phase %d\n",ph);
+ for (j = -FWIDTH; j <= FWIDTH; j++)
+ printf(" %d -> min %f, avg %f, max %f\n",j,filtn[1][ph][j],filt[1][ph][j],filtx[1][ph][j]);
+ }
+#endif
+
+ for (vv = 0; vv < 2; vv++) {
+ compute2(vv);
+ check(vv);
+ }
+
+ /* Output code to stdout */
+ for (ph = 0; ph < 2; ph++) {
+ fprintf(stderr,"/* Weightings [horiz/vert] */\n");
+ fprintf(stderr,"double filt_v_%s[2][%s_WIDTH] = {\n",ph ? "od" : "ev", ph ? "OD" : "EV");
+
+ for (vv = 0; vv < 2; vv++) {
+ int fix; /* Filter index */
+
+ fprintf(stderr,"{ ");
+ for (fix = fneg[vv][ph]; fix <= fpos[vv][ph]; fix++) {
+ if (fix > fneg[vv][ph])
+ fprintf(stderr,", ");
+ fprintf(stderr,"%f", filt[vv][ph][fix]);
+ }
+ fprintf(stderr," }%s\n",vv == 0 ? "," : "");
+ }
+ fprintf(stderr,"};\n\n");
+ }
+
+ return 0;
+}
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-19 06:29:01 +0000