From 3db384424bd7398ffbb7a355cab8f15f3add009f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B6rg=20Frings-F=C3=BCrst?= Date: Sun, 2 Oct 2016 19:24:58 +0200 Subject: New upstream version 1.9.1+repack --- spectro/average.c | 438 +++++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 321 insertions(+), 117 deletions(-) (limited to 'spectro/average.c') diff --git a/spectro/average.c b/spectro/average.c index 5a8fb3f..e961262 100644 --- a/spectro/average.c +++ b/spectro/average.c @@ -23,6 +23,7 @@ */ + #undef DEBUG #define verbo stdout @@ -34,10 +35,15 @@ #include "copyright.h" #include "aconfig.h" #include "numlib.h" +#include "sort.h" #include "cgats.h" #include "xicc.h" #include "insttypes.h" +static double average(double *vals, int nvals); +static double median(double *vals, int nvals); +static void geommed(double res[3], double vals[][3], int nvals); + void usage(char *diag, ...) { int i; fprintf(stderr,"Average or merge values in .ti3 like files, Version %s\n",ARGYLL_VERSION_STR); @@ -52,6 +58,10 @@ void usage(char *diag, ...) { fprintf(stderr,"Author: Graeme W. Gill, licensed under the AGPL Version 3\n"); fprintf(stderr,"usage: average [-options] input1.ti3 input2.ti3 ... output.ti3\n"); fprintf(stderr," -v Verbose\n"); + fprintf(stderr," -e Median rather than average\n"); + fprintf(stderr," -g Geometric Median of PCS in encoded space\n"); + fprintf(stderr," -L Geometric Median of PCS in L*a*b* space\n"); + fprintf(stderr," -X Geometric Median of PCS in XYZ space\n"); fprintf(stderr," -m Merge rather than average\n"); fprintf(stderr," input1.ti3 First input file\n"); fprintf(stderr," input2.ti3 Second input file\n"); @@ -69,6 +79,8 @@ struct _inpinfo { int main(int argc, char *argv[]) { int fa,nfa; /* current argument we're looking at */ int verb = 0; + int domedian = 0; /* Median rather than average */ + int dogeom = 0; /* Do geometric median of PCS, 2 = Lab, 3 = PCS */ int domerge = 0; /* Merge rather than average */ int ninps = 0; /* Number of input files */ @@ -78,10 +90,12 @@ int main(int argc, char *argv[]) { cgats_set_elem *setel; /* Array of set value elements */ int *flags; /* Point to destination of set */ - int nchan; /* Number of device channels */ + int nchan = 0; /* Number of device channels */ int chix[ICX_MXINKS]; /* Device chanel indexes */ - int pcsix[3]; /* PCS chanel indexes */ - int isLab = 0; + int npcs = 0; + + int haspcs[2] = { 0 }; /* Has Lab, XYZ */ + int pcsix[3][3]; /* Lab, XYZ chanel indexes */ int i, j, n; @@ -110,13 +124,33 @@ int main(int argc, char *argv[]) { if (argv[fa][1] == '?') usage("Usage requested"); + /* Median */ + else if (argv[fa][1] == 'e') { + domedian = 1; + } + + /* Geometric Median of PCS */ + else if (argv[fa][1] == 'g') { + dogeom = 1; + } + + /* Geometric Median of PCS in L*a*b* */ + else if (argv[fa][1] == 'L') { + dogeom = 2; + } + + /* Geometric Median of PCS in XYZ */ + else if (argv[fa][1] == 'X') { + dogeom = 3; + } + /* Merge */ - else if (argv[fa][1] == 'm' || argv[fa][1] == 'M') { + else if (argv[fa][1] == 'm') { domerge = 1; } /* Verbosity */ - else if (argv[fa][1] == 'v' || argv[fa][1] == 'V') { + else if (argv[fa][1] == 'v') { verb = 1; } @@ -192,8 +226,9 @@ int main(int argc, char *argv[]) { ocg->add_field(ocg, n, inps[0].c->t[n].fsym[i], inps[0].c->t[n].ftype[i]); } + /* If more than one file, must be merging or averaging between files */ if (ninps > 1) { - /* Duplicate all of the data */ + /* Duplicate all of the data or first file to output file */ if ((setel = (cgats_set_elem *)malloc( sizeof(cgats_set_elem) * inps[0].c->t[n].nfields)) == NULL) error("Malloc failed!"); @@ -207,11 +242,11 @@ int main(int argc, char *argv[]) { } /* Figure out the indexes of the device channels */ - { + if (inps[0].c->find_kword(inps[0].c, 0, "COLOR_REP") < 0) { + warning("Input file '%s' doesn't contain keyword COLOR_REP", inps[0].name); + } else { int ti; char *buf; - char *xyzfname[3] = { "XYZ_X", "XYZ_Y", "XYZ_Z" }; - char *labfname[3] = { "LAB_L", "LAB_A", "LAB_B" }; char *outc; int nmask; char *bident; @@ -227,13 +262,6 @@ int main(int argc, char *argv[]) { error("COLOR_REP '%s' invalid", inps[0].c->t[0].kdata[ti]); *outc++ = '\000'; - if (strcmp(outc, "XYZ") == 0) { - isLab = 0; - } else if (strcmp(outc, "LAB") == 0) { - isLab = 1; - } else - error("COLOR_REP '%s' invalid (Neither XYZ nor LAB)", inps[0].c->t[0].kdata[ti]); - if ((nmask = icx_char2inkmask(buf)) == 0) { error ("File '%s' keyword COLOR_REP has unknown device value '%s'",inps[0].name,buf); } @@ -256,19 +284,32 @@ int main(int argc, char *argv[]) { error ("Field %s is wrong type",fname); chix[j] = ii; } - - /* Find PCS fields */ - for (j = 0; j < 3; j++) { - int ii; - - if ((ii = inps[0].c->find_field(inps[0].c, 0, isLab ? labfname[j] : xyzfname[j])) < 0) - error ("Input file doesn't contain field %s",isLab ? labfname[j] : xyzfname[j]); - if (inps[0].c->t[0].ftype[ii] != r_t) - error ("Field %s is wrong type",isLab ? labfname[j] : xyzfname[j]); - pcsix[j] = ii; + } + + /* Figure out the indexes of the PCS channels, if any */ + { + int npcs; + char *fname[2][3] = { { "LAB_L", "LAB_A", "LAB_B" }, + { "XYZ_X", "XYZ_Y", "XYZ_Z" } }; + + /* For Lab and XYZ */ + for (j = 0; j < 2; j++) { + for (npcs = 0; npcs < 3; npcs++) { + int ii; + + if ((ii = inps[0].c->find_field(inps[0].c, 0, fname[j][npcs])) < 0) + break; /* Try next or give up */ + + if (inps[0].c->t[0].ftype[ii] != r_t) + error ("Field %s is wrong type",fname[j][npcs]); + pcsix[j][npcs] = ii; + } + if (npcs == 3) + haspcs[j] = 1; } - free(bident); } + if (!haspcs[0] && !haspcs[1]) + warning("No PCS fields found - hope that's OK!"); if (!domerge && verb) { printf("Averaging the following fields:"); @@ -294,7 +335,8 @@ int main(int argc, char *argv[]) { printf("\n"); } - /* Get ready to add more values to output */ + /* Get ready to add more values to output, */ + /* for merging or averaging within one file. */ if ((setel = (cgats_set_elem *)malloc( sizeof(cgats_set_elem) * inps[0].c->t[0].nfields)) == NULL) error("Malloc failed!"); @@ -302,13 +344,24 @@ int main(int argc, char *argv[]) { /* If averaging values within the one file */ if (ninps == 1) { int *valdone; - double npatches; - int k; + int npat; + double *vlist; + double (*v3list)[3] = NULL; + int k, e; n = 0; /* Output set index */ if ((valdone = (int *)calloc(inps[0].c->t[0].nsets, sizeof(int))) == NULL) error("Malloc failed!"); + if ((vlist = (double *)calloc(inps[0].c->t[0].nsets, sizeof(double))) == NULL) + error("Malloc failed!"); + + if (dogeom && (haspcs[0] || haspcs[1])) { + if ((v3list = (double (*)[3])calloc(inps[0].c->t[0].nsets, 3 * sizeof(double))) == NULL) + error("Malloc failed!"); + } + + /* For each patch */ for (i = 0; i < inps[0].c->t[0].nsets; i++) { if (valdone[i]) @@ -316,49 +369,8 @@ int main(int argc, char *argv[]) { inps[0].c->get_setarr(inps[0].c, 0, i, setel); ocg->add_setarr(ocg, 0, setel); - npatches = 1.0; - - /* Locate and patches with matching device values */ - for (k = i+1; k < inps[0].c->t[0].nsets; k++) { - - /* Check if the device values match */ - for (j = 0; j < nchan; j++) { - double diff; - - diff = *((double *)inps[0].c->t[0].fdata[i][chix[j]]) - - *((double *)inps[0].c->t[0].fdata[k][chix[j]]); - - if (fabs(diff) > 0.001) { - break; - } - } - if (j < nchan) { - continue; - } - /* Add all the non-device real field values */ - for (j = 0; j < inps[0].c->t[0].nfields; j++) { - int jj; - - /* Only real types */ - if (inps[0].c->t[0].ftype[j] != r_t) - continue; - - /* Not device channels */ - for (jj = 0; jj < nchan; jj++) { - if (chix[jj] == j) - break; - } - if (jj < nchan) - continue; - - *((double *)ocg->t[0].fdata[n][j]) - += *((double *)inps[0].c->t[0].fdata[k][j]); - } - npatches++; - valdone[k] = 1; - } - /* Average them out */ + /* For each non-device real field values */ for (j = 0; j < inps[0].c->t[0].nfields; j++) { int jj; @@ -374,20 +386,102 @@ int main(int argc, char *argv[]) { if (jj < nchan) continue; - *((double *)ocg->t[0].fdata[n][j]) /= npatches; + /* Locate any patches (including starting patch) with matching device values */ + npat = 0; + for (k = i; k < inps[0].c->t[0].nsets; k++) { + + /* Check if the device values match */ + for (e = 0; e < nchan; e++) { + double diff; + + diff = *((double *)inps[0].c->t[0].fdata[i][chix[e]]) + - *((double *)inps[0].c->t[0].fdata[k][chix[e]]); + + if (fabs(diff) > 0.001) { + break; + } + } + if (e < nchan) { + continue; + } + + vlist[npat++] = *((double *)inps[0].c->t[0].fdata[k][j]); + valdone[k] = 1; + } + if (domedian) + *((double *)ocg->t[0].fdata[n][j]) = median(vlist, npat); + else + *((double *)ocg->t[0].fdata[n][j]) = average(vlist, npat); + } + + /* Override per-component average/median if PCS Geometric Median */ + if (dogeom && (haspcs[0] || haspcs[1])) { + double res[3]; + + /* For Lab and XYZ */ + for (j = 0; j < 2; j++) { + + if (haspcs[j] == 0) + continue; + + /* Locate any patches (including starting patch) with matching device values */ + npat = 0; + for (k = i; k < inps[0].c->t[0].nsets; k++) { + + /* Check if the device values match */ + for (e = 0; e < nchan; e++) { + double diff; + + diff = *((double *)inps[0].c->t[0].fdata[i][chix[e]]) + - *((double *)inps[0].c->t[0].fdata[k][chix[e]]); + + if (fabs(diff) > 0.001) { + break; + } + } + if (e < nchan) { + continue; + } + + v3list[npat][0] = *((double *)inps[0].c->t[0].fdata[k][pcsix[j][0]]); + v3list[npat][1] = *((double *)inps[0].c->t[0].fdata[k][pcsix[j][1]]); + v3list[npat][2] = *((double *)inps[0].c->t[0].fdata[k][pcsix[j][2]]); + + if (j == 0 && dogeom == 3) /* Lab and want XYZ */ + icmLab2XYZ(&icmD50_100, v3list[npat], v3list[npat]); + else if (j == 1 && dogeom == 2) /* XYZ and want Lab */ + icmXYZ2Lab(&icmD50_100, v3list[npat], v3list[npat]); + + npat++; + } + geommed(res, v3list, npat); + + if (j == 0 && dogeom == 3) + icmXYZ2Lab(&icmD50_100, res, res); + else if (j == 1 && dogeom == 2) + icmLab2XYZ(&icmD50_100, res, res); + + *((double *)ocg->t[0].fdata[n][pcsix[j][0]]) = res[0]; + *((double *)ocg->t[0].fdata[n][pcsix[j][1]]) = res[1]; + *((double *)ocg->t[0].fdata[n][pcsix[j][2]]) = res[2]; + } } n++; /* One more output set */ } + if (v3list != NULL) + free(v3list); + free(vlist); free(valdone); /* Averaging patches between identical files, */ /* or concatenating (merging) several files */ } else { - /* Process all the other input files */ + + /* Check/process all the other input files */ for (n = 1; n < ninps; n++) { - /* Check all the fields match */ + /* Check all the fields match the first file */ if (inps[0].c->t[0].nfields != inps[n].c->t[0].nfields) error ("File '%s' has %d fields, file '%s has %d", inps[n].name, inps[n].c->t[0].nfields, inps[0].name, inps[0].c->t[0].nfields); @@ -405,15 +499,15 @@ int main(int argc, char *argv[]) { } } else { /* Averaging */ - /* Check the number of values matches */ + + /* Check the number of patches matches the first file */ if (inps[0].c->t[0].nsets != inps[n].c->t[0].nsets) error ("File '%s' has %d sets, file '%s has %d", inps[n].name, inps[n].c->t[0].nsets, inps[0].name, inps[0].c->t[0].nsets); - - /* Add the numeric field values to corresponding output */ + /* For all the patches: */ for (i = 0; i < inps[n].c->t[0].nsets; i++) { - /* Check that the device values match */ + /* Check that the device values match the first file */ for (j = 0; j < nchan; j++) { double diff; diff = *((double *)inps[0].c->t[0].fdata[i][chix[j]]) @@ -423,53 +517,100 @@ int main(int argc, char *argv[]) { error ("File '%s' set %d has field '%s' value that differs from '%s'", inps[n].name, i+1, inps[n].c->t[0].fsym[j], inps[0].name); } - - /* Add all the non-device real field values */ - for (j = 0; j < inps[0].c->t[0].nfields; j++) { - int jj; - - /* Only real types */ - if (inps[0].c->t[0].ftype[j] != r_t) - continue; - - /* Not device channels */ - for (jj = 0; jj < nchan; jj++) { - if (chix[jj] == j) - break; - } - if (jj < nchan) - continue; - - *((double *)ocg->t[0].fdata[i][j]) - += *((double *)inps[n].c->t[0].fdata[i][j]); - } } } } - - /* If averaging, divide out the number of files */ + + /* If averaging */ if (!domerge) { + int npat; + double *vlist; + double (*v3list)[3] = NULL; + + if ((vlist = (double *)calloc(ninps, sizeof(double))) == NULL) + error("Malloc failed!"); + + if (dogeom && (haspcs[0] || haspcs[1])) { + if ((v3list = (double (*)[3])calloc(inps[0].c->t[0].nsets, 3 * sizeof(double))) == NULL) + error("Malloc failed!"); + } + + /* For all the non-device real field values */ + for (j = 0; j < inps[0].c->t[0].nfields; j++) { + int jj; + + /* Only real types */ + if (inps[0].c->t[0].ftype[j] != r_t) + continue; + + /* Not device channels */ + for (jj = 0; jj < nchan; jj++) { + if (chix[jj] == j) + break; + } + if (jj < nchan) + continue; - for (i = 0; i < inps[n].c->t[0].nsets; i++) { - - for (j = 0; j < inps[0].c->t[0].nfields; j++) { - int jj; - - /* Only real types */ - if (inps[0].c->t[0].ftype[j] != r_t) - continue; - - /* Not device channels */ - for (jj = 0; jj < nchan; jj++) { - if (chix[jj] == j) - break; + /* For each patch */ + for (i = 0; i < inps[n].c->t[0].nsets; i++) { + + /* For all input files */ + npat = 0; + for (n = 0; n < ninps; n++) { + vlist[npat++] = *((double *)inps[n].c->t[0].fdata[i][j]); } - if (jj < nchan) - continue; + + if (domedian) + *((double *)ocg->t[0].fdata[i][j]) = median(vlist, npat); + else + *((double *)ocg->t[0].fdata[i][j]) = average(vlist, npat); + } + } + + /* Override per-component average/median if PCS Geometric Median */ + if (dogeom && (haspcs[0] || haspcs[1])) { + double res[3]; + + /* For each patch */ + for (i = 0; i < inps[n].c->t[0].nsets; i++) { + + /* For Lab and XYZ */ + for (j = 0; j < 2; j++) { + + if (haspcs[j] == 0) + continue; - *((double *)ocg->t[0].fdata[i][j]) /= (double)ninps; + /* For all input files */ + npat = 0; + for (n = 0; n < ninps; n++) { + v3list[npat][0] = *((double *)inps[n].c->t[0].fdata[i][pcsix[j][0]]); + v3list[npat][1] = *((double *)inps[n].c->t[0].fdata[i][pcsix[j][1]]); + v3list[npat][2] = *((double *)inps[n].c->t[0].fdata[i][pcsix[j][2]]); + + if (j == 0 && dogeom == 3) /* Lab and want XYZ */ + icmLab2XYZ(&icmD50_100, v3list[npat], v3list[npat]); + else if (j == 1 && dogeom == 2) /* XYZ and want Lab */ + icmXYZ2Lab(&icmD50_100, v3list[npat], v3list[npat]); + + npat++; + } + geommed(res, v3list, npat); + + if (j == 0 && dogeom == 3) + icmXYZ2Lab(&icmD50_100, res, res); + else if (j == 1 && dogeom == 2) + icmLab2XYZ(&icmD50_100, res, res); + + *((double *)ocg->t[0].fdata[i][pcsix[j][0]]) = res[0]; + *((double *)ocg->t[0].fdata[i][pcsix[j][1]]) = res[1]; + *((double *)ocg->t[0].fdata[i][pcsix[j][2]]) = res[2]; + } } } + + if (v3list != NULL) + free(v3list); + free(vlist); } } @@ -490,6 +631,69 @@ int main(int argc, char *argv[]) { } +static double average(double *vals, int nvals) { + double rv; + int i; + + for (rv = 0.0, i = 0; i < nvals; i++) + rv += vals[i]; + + if (nvals > 0) + rv /= (double)nvals; + + return rv; +} + +static double median(double *vals, int nvals) { + if (nvals < 3) + return average(vals, nvals); + +#define HEAP_COMPARE(A,B) (A < B) + HEAPSORT(double,vals,nvals); + + if ((nvals & 1) != 0) + return vals[nvals/2]; + else + return 0.5 * (vals[nvals/2] + vals[nvals/2-1]); +} +/* Compute Geometric Median of PCS values */ +/* using Weiszfeld's algorithm. */ +static void geommed(double res[3], double vals[][3], int nvals) { + int i, j; + + /* Start with mean value */ + icmSet3(res, 0.0); + for (i = 0; i < nvals; i++) + icmAdd3(res, res, vals[i]); + icmScale3(res, res, 1.0/(double)nvals); + +//printf("\n~1 average = %f %f %f\n", res[0], res[1], res[2]); + + /* Itterate to approach Geometric Mean */ + for (j = 0; j < 20; j++) { + double tv[3], tt; + int k; + + icmSet3(tv, 0.0); + tt = 0.0; + for (k = i = 0; i < nvals; i++) { + double norm = icmNorm33(vals[i], res); + if (norm < 1e-6) + continue; + tv[0] += vals[i][0]/norm; + tv[1] += vals[i][1]/norm; + tv[2] += vals[i][2]/norm; + tt += 1.0/norm; + k++; +//printf("Norm = %f, tv = %f %f %f, tt = %f\n",norm, tv[0], tv[1], tv[2], tt); + } + if (k > 0) + icmScale3(res, tv, 1.0/tt); +//printf("~1 res = %f %f %f\n", res[0], res[1], res[2]); + } + +//printf("~1 geomm = %f %f %f\n", res[0], res[1], res[2]); +} -- cgit v1.2.3