summaryrefslogtreecommitdiff
path: root/spectro/average.c
diff options
context:
space:
mode:
Diffstat (limited to 'spectro/average.c')
-rw-r--r--spectro/average.c438
1 files changed, 321 insertions, 117 deletions
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]);
+}