gronod/argyllcms
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 11 Wiki Activity
Files
fix/win32-pipe-peek-named-pipe
argyllcms/spectro/average.c
T
gronod d0209ec246
Linux x86_64 Release Build / build (push) Failing after 1m5s
Details
Windows x86_64 Release Build / build (push) Has been cancelled
Details
macOS Release Build / build (arm64, macOS_arm64_bin, arm64-apple-darwin, macos-latest) (push) Has been cancelled
Details
macOS Release Build / build (x86_64, osx10.6_x86_64_bin, x86_64-apple-darwin, macos-13) (push) Has been cancelled
Details
Fix macOS UI link errors by adding missing ui.h includes
2026-08-22 21:27:56 +01:00

920 lines
24 KiB
C
Raw Permalink Blame History

/*
* Argyll Color Management System
* Average one or more .ti3 (or other CGATS like) file values together.
* If just one file is supplied, all patches within it with
* the same device value are averaged together.
*
* Author: Graeme W. Gill
* Date: 18/1/2011
*
* Copyright 2005, 2010, 2011, 2025 Graeme W. Gill
* All rights reserved.
*
* This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
* see the License.txt file for licencing details.
*
* (based on splitti3.c)
*/
/*
* TTBD:
Should probably re-index SAMPLE_ID field
Geometric mean is limited to just 3 channels.
*/
#undef DEBUG
#define verbo stdout
#include <stdio.h>
#include <string.h>
#include <sys/types.h>
#include <time.h>
#include "copyright.h"
#include "aconfig.h"
#include "numlib.h"
#include "sort.h"
#include "cgats.h"
#include "conv.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);
if (diag != NULL) {
va_list args;
fprintf(stderr," Diagnostic: ");
va_start(args, diag);
vfprintf(stderr, diag, args);
va_end(args);
fprintf(stderr,"\n");
}
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");
fprintf(stderr," ... etc.\n");
fprintf(stderr," output.ti3 Resulting averaged or merged output file\n");
exit(1);
}
/* Information about each file */
struct _inpinfo {
char name[MAXNAMEL+1];
cgats *c;
}; typedef struct _inpinfo inpinfo;
#include "ui.h"
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 or Dev, 2 = Lab, 3 = PCS */
int domerge = 0; /* Merge rather than average */
int ninps = 0; /* Number of input files */
inpinfo *inps; /* Input file info. inp[ninp] == output file */
cgats *ocg; /* Copy of output cgats * */
cgats_set_elem *setel; /* Array of set value elements */
int *flags; /* Point to destination of set */
int inp = -1; /* 0 for output, 1 for input device files */
int nchan = 0; /* Number of device channels */
int chix[ICX_MXINKS]; /* Device chanel indexes */
int haspcs[2] = { 0 }; /* Has [Lab, XYZ] */
int pcsix[2][3]; /* Lab, XYZ chanel indexes */
int dti;
int i, j, n;
error_program = "average";
if (argc < 2)
usage("Too few arguments (%d, minimum is 2)",argc-1);
/* Process the arguments */
for (fa = 1;fa < argc;fa++) {
nfa = fa; /* skip to nfa if next argument is used */
if (argv[fa][0] == '-') { /* Look for any flags */
char *na = NULL; /* next argument after flag, null if none */
if (argv[fa][2] != '\000')
na = &argv[fa][2]; /* next is directly after flag */
else {
if ((fa+1) < argc) {
if (argv[fa+1][0] != '-') {
nfa = fa + 1;
na = argv[nfa]; /* next is seperate non-flag argument */
}
}
}
if (argv[fa][1] == '?')
usage("Usage requested");
/* Median */
else if (argv[fa][1] == 'e') {
domedian = 1;
}
/* Geometric Median of PCS or Dev */
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') {
domerge = 1;
}
/* Verbosity */
else if (argv[fa][1] == 'v') {
verb = 1;
}
else {
usage("Unknown flag '%c'",argv[fa][1]);
}
} else if (argv[fa][0] != '\000') {
/* Get the next filename */
if (ninps == 0)
inps = (inpinfo *)malloc(sizeof(inpinfo));
else
inps = (inpinfo *)realloc(inps, (ninps+1) * sizeof(inpinfo));
if (inps == NULL)
error("Malloc failed in allocating space for file info.");
memset((void *)&inps[ninps], 0, sizeof(inpinfo));
strncpy(inps[ninps].name,argv[fa],MAXNAMEL);
inps[ninps].name[MAXNAMEL] = '\000';
ninps++;
} else {
break;
}
}
if (ninps < 2)
error("Must be at least one input and one output file specified");
ninps--; /* Number of inputs */
/* Open and read each input file, and create output file */
for (n = 0; n <= ninps; n++) {
if ((inps[n].c = new_cgats()) == NULL)
error("Failed to create cgats object for file '%s'",inps[n].name);
if (n < ninps) { /* If input file, read it */
inps[n].c->add_other(inps[n].c, ""); /* Allow any signature file */
if (inps[n].c->read_name(inps[n].c, inps[n].name))
error("CGATS file '%s' read error : %s",inps[n].name,inps[n].c->e.m);
if (inps[n].c->ntables < 1)
error ("Input file '%s' doesn't contain at least one table",inps[n].name);
}
}
ocg = inps[ninps].c; /* Alias for output file */
/* Check the device/profile type tag to see if input type */
for (n = 0; n < inps[0].c->ntables; n++) {
if ((dti = inps[n].c->find_kword(inps[n].c, 0, "DEVICE_CLASS")) < 0) {
if (n == 0)
inp = 0; /* Default to output */
continue;
}
if (strcmp(inps[n].c->t[0].kdata[dti],"INPUT") == 0
|| strcmp(inps[n].c->t[0].kdata[dti],"EMISINPUT") == 0) {
if (n == 0)
inp = 1;
else if (inp != 1)
error ("Input file '%s' has conflicting DEVICE_CLASS",inps[n].name);
} else {
if (n == 0)
inp = 0;
else if (inp != 0)
error ("Input file '%s' has conflicting DEVICE_CLASS",inps[n].name);
}
}
/* Duplicate everything from the first input file into the output file. */
for (n = 0; n < inps[0].c->ntables; n++) {
if (inps[0].c->t[n].tt == cgats_X) {
ocg->add_other(ocg, inps[0].c->cgats_type);
ocg->add_table(ocg, tt_other, 0);
} else if (inps[0].c->t[n].tt == tt_other) {
int oi;
oi = ocg->add_other(ocg, inps[0].c->others[inps[0].c->t[n].oi]);
ocg->add_table(ocg, tt_other, oi);
} else {
ocg->add_table(ocg, inps[0].c->t[n].tt, 0);
}
/* Duplicate all the keywords */
for (i = 0; i < inps[0].c->t[n].nkwords; i++) {
//printf("~1 table %d, adding keyword '%s'\n",n,inps[0].c->t[n].ksym[i]);
ocg->add_kword(ocg, n, inps[0].c->t[n].ksym[i], inps[0].c->t[n].kdata[i], NULL);
}
/* Duplicate all of the fields */
for (i = 0; i < inps[0].c->t[n].nfields; i++) {
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 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!");
for (i = 0; i < inps[0].c->t[n].nsets; i++) {
inps[0].c->get_setarr(inps[0].c, n, i, setel);
ocg->add_setarr(ocg, n, setel);
}
free(setel);
}
}
/* 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 *outc;
int nmask;
char *bident;
if ((ti = inps[0].c->find_kword(inps[0].c, 0, "COLOR_REP")) < 0)
error("Input file '%s' doesn't contain keyword COLOR_REP", inps[0].name);
if ((buf = strdup(inps[0].c->t[0].kdata[ti])) == NULL)
error("Malloc failed");
/* Split COLOR_REP into device and PCS space */
if ((outc = strchr(buf, '_')) == NULL)
error("COLOR_REP '%s' invalid", inps[0].c->t[0].kdata[ti]);
*outc++ = '\000';
if ((nmask = icx_char2inkmask(buf)) == 0) {
/* Hmm. This might be an input reference. */
if (strcmp(buf, "XYZ") == 0 || strcmp(buf, "LAB") == 0) {
if ((nmask = icx_char2inkmask(outc)) == 0) {
error ("File '%s' keyword COLOR_REP has unknown device value '%s'",inps[0].name,outc);
}
} else {
error ("File '%s' keyword COLOR_REP has unknown device value '%s'",inps[0].name,buf);
}
}
nchan = icx_noofinks(nmask);
bident = icx_inkmask2char(nmask, 0);
/* Find device fields */
for (j = 0; j < nchan; j++) {
int ii, imask;
char fname[100];
imask = icx_index2ink(nmask, j);
sprintf(fname,"%s_%s",nmask == ICX_W || nmask == ICX_K ? "GRAY" : bident,
icx_ink2char(imask));
if ((ii = inps[0].c->find_field(inps[0].c, 0, fname)) < 0)
error ("Input file doesn't contain field %s",fname);
if (inps[0].c->t[0].ftype[ii] != r_t)
error ("Field %s is wrong type",fname);
chix[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;
}
}
if (!haspcs[0] && !haspcs[1]) {
if (inp)
error("Need a PCS to average for input device chart");
else
warning("No PCS fields found - hope that's OK!");
}
if (inp && dogeom != 0) {
if (nchan != 3)
error ("Can't do geometric mean for %d device channels",nchan);
if (dogeom > 1) {
warning("Geometric mean of XYZ or Lab ignored for input type since device values are being averaged");
dogeom = 1;
}
}
if (!domerge && verb) {
printf("Averaging the following fields:");
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;
}
if (inp) {
/* Not PCS channels */
int k;
for (k = 0; k < 2; k++) {
if (haspcs[k] == 0)
continue;
for (jj = 0; jj < 3; jj++) {
if (pcsix[k][jj] == j)
break;
}
if (jj < 3)
break;
}
if (k < 2)
continue;
} else {
/* Not device channels */
for (jj = 0; jj < nchan; jj++) {
if (chix[jj] == j)
break;
}
if (jj < nchan)
continue;
}
printf(" %s",inps[0].c->t[0].fsym[j]);
}
printf("\n");
}
/* 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!");
/* If averaging values within the one file */
if (ninps == 1) {
int *valdone;
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])
continue;
inps[0].c->get_setarr(inps[0].c, 0, i, setel);
ocg->add_setarr(ocg, 0, setel);
/* For each 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;
if (inp) {
/* Not PCS channels */
int k;
for (k = 0; k < 2; k++) {
if (haspcs[k] == 0)
continue;
for (jj = 0; jj < 3; jj++) {
if (pcsix[k][jj] == j)
break;
}
if (jj < 3)
break;
}
if (k < 2)
continue;
/* Locate any patches (including starting patch) with matching PCS values */
npat = 0;
for (k = i; k < inps[0].c->t[0].nsets; k++) {
int mm;
/* Check if the PCS values match */
for (mm = 0; mm < 2; mm++) {
if (haspcs[mm] == 0)
continue;
for (e = 0; e < nchan; e++) {
double diff;
diff = *((double *)inps[0].c->t[0].fdata[i][pcsix[mm][e]])
- *((double *)inps[0].c->t[0].fdata[k][pcsix[mm][e]]);
if (fabs(diff) > 0.001) {
break;
}
}
if (e < nchan)
break;
}
if (mm < 2)
continue; /* They don't match */
vlist[npat++] = *((double *)inps[0].c->t[0].fdata[k][j]);
valdone[k] = 1;
}
} else {
/* Not device channels */
for (jj = 0; jj < nchan; jj++) {
if (chix[jj] == j)
break;
}
if (jj < nchan)
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;
}
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 Device Geometric Median */
if (inp && dogeom && (haspcs[0] || haspcs[1])) {
double res[3];
/* Locate any patches (including starting patch) with matching PCS values */
npat = 0;
for (k = i; k < inps[0].c->t[0].nsets; k++) {
int mm;
/* Check if the PCS values match */
for (mm = 0; mm < 2; mm++) {
if (haspcs[mm] == 0)
continue;
for (e = 0; e < nchan; e++) {
double diff;
diff = *((double *)inps[0].c->t[0].fdata[i][pcsix[mm][e]])
- *((double *)inps[0].c->t[0].fdata[k][pcsix[mm][e]]);
if (fabs(diff) > 0.001) {
break;
}
}
if (e < nchan)
break;
}
if (mm < 2)
continue; /* They don't match */
/* (nchan == 3 when inp and dogeom) */
v3list[npat][0] = *((double *)inps[0].c->t[0].fdata[k][chix[0]]);
v3list[npat][1] = *((double *)inps[0].c->t[0].fdata[k][chix[1]]);
v3list[npat][2] = *((double *)inps[0].c->t[0].fdata[k][chix[2]]);
npat++;
}
geommed(res, v3list, npat);
*((double *)ocg->t[0].fdata[n][chix[0]]) = res[0];
*((double *)ocg->t[0].fdata[n][chix[1]]) = res[1];
*((double *)ocg->t[0].fdata[n][chix[2]]) = res[2];
/* Override per-component average/median if PCS Geometric Median */
} else if (!inp && 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 {
/* Check/process all the other input files */
for (n = 1; n < ninps; n++) {
/* 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);
for (j = 0; j < inps[0].c->t[0].nfields; j++) {
if (inps[0].c->t[0].ftype[j] != inps[n].c->t[0].ftype[j])
error ("File '%s' field no. %d named '%s' doesn't match file '%s' field '%s'",
inps[n].name, j, inps[n].c->t[0].fsym[j], inps[0].name, inps[0].c->t[0].fsym[j]);
}
/* If merging, append all the values */
if (domerge) {
for (i = 0; i < inps[n].c->t[0].nsets; i++) {
inps[n].c->get_setarr(inps[n].c, 0, i, setel);
ocg->add_setarr(ocg, 0, setel);
}
} else { /* Averaging */
/* 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);
/* For all the patches: */
for (i = 0; i < inps[n].c->t[0].nsets; i++) {
if (inp) {
int k;
/* Check that the PCS values match the first file */
for (k = 0; k < 2; k++) {
if (haspcs[k] == 0)
continue;
for (j = 0; j < 3; j++) {
double diff;
diff = *((double *)inps[0].c->t[0].fdata[i][pcsix[k][j]])
- *((double *)inps[n].c->t[0].fdata[i][pcsix[k][j]]);
if (fabs(diff) > 0.001)
error ("File '%s' set %d has field '%s' value that differs from '%s'",
inps[n].name, i+1, inps[n].c->t[0].fsym[pcsix[k][j]], inps[0].name);
}
}
} else {
/* 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]])
- *((double *)inps[n].c->t[0].fdata[i][chix[j]]);
if (fabs(diff) > 0.001)
error ("File '%s' set %d has field '%s' value that differs from '%s'",
inps[n].name, i+1, inps[n].c->t[0].fsym[chix[j]], inps[0].name);
}
}
}
}
}
/* 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;
if (inp) {
/* Not PCS channels */
int k;
for (k = 0; k < 2; k++) {
if (haspcs[k] == 0)
continue;
for (jj = 0; jj < 3; jj++) {
if (pcsix[k][jj] == j)
break;
}
if (jj < 3)
break;
}
if (k < 2)
continue;
} else {
/* Not device channels */
for (jj = 0; jj < nchan; jj++) {
if (chix[jj] == j)
break;
}
if (jj < nchan)
continue;
}
/* For each patch */
for (i = 0; i < inps[0].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 (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 (inp && dogeom && (haspcs[0] || haspcs[1])) {
double res[3];
/* For each patch */
for (i = 0; i < inps[0].c->t[0].nsets; i++) {
/* For all input files */
npat = 0;
for (n = 0; n < nchan; n++) {
v3list[npat][0] = *((double *)inps[n].c->t[0].fdata[i][chix[0]]);
v3list[npat][1] = *((double *)inps[n].c->t[0].fdata[i][chix[1]]);
v3list[npat][2] = *((double *)inps[n].c->t[0].fdata[i][chix[2]]);
npat++;
}
geommed(res, v3list, npat);
*((double *)ocg->t[0].fdata[i][chix[0]]) = res[0];
*((double *)ocg->t[0].fdata[i][chix[1]]) = res[1];
*((double *)ocg->t[0].fdata[i][chix[2]]) = res[2];
}
/* Override per-component average/median if PCS Geometric Median */
} else if (!inp && dogeom && (haspcs[0] || haspcs[1])) {
double res[3];
/* For each patch */
for (i = 0; i < inps[0].c->t[0].nsets; i++) {
/* For Lab and XYZ */
for (j = 0; j < 2; j++) {
if (haspcs[j] == 0)
continue;
/* 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);
}
}
/* Write out the output and free the cgats * */
for (n = 0; n <= ninps; n++) {
if (n >= ninps) { /* If output file, write it */
if (inps[n].c->write_name(inps[n].c, inps[ninps].name))
error("CGATS file '%s' write error : %s",inps[n].name,inps[n].c->e.m);
}
inps[n].c->del(inps[n].c);
}
free(setel);
free(inps);
return 0;
}
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 or Dev3 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]);
}
Reference in New Issue View Git Blame Copy Permalink