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/dpat.c | 952 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 952 insertions(+) create mode 100644 ccast/dpat.c (limited to 'ccast/dpat.c') diff --git a/ccast/dpat.c b/ccast/dpat.c new file mode 100644 index 0000000..4962ae1 --- /dev/null +++ b/ccast/dpat.c @@ -0,0 +1,952 @@ + +/* + * Argyll Color Correction System + * ChromCast dither pattern code + * + * Author: Graeme W. Gill + * Date: 11/12/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. + * + */ + +/* + * Locates minimal list of surrounders + * Find initial weight for them using minimizer + Itterates over + * Locate cell changes that match desired correction + with increasing randomness of choosing a poor match + * Reset to starting pattern when accumulated error over + best current is exceeded. + + + Problem is this is unreliable is locating good matches + for some cases. + + Problem is that it doesn't seem to work - ChromeCast doesn't + decode the way we model it :-( + + */ + +#include +#include +#include +#include +#include +#include "copyright.h" +#include "aconfig.h" +#include "counters.h" +#ifndef SALONEINSTLIB +#include "numlib.h" +#else +#include "numsup.h" +#endif +#include "conv.h" +#include "base64.h" +#include "yajl.h" +#include "ccmdns.h" +#include "ccpacket.h" +#include "ccmes.h" +#include "ccast.h" + +#ifdef STANDALONE_TEST +# define DIAGNOSTICS +#endif + +#undef TEST_POINT + +/* Even/Odd filter filter index ranges (inclusive) and total size */ +#define EV_NEG -4 +#define EV_POS 4 +#define EV_WIDTH (-(EV_NEG) + 1 + EV_POS) +#define OD_NEG -4 +#define OD_POS 3 +#define OD_WIDTH (-(OD_NEG) + 1 + OD_POS) + +/* Weightings [horiz/vert] */ +double filt_v_ev[2][EV_WIDTH] = { +{ -0.003467, -0.048341, 0.028688, 0.418873, 0.704674, 0.427551, 0.032480, -0.050146, -0.003793 }, +{ -0.013477, 0.000081, -0.087119, 0.357627, 1.000252, 0.378473, -0.083496, -0.000155, -0.009850 } +}; + +/* Weightings [horiz/vert] */ +double filt_v_od[2][OD_WIDTH] = { +{ -0.026810, -0.045190, 0.186825, 0.614579, 0.623721, 0.206987, -0.040217, -0.026419 }, +{ 0.008396, -0.078987, -0.013733, 0.796594, 0.813578, 0.013555, -0.086466, 0.004781 } +}; + + +/* [horiz/vert][phase] */ +double *filt[2][2] = { { &filt_v_ev[0][-(EV_NEG)], &filt_v_od[0][-(OD_NEG)] }, + { &filt_v_ev[1][-(EV_NEG)], &filt_v_od[1][-(OD_NEG)] } }; + +/* [phase] */ +int fneg[2] = { EV_NEG, OD_NEG }; +int fpos[2] = { EV_POS, OD_POS }; + +/* Dither input size */ +#ifdef TEST_POINT +# define DISIZE 8 +#else +# define DISIZE CCDITHSIZE +#endif +#define DOSIZE ((DISIZE * 3)/2) + +#ifdef DIAGNOSTICS +int vv = 0; +#endif + +/* Upsample ipat to opat */ +/* And return the average RGB value */ +static void upsample(double ret[3], double iipat[DISIZE][DISIZE][3]) { + double ipat[DISIZE][DISIZE][3]; + double tpat[DOSIZE][DISIZE][3]; /* Intermediate pattern - horizontal up-sampled */ + double opat[DOSIZE][DOSIZE][3]; /* Output pattern */ + int x, y; + int i, j; + +//printf("~1 doing conversion\n"); + /* Convert from RGB to YCbCr */ + for (x = 0; x < DISIZE; x++) { + for (y = 0; y < DISIZE; y++) { + ccast2YCbCr(NULL, ipat[x][y], iipat[x][y]); +// ipat[x][y][0] = iipat[x][y][0]; +// ipat[x][y][1] = iipat[x][y][1]; +// ipat[x][y][2] = iipat[x][y][2]; +//printf("[%d][%d] %f %f %f -> %f %f %f\n",x,y, iipat[x][y][0], iipat[x][y][1], iipat[x][y][2], ipat[x][y][0], ipat[x][y][1], ipat[x][y][2]); + } + } + +//printf("~1 doing horiz\n"); + /* Up-sample in horizontal direction */ + for (y = 0; y < DISIZE; y++) { + + /* Zero the intermediate pattern */ + for (i = 0; i < DOSIZE; i++) { + tpat[i][y][0] = 0.0; + tpat[i][y][1] = 0.0; + tpat[i][y][2] = 0.0; + } + + /* Distribute the input according to the filter */ + for (x = 0; x < DISIZE; x++) { + int ph, ii, k; + + ph = x & 1; + ii = (int)floor(x * 1.5 + 0.5); + + /* For all the filter cooeficients */ + for (k = fneg[ph]; k <= fpos[ph]; k++) { + i = (ii + k); + while (i < 0) + i += DOSIZE; + while (i >= DOSIZE) + i -= DOSIZE; + tpat[i][y][0] += filt[0][ph][k] * ipat[x][y][0]; + tpat[i][y][1] += filt[0][ph][k] * ipat[x][y][1]; + tpat[i][y][2] += filt[0][ph][k] * ipat[x][y][2]; + } + } + } + +//printf("~1 doing vert\n"); + /* Up-sample in vertical direction */ + for (i = 0; i < DOSIZE; i++) { + + /* Zero the output pattern */ + for (j = 0; j < DOSIZE; j++) { + opat[i][j][0] = 0.0; + opat[i][j][1] = 0.0; + opat[i][j][2] = 0.0; + } + + /* Distribute the input according to the filter */ + for (y = 0; y < DISIZE; y++) { + int ph, jj, k; + + ph = y & 1; + jj = (int)floor(y * 1.5 + 0.5); + + /* For all the filter cooeficients */ + for (k = fneg[ph]; k <= fpos[ph]; k++) { + j = (jj + k); + while (j < 0) + j += DOSIZE; + while (j >= DOSIZE) + j -= DOSIZE; + opat[i][j][0] += filt[1][ph][k] * tpat[i][y][0]; + opat[i][j][1] += filt[1][ph][k] * tpat[i][y][1]; + opat[i][j][2] += filt[1][ph][k] * tpat[i][y][2]; + } + } + } +//printf("~1 doing conversion\n"); + + /* Convert from YCbCr to RGB */ + for (i = 0; i < DOSIZE; i++) { + for (j = 0; j < DOSIZE; j++) { + YCbCr2ccast(NULL, opat[i][j], opat[i][j]); + } + } +//printf("~1 done\n"); + + if (ret != NULL) { + ret[0] = ret[1] = ret[2] = 0.0; + for (j = 0; j < DOSIZE; j++) { + for (i = 0; i < DOSIZE; i++) { + ret[0] += opat[i][j][0]; + ret[1] += opat[i][j][1]; + ret[2] += opat[i][j][2]; + } + } + ret[0] /= (double)(DOSIZE * DOSIZE); + ret[1] /= (double)(DOSIZE * DOSIZE); + ret[2] /= (double)(DOSIZE * DOSIZE); + } + +#ifdef DIAGNOSTICS + if (vv) { + printf("Result\n"); + for (j = 0; j < DOSIZE; j++) { + for (i = 0; i < DOSIZE; i++) { + if (i > 0) + printf(", "); + printf("% 5.1f % 5.1f % 5.1f",opat[i][j][0],opat[i][j][1], opat[i][j][2]); + } + printf("\n"); + } + } +#endif +} + +/* ------------------------------------------------------------- */ + +/* Given a quantized RGB target, return a quantized RGB that either exactly */ +/* maps to it through YCbCr conversion, or is the closest in the */ +/* direction away from the target value, or the closest clipped value */ +static void quant_rgb(int n, double out[3], double rgb[3]) { + double ycc[3], base[3]; + double tmp[3], chval[3]; + double dist, bdist = 1e6; + double brgb[3], borgb[3]; + int ix, k; + +//printf("Quant RGB %f %f %f n %d\n", rgb[0], rgb[1], rgb[2], n); + + /* Search current rgb and surround in the same direction */ + /* it is from the target value, for a quantized rgb */ + /* that is in that direction */ + for (k = 0; k < 3; k++) + base[k] = rgb[k]; + + for (ix = 0; ix < 8; ix++) { + double dist; + + /* Comp trial RGB */ + for (k = 0; k < 3; k++) { + tmp[k] = base[k]; + if (ix & (1 << k)) { + if (n & (1 << k)) /* Move in direction base point is in */ + tmp[k] += 1.0; + else + tmp[k] -= 1.0; + if (tmp[k] < 0.0 || tmp[k] > 255.0) + break; + } + } + if (k < 3) /* Trial is out of gamut */ + continue; + + /* Quantize it */ + ccast2YCbCr(NULL, ycc, tmp); + YCbCr2ccast(NULL, chval, ycc); +//printf("Trial RGB %f %f %f\n", tmp[0], tmp[1], tmp[2]); +//printf(" result %f %f %f\n", chval[0], chval[1], chval[2]); + + /* It's OK if it is eual or greater in the desired */ + /* direction than the input rgb. */ + /* Best would be closest to input. */ + /* Least worst would be what ? */ + + dist = 0.0; + for (k = 0; k < 3; k++) { + double tt; + if (n & (1 << k)) { + tt = chval[k] - base[k]; + } else { + tt = base[k] - chval[k]; + } + if (tt >= 0.0) + dist += 0.1 * tt; + else + dist += 2.0 * -tt; + } +//printf(" dist %f\n", dist); + + /* Pick it if it is at least in the right direction */ + if (dist < bdist) { + for (k = 0; k < 3; k++) { + rgb[k] = tmp[k]; + out[k] = chval[k]; + } + bdist = dist; + } + } + +#ifdef NEVER + if (bdist > 0.0 + && rgb[0] != 0.0 && rgb[0] != 1.0 && rgb[0] != 254.0 && rgb[0] != 255.0 + && rgb[1] != 0.0 && rgb[1] != 1.0 && rgb[1] != 254.0 && rgb[1] != 255.0 + && rgb[2] != 0.0 && rgb[2] != 1.0 && rgb[2] != 254.0 && rgb[2] != 255.0) { + printf("quant_rgb failed with bdist %f\n",bdist); + printf(" rgb in %f %f %f\n",base[0],base[1],base[2]); + printf(" returning rgb %f %f %f\n",rgb[0],rgb[1],rgb[2]); + printf(" resrgb %f %f %f\n",out[0],out[1],out[2]); + } +#endif /* NEVER */ +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - */ + +typedef struct optcntx { + int di; + double *val; /* Target */ + double (*ressur)[8][3]; /* resulting RGB surrounding values */ +} optcntx; + +static double optfunc(void *fdata, double tp[]) { + optcntx *cntx = (optcntx *)fdata; + int i, k; + double iw, tmp[3]; + double err = 0.0; + + /* Compute interpolated result */ + for (k = 0; k < 3; k++) + tmp[k] = 0.0; + + iw = 1.0; + for (i = 0; i < cntx->di; i++) { + if (tp[i] < 0.0) + err += 1000.0 * tp[i] * tp[i]; + else if (tp[i] > 1.0) + err += 1000.0 * (tp[i] - 1.0) * (tp[i] - 1.0); + + for (k = 0; k < 3; k++) + tmp[k] += tp[i] * (*cntx->ressur)[i][k]; + iw -= tp[i]; + } + for (k = 0; k < 3; k++) + tmp[k] += iw * (*cntx->ressur)[i][k]; + + /* Compute error */ + for (k = 0; k < 3; k++) { + double tt = tmp[k] - cntx->val[k]; + err += tt * tt; + } +//printf("Returning %f from %f %f\n",err,tp[0],tp[1]); + return err; +} + +/* Compute pattern */ +/* return the delta to the target */ +double get_ccast_dith(double ipat[DISIZE][DISIZE][3], double val[3]) { + double itpat[DISIZE][DISIZE][3], tpat[DISIZE][DISIZE][3]; + double irtpat[DISIZE][DISIZE][3], rtpat[DISIZE][DISIZE][3]; /* resulting for tpat */ + double berr = 0.0; + int n, k; + int x, y; + int i, j; + int ii; + struct { + int x, y; + } order[16] = { +// Dispersed: + { 3, 1 },{ 1, 3 },{ 1, 1 },{ 3, 0 },{ 3, 3 },{ 1, 2 },{ 3, 2 },{ 2, 0 }, + { 1, 0 },{ 0, 3 },{ 0, 1 },{ 2, 1 },{ 2, 2 },{ 0, 0 },{ 0, 2 },{ 2, 3 } +// Clustered: +// { 0, 0 },{ 0, 1 },{ 1, 1 },{ 1, 0 },{ 2, 0 },{ 3, 0 },{ 3, 1 },{ 2, 1 }, +// { 2, 2 },{ 3, 2 },{ 3, 3 },{ 2, 3 },{ 1, 3 },{ 1, 2 },{ 0, 2 },{ 0, 3 } + }; + int cix = 0; + int nsur = 8, ncomb = 4; + double sur[8][3]; /* RGB surrounding values to use */ + double ressur[8][3]; /* resulting RGB surrounding values */ + int bcc[8]; /* Best combination */ + double bw[8]; /* Best weight */ + int biw[8]; /* Best integer weight/16 */ + double dval[3]; /* Dithered value */ + double err[3], werr; + double errxs; + /* Tuning params */ + int nitters = 300; /* No itters */ + int rand_count = 150; + double rand_start = 4.5; /* Ramp with itters */ + double rand_end = 0.2; + double rand_pow = 1.5; + double mxerrxs = 2.5; /* accumulated error reset threshold */ + unsigned int randv = 0x1234; + + /* 32 bit pseudo random sequencer based on XOR feedback */ + /* generates number between 1 and 4294967295 */ +#define PSRAND32(S) (((S) & 0x80000000) ? (((S) << 1) ^ 0xa398655d) : ((S) << 1)) + + /* Locate the 8 surrounding RGB verticies */ + for (n = 0; n < 8; n++) { + for (k = 0; k < 3; k++) { + if (n & (1 << k)) + sur[n][k] = ceil(val[k]); + else + sur[n][k] = floor(val[k]); + } + +//printf("Input sur %d: %f %f %f\n",n,sur[n][0], sur[n][1], sur[n][2], sur[n][3]); + /* Figure out what RGB values to use to surround the point, */ + /* and what actual RGB values to expect from using them. */ + quant_rgb(n, ressur[n], sur[n]); +//printf("Quant sur %d: %f %f %f\n",n,sur[n][0], sur[n][1], sur[n][2], sur[n][3]); +//printf(" ressur %d: %f %f %f\n",n,ressur[n][0], ressur[n][1], ressur[n][2], ressur[n][3]); +// printf("\n"); + } + + /* Reduce this to unique surrounders */ + for (nsur = 0, i = 0; i < 8; i++) { + + /* Check if the i'th entry is already in the list */ + for (j = 0; j < nsur; j++) { + for (k = 0; k < 3; k++) { + if (ressur[i][k] != ressur[j][k]) + break; + } + if (k < 3) + continue; /* Unique */ + break; /* Duplicate */ + } + if (j < nsur) /* Duplicate */ + continue; + + /* Copy i'th to nsur */ + for (k = 0; k < 3; k++) { + sur[nsur][k] = sur[i][k]; + ressur[nsur][k] = ressur[i][k]; + } + nsur++; + } + +#ifdef DIAGNOSTICS + if (vv) { + printf("There are %d unique surrounders:\n",nsur); + for (n = 0; n < nsur; n++) { + printf("sur %f %f %f\n",ressur[n][0], ressur[n][1], ressur[n][2], ressur[n][3]); + } + } +#endif + + /* Use an optimzer to set the initial values using all the unique */ + /* surrounders. */ + { + double s[8]; + optcntx cntx; + + ncomb = nsur; + for (n = 0; n < nsur; n++) { + bcc[n] = n; + bw[n] = 1.0/nsur; + s[n] = 0.1; + } + + cntx.di = nsur-1; + cntx.val = val; + cntx.ressur = &ressur; + + powell(NULL, cntx.di, bw, s, 1e-4, 1000, optfunc, &cntx, NULL, NULL); + + /* Compute baricentric values */ + bw[nsur-1] = 0.0; + for (n = 0; n < (nsur-1); n++) { + if (bw[n] < 0.0) + bw[n] = 0.0; + else if (bw[n] > 1.0) + bw[n] = 1.0; + bw[nsur-1] += bw[n]; + } + if (bw[nsur-1] > 1.0) { /* They summed to over 1.0 */ + for (n = 0; n < (nsur-1); n++) + bw[n] *= 1.0/bw[nsur-1]; /* Scale them down */ + bw[nsur-1] = 1.0; + } + bw[nsur-1] = 1.0 - bw[nsur-1]; /* Remainder */ + } + + /* Check the result */ +#ifdef DIAGNOSTICS + if (vv) { + double tmp[3], err; + + /* Compute interpolated result */ + for (k = 0; k < 3; k++) + tmp[k] = 0.0; + for (n = 0; n < ncomb; n++) { + for (k = 0; k < 3; k++) + tmp[k] += bw[n] * ressur[bcc[n]][k]; + } + /* Compute error */ + err = 0.0; + for (k = 0; k < 3; k++) { + tmp[k] -= val[k]; + err += tmp[k] * tmp[k]; + } + err = sqrt(err); + for (n = 0; n < ncomb; n++) + printf("Comb %d weight %f rgb %f %f %f\n",bcc[n],bw[n],ressur[bcc[n]][0],ressur[bcc[n]][1],ressur[bcc[n]][2]); + printf("Error %f %f %f rms %f\n",tmp[0], tmp[1], tmp[2], err); + printf("\n"); + } +#endif + + /* Compute the number of pixels for each surounder value */ + { + int sw[8], rem; + + /* Sort the weightings from smallest to largest */ + for (n = 0; n < 8; n++) + sw[n] = n; + + for (i = 0; i < (ncomb-1); i++) { + for (j = i+1; j < ncomb; j++) { + if (bw[sw[j]] < bw[sw[i]]) { + int tt = sw[i]; /* Swap them */ + sw[i] = sw[j]; + sw[j] = tt; + } + } + } + + /* Compute the nearest integer weighting out of 16 */ + rem = 16; + for (i = 0; i < (ncomb-1) && rem > 0; i++) { + n = sw[i]; + biw[n] = (int)(16.0 * bw[n] + 0.5); + rem -= biw[n]; + if (rem <= 0) + rem = 0; + } + for (; i < ncomb; i++) { + n = sw[i]; + biw[n] = rem; + } +#ifdef DIAGNOSTICS + if (vv) { + for (n = 0; n < ncomb; n++) + printf("Comb %d iweight %i rgb %f %f %f\n",bcc[n],biw[n],ressur[bcc[n]][0],ressur[bcc[n]][1],ressur[bcc[n]][2]); + } +#endif + } + + /* Set the initial pattern according to the integer weighting */ + for (cix = 0, n = 0; n < ncomb; n++) { + for (i = 0; i < biw[n]; i++) { + x = order[cix].x; + y = order[cix].y; + cix++; + for (k = 0; k < 3; k++) { + tpat[x][y][k] = itpat[x][y][k] = sur[bcc[n]][k]; + rtpat[x][y][k] = irtpat[x][y][k] = ressur[bcc[n]][k]; + } + } + } + +#ifdef DIAGNOSTICS + /* Check initial pattern error */ + if (vv) { + printf("Input pat:\n"); + for (x = 0; x < DISIZE; x++) { + for (y = 0; y < DISIZE; y++) { + if (y > 0) + printf(", "); + printf("%3.0f %3.0f %3.0f",rtpat[x][y][0], rtpat[x][y][1], rtpat[x][y][2]); + } + printf("\n"); + } + } +#endif + + upsample(dval, rtpat); + + werr = 0.0; + for (k = 0; k < 3; k++) { + err[k] = dval[k] - val[k]; + if (fabs(err[k]) > werr) + werr = fabs(err[k]); + } + +#ifdef DIAGNOSTICS + if (vv) { + printf("Target %f %f %f -> %f %f %f\n", val[0], val[1], val[2], dval[0], dval[1], dval[2]); + printf("Error %f %f %f werr %f\n", err[0], err[1], err[2], werr); + } +#endif + + berr = werr; + for (x = 0; x < DISIZE; x++) { + for (y = 0; y < DISIZE; y++) { + for (k = 0; k < 3; k++) + ipat[x][y][k] = tpat[x][y][k]; + } + } + + /* Improve fit if needed */ + /* This is a bit stocastic */ + errxs = 0.0; + for (ii = 0; ii < nitters; ii++) { /* Until we give up */ + double corr[3]; /* Correction direction needed */ + double bdot; + int bx, by; + + double wycc, mm; + double cell[3]; /* Cell being modified value */ + double ccell[3]; /* Corrected cell */ + double pcell[3]; /* Proposed new cell value */ + + for (k = 0; k < 3; k++) + corr[k] = val[k] - dval[k]; +#ifdef DIAGNOSTICS + if (vv) + printf("corr needed %f %f %f\n", corr[0], corr[1], corr[2]); +#endif + + /* Scale it and limit it */ + for (k = 0; k < 3; k++) { + double dd = 16.0 * corr[k]; + if (dd >= 1.0) + dd = 1.0; + else if (dd <= -1.0) + dd = -1.0; + else + dd = 0.0; + corr[k] = dd; + } + + if (corr[0] == 0.0 && corr[1] == 0 && corr[2] == 0.0) { +#ifdef DIAGNOSTICS + if (vv) + printf("No correction possible - done\n"); +#endif + break; + } + +#ifdef DIAGNOSTICS + if (vv) + printf("scaled corr %f %f %f\n", corr[0], corr[1], corr[2]); +#endif + + /* Search dither cell and surrounder for a combination */ + /* that is closest to the change we want to make. */ + bdot = 1e6; + bx = by = n = 0; + for (x = 0; x < DISIZE; x++) { + double rlevel = rand_start + (rand_end - rand_start) + * pow((ii % rand_count)/rand_count, rand_pow); + + for (y = 0; y < DISIZE; y++) { + for (i = 0; i < ncomb; i++) { + double dot = 0.0; + for (k = 0; k < 3; k++) + dot += (ressur[bcc[i]][k] - rtpat[x][y][k]) * corr[k]; + + /* Ramp the randomness up */ +// dot += d_rand(0.0, 0.1 + (2.5-0.1) * ii/nitters); +// dot += d_rand(-rlevel, rlevel); + /* use a deterministic random element, so that */ + /* the dither patterns are repeatable. */ + randv = PSRAND32(randv); + dot += rlevel * 2.0 * ((randv - 1)/4294967294.0 - 0.5); + + if (dot <= 0.0) + dot = 1e7; + else { + dot = (dot - 1.0) * (dot - 1.0); + } +//printf("dot %f from sur %f %f %f to pat %f %f %f\n", dot, ressur[bcc[i]][0], ressur[bcc[i]][1], ressur[bcc[i]][2], rtpat[x][y][0], rtpat[x][y][1], rtpat[x][y][2]); + + if (dot < bdot) { + bdot = dot; + bx = x; + by = y; + n = i; + } + } + } + } + +#ifdef DIAGNOSTICS + if (vv) { + printf("Changing cell [%d][%d] %f %f %f with dot %f\n",bx,by,rtpat[bx][by][0],rtpat[bx][by][1],rtpat[bx][by][2],bdot); + printf(" to sur %d: %f %f %f\n",n, ressur[bcc[n]][0], ressur[bcc[n]][1], ressur[bcc[n]][2]); + } +#endif + + /* Substitute the best correction for this cell */ + for (k = 0; k < 3; k++) { + tpat[bx][by][k] = sur[bcc[n]][k]; + rtpat[bx][by][k] = ressur[bcc[n]][k]; + } + +#ifdef DIAGNOSTICS + if (vv) { + printf("Input pat:\n"); + for (x = 0; x < DISIZE; x++) { + for (y = 0; y < DISIZE; y++) { + if (y > 0) + printf(", "); + printf("%3.0f %3.0f %3.0f",rtpat[x][y][0], rtpat[x][y][1], rtpat[x][y][2]); + } + printf("\n"); + } + } +#endif + + upsample(dval, rtpat); + + werr = 0.0; + for (k = 0; k < 3; k++) { + err[k] = dval[k] - val[k]; + if (fabs(err[k]) > werr) + werr = fabs(err[k]); + } + + if (werr > berr) { + errxs += werr - berr; + } + + /* New best */ + if (werr < berr) { + berr = werr; + errxs = 0.0; + for (x = 0; x < DISIZE; x++) { + for (y = 0; y < DISIZE; y++) { + for (k = 0; k < 3; k++) { + ipat[x][y][k] = tpat[x][y][k]; + + itpat[x][y][k] = tpat[x][y][k]; + irtpat[x][y][k] = rtpat[x][y][k]; + } + } + } + } + +#ifdef DIAGNOSTICS + if (vv) { + printf("Target %f %f %f -> %f %f %f\n", val[0], val[1], val[2], dval[0], dval[1], dval[2]); + printf("Error %f %f %f werr %f\n", err[0], err[1], err[2], werr); + } +#endif + + if (berr < 0.11) { +#ifdef DIAGNOSTICS + if (vv) + printf("best error %f < 0.11 - give up\n",berr); +#endif + break; + } + + /* If we're not making progress, reset to the last best */ + if (errxs > mxerrxs) { +#ifdef DIAGNOSTICS + if (vv) + printf("Restarting at ii %d \n",ii); +#endif + errxs = 0.0; + for (x = 0; x < DISIZE; x++) { + for (y = 0; y < DISIZE; y++) { + for (k = 0; k < 3; k++) { + tpat[x][y][k] = itpat[x][y][k]; + rtpat[x][y][k] = irtpat[x][y][k]; + } + } + } + } + } + +#ifdef DIAGNOSTICS + if (vv) { + printf("Returning best error %f pat:\n",berr); + for (y = 0; y < DISIZE; y++) { + for (x = 0; x < DISIZE; x++) { + if (x > 0) + printf(", "); + printf("%3.0f %3.0f %3.0f",ipat[x][y][0], ipat[x][y][1], ipat[x][y][2]); + } + printf("\n"); + } + } +#endif + + return berr; +} + +/* ====================================================================================== */ + +#ifdef STANDALONE_TEST + +int +main(int argc, + char *argv[] +) { + double val[3], out[3], err; + double ipat[DISIZE][DISIZE][3]; + double aerr, acount, xerr; + int x, y; + int i, j; + int k, nn; + + printf("Hi there\n"); + + rand32(time(NULL)); + +#ifdef TEST_POINT + for (x = 0; x < DISIZE; x++) { + for (y = 0; y < DISIZE; y++) { + ipat[x][y][0] = 0.0; + ipat[x][y][1] = 0.0; + ipat[x][y][2] = 0.0; + } + } + + ipat[5][4][0] = 255.0; + ipat[5][4][1] = 255.0; + ipat[5][4][2] = 255.0; + + upsample(NULL, ipat); + +#else + +#ifdef NEVER + +// val[0] = 201.500000; +// val[1] = 115.533403; +// val[2] = 76.300000; + +// val[0] = 255.000000; +// val[1] = 115.533403; +// val[2] = 255.000000; + +// val[0] = 221.689875; +// val[1] = 29.593255; +// val[2] = 140.820878; + +// val[0] = 212.377797; +// val[1] = 228.338903; +// val[2] = 70.153296; + +// val[0] = 231.554511; +// val[1] = 0.000000; +// val[2] = 51.958048; + +// val[0] = 255.000000; +// val[1] = 144.768052; +// val[2] = 179.737212; + +// val[0] = 194.854956; +// val[1] = 41.901887; +// val[2] = 20.434793; + +// val[0] = 250.100121; +// val[1] = 83.484217; +// val[2] = 42.867603; + +// val[0] = 255.000000; +// val[1] = 255.000000; +// val[2] = 228.534759; + + // 1.71 -> 1.58, rand 2.2 +// val[0] = 255.000000; +// val[1] = 176.894769; +// val[2] = 8.932806; + + // 1.54 -> 0.762592, rand 0.3 +// val[0] = 216.873703; +// val[1] = 250.908094; + + // 1.05 +// val[0] = 167.284458; +// val[1] = 248.945210; +// val[2] = 199.023452; + + // 1.07 +// val[0] = 211.045184; +// val[1] = 27.825141; +// val[2] = 63.883148; + + // 1.928 +// val[0] = 255.000000; +// val[1] = 0.439284; +// val[2] = 210.928135; + + // 1.278 +// val[0] = 218.693614; +// val[1] = 222.890101; +// val[2] = 174.779727; + + // 1.334501 +// val[0] = 253.931573; +// val[1] = 230.278945; +// val[2] = 185.677389; + +// printf("In RGB %f %f %f\n",val[0],val[1],val[2]); +// ccast2YCbCr(NULL, out, val); +// printf("YCbCr %f %f %f\n",out[0],out[1],out[2]); +// YCbCr2ccast(NULL, out, out); +// printf("RGB %f %f %f\n",out[0],out[1],out[2]); + + + vv = 1; + err = get_ccast_dith(ipat, val); + + printf("Got pat with err %f\n",err); + +#else + aerr = 0.0; + acount = 0.0; + xerr = 0.0; + for (nn = 0; nn < 200000; nn++) { + + for (k = 0; k < 3; k++) { + val[k] = d_rand(-5.0, 255.0 + 5.0); + if (val[k] < 0.0) + val[k] = 0.0; + else if (val[k] > 255.0) + val[k] = 255.0; + } + + err = get_ccast_dith(ipat, val); + + if (err >= 1.5 || + ( err >= 1.0 + && val[0] != 0.0 && val[0] != 255.0 + && val[1] != 0.0 && val[1] != 255.0 + && val[2] != 0.0 && val[2] != 255.0)) { + printf("Target RGB %f %f %f, err %f\n", val[0], val[1], val[2],err); +// vv = 1; +// comput_pat(ipat, val); +// break; + } + + aerr += err; + acount++; + if (err > xerr) + xerr = err; + } + + aerr /= acount; + printf("After %d trials, aerr = %f, maxerr = %f\n",nn,aerr,xerr); +#endif +#endif + + return 0; +} + +#endif /* STANDALONE_TEST */ -- cgit v1.2.3