From 094535c010320967639e8e86f974d878e80baa72 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B6rg=20Frings-F=C3=BCrst?= Date: Fri, 1 May 2015 16:13:57 +0200 Subject: Imported Upstream version 1.7.0 --- ccast/filt.c | 531 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 531 insertions(+) create mode 100644 ccast/filt.c (limited to 'ccast/filt.c') 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 +#include +#include +#include +#include +#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; +} + -- cgit v1.2.3