gronod/argyllcms
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 11 Wiki Activity
Files
development
argyllcms/spectro/rspec.c
T
gronod 45bdf51201 Initial commit of v3.5.0
2026-08-20 20:28:40 +01:00

1302 lines
31 KiB
C
Raw Permalink Blame History

/*
* Argyll Color Management System
*
* Author: Graeme W. Gill
* Date: 20015
*
* Copyright 2006 - 2015 Graeme W. Gill
* All rights reserved.
*
* This material is licenced under the GNU GENERAL PUBLIC LICENSE Version 2 or later :-
* see the License2.txt file for licencing details.
*
* Derived from i1pro_imp.c & munki_imp.c
*/
/*
* A library for processing raw spectrometer values.
*
* Currently this is setup for the EX1 spectrometer,
* but the longer term plan is to expand the functionality
* so that it becomes more generic, and can replace a lot
* of common code in i1pro_imp.c & munki_imp.c.
*/
//#define EN_PLOT /* Enable plot support. Set in Jamfile */
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <time.h>
#include <fcntl.h>
#if defined(UNIX)
# include <utime.h>
#else
# include <sys/utime.h>
#endif
#include <sys/stat.h>
#include <stdarg.h>
#ifndef SALONEINSTLIB
#include "copyright.h"
#include "aconfig.h"
#include "numlib.h"
#else /* !SALONEINSTLIB */
#include "sa_config.h"
#include "numsup.h"
#endif /* !SALONEINSTLIB */
# ifndef SALONEINSTLIB
# ifdef EN_PLOT
# include "plot.h"
# endif
# endif
#include "cgats.h"
#include "xspect.h"
#include "insttypes.h"
#include "conv.h"
#include "icoms.h"
#include "inst.h"
#include "rspec.h"
#define BOX_INTEGRATE /* [und] Integrate raw samples as if they were +/-0.5 boxes */
/* (This improves coeficient consistency a bit ?) */
# ifndef SALONEINSTLIB
# ifndef EN_PLOT
# pragma message("###### EN_PLOT is not defined ######")
# endif
# endif
/* -------------------------------------------------- */
#if defined(__APPLE__) && defined(__POWERPC__)
/* Workaround for a ppc gcc 3.3 optimiser bug... */
static int gcc_bug_fix(int i) {
static int nn;
nn += i;
return nn;
}
#endif /* APPLE */
/* -------------------------------------------------- */
/* Setup code */
/* Fit a wavelength polynomial to a set of mapping points */
// ~~~~9999
/* Completely clear an rspec_inf. */
void clear_rspec_inf(rspec_inf *inf) {
memset(inf, 0, sizeof(rspec_inf));
}
/* Completely free contesnt of rspec_inf. */
void free_rspec_inf(rspec_inf *inf) {
if (inf != NULL) {
if (inf->straylight != NULL) {
error("rspec_inf: help - don't know how to free straylight!");
}
if (inf->wlcal)
free(inf->wlcal);
if (inf->findex != NULL)
free(inf->findex);
if (inf->fnocoef != NULL)
free(inf->fnocoef);
if (inf->fcoef != NULL)
free(inf->fcoef);
if (inf->lin != NULL)
free(inf->lin);
if (inf->idark[0] != NULL)
del_rspec(inf->idark[0]);
if (inf->idark[1] != NULL)
del_rspec(inf->idark[1]);
if (inf->ecal != NULL)
free(inf->ecal);
clear_rspec_inf(inf); /* In case it gets reused */
}
}
/* return the number of samples for the given spectral type */
int rspec_typesize(rspec_inf *inf, rspec_type ty) {
int no;
if (ty == rspec_sensor)
no = inf->nsen;
else if (ty == rspec_raw)
no = inf->nraw;
else if (ty == rspec_wav)
no = inf->nwav;
else
error("rspec_typesize type %d unknown",ty);
return no;
}
/* Compute the valid raw range from the calibration information */
void rspec_comp_raw_range_from_ecal(rspec_inf *inf) {
int i;
if (inf->ecaltype != rspec_raw)
error("rspec_comp_raw_range_from_ecal: ecaltype not raw");
for (i = 0; i < inf->nraw; i++) {
if (inf->ecal[i] != 0.0) {
inf->rawrange.off = i;
break;
}
}
if (i >= inf->nraw)
error("rspec_comp_raw_range_from_ecal: ecal is zero");
for (i = inf->rawrange.off; i < inf->nraw; i++) {
if (inf->ecal[i] == 0.0) {
break;
}
}
inf->rawrange.num = i - inf->rawrange.off;
}
/* Convert a raw index to nm using polynomial */
double rspec_raw2nm(rspec_inf *inf, double rix) {
int k;
double wl;
if (inf->nwlcal == 0)
error("rspec_raw2nm: nwlcal == 0");
/* Compute polinomial */
for (wl = inf->wlcal[inf->nwlcal-1], k = inf->nwlcal-2; k >= 0; k--)
wl = wl * rix + inf->wlcal[k];
return wl;
}
/* Convert a cooked index to nm */
double rspec_wav2nm(rspec_inf *inf, double ix) {
return inf->wl_short + ix * inf->wl_space;
}
/* -------------------------------------------------- */
/* Create a new rspec from scratch. */
/* Don't allocate samp if nmeas == 0 */
/* This always succeeds (i.e. application bombs if malloc fails) */
rspec *new_rspec(rspec_inf *inf, rspec_type ty, int nmeas) {
rspec *p;
int no;
if ((p = (rspec *)calloc(1, sizeof(rspec))) == NULL) {
error("Malloc failure in rspec()");
}
p->inf = inf;
p->stype = ty;
p->nmeas = nmeas;
p->nsamp = rspec_typesize(inf, p->stype);
if (nmeas > 0)
p->samp = dmatrix(0, p->nmeas-1, 0, p->nsamp-1);
return p;
}
/* Create a new rspec based on an existing prototype */
/* If nmeas == 0, create space for the same number or measurements */
rspec *new_rspec_proto(rspec *rs, int nmeas) {
rspec *p;
if ((p = (rspec *)calloc(1, sizeof(rspec))) == NULL) {
error("Malloc failure in rspec()");
}
p->inf = rs->inf;
p->stype = rs->stype;
p->mtype = rs->mtype;
p->mcond = rs->mcond;
p->state = rs->state;
p->inttime = rs->inttime;
if (nmeas == 0)
p->nmeas = rs->nmeas;
else
p->nmeas = nmeas;
p->nsamp = rs->nsamp;
p->samp = dmatrix(0, p->nmeas-1, 0, p->nsamp-1);
return p;
}
/* Create a new rspec by cloning an existing one */
rspec *new_rspec_clone(rspec *rs) {
rspec *p;
int i, j;
if ((p = (rspec *)calloc(1, sizeof(rspec))) == NULL) {
error("Malloc failure in rspec()");
}
p->inf = rs->inf;
p->stype = rs->stype;
p->mtype = rs->mtype;
p->mcond = rs->mcond;
p->state = rs->state;
p->inttime = rs->inttime;
p->nmeas = rs->nmeas;
p->nsamp = rs->nsamp;
p->samp = dmatrix(0, p->nmeas-1, 0, p->nsamp-1);
for (i = 0; i < p->nmeas; i++) {
for (j = 0; j < p->nsamp; j++) {
p->samp[i][j] = rs->samp[i][j];
}
}
return p;
}
/* Free a rspec */
void del_rspec(rspec *p) {
if (p != NULL) {
if (p->samp != NULL)
free_dmatrix(p->samp, 0, p->nmeas-1, 0, p->nsamp-1);
free(p);
}
}
/* Plot the first rspec */
void plot_rspec1(rspec *p) {
#ifndef SALONEINSTLIB
# ifdef EN_PLOT
int i, no;
double xx[RSPEC_MAXSAMP];
double yy[RSPEC_MAXSAMP];
no = rspec_typesize(p->inf, p->stype);
for (i = 0; i < no; i++) {
if (p->stype == rspec_wav)
xx[i] = rspec_wav2nm(p->inf, (double)i);
else
xx[i] = (double)i;
yy[i] = p->samp[0][i];
}
do_plot(xx, yy, NULL, NULL, no);
# endif
#endif
}
/* Plot the first rspec of 2 */
void plot_rspec2(rspec *p1, rspec *p2) {
#ifndef SALONEINSTLIB
# ifdef EN_PLOT
int i, no;
double xx[RSPEC_MAXSAMP];
double y1[RSPEC_MAXSAMP];
double y2[RSPEC_MAXSAMP];
// Should check p1 & p2 are compatible ??
no = rspec_typesize(p1->inf, p1->stype);
for (i = 0; i < no; i++) {
if (p1->stype == rspec_wav)
xx[i] = rspec_wav2nm(p1->inf, (double)i);
else
xx[i] = (double)i;
y1[i] = p1->samp[0][i];
y2[i] = p2->samp[0][i];
}
do_plot(xx, y1, y2, NULL, no);
# endif
#endif
}
void plot_ecal(rspec_inf *inf) {
#ifndef SALONEINSTLIB
# ifdef EN_PLOT
int i, no;
double xx[RSPEC_MAXSAMP];
double yy[RSPEC_MAXSAMP];
no = rspec_typesize(inf, inf->ecaltype);
for (i = 0; i < no; i++) {
if (inf->ecaltype == rspec_wav)
xx[i] = rspec_wav2nm(inf, (double)i);
else
xx[i] = (double)i;
yy[i] = inf->ecal[i];
}
do_plot(xx, yy, NULL, NULL, no);
# endif
#endif
}
/* -------------------------------------------------- */
/* Return the largest value */
/* Optionally return the measurement and sample idex of that sample */
double largest_val_rspec(int *pmix, int *psix, rspec *raw) {
double mx = -1e38;
int mi = -1, mj = -1;
int i, j;
if (raw->nmeas <= 0)
error("largest_val_rspec: raw has zero measurements");
for (i = 0; i < raw->nmeas; i++) {
for (j = 0; j < raw->nsamp; j++) {
if (raw->samp[i][j] > mx) {
mx = raw->samp[i][j];
mi = i;
mj = j;
}
}
}
if (pmix != NULL)
*pmix = mi;
if (psix != NULL)
*psix = mj;
return mx;
}
/* return a raw rspec from a sensor rspec */
/* (This does not make any adjustments to the values) */
rspec *extract_raw_from_sensor_rspec(rspec *sens) {
rspec *raw;
int off, i, j;
if (sens->stype != rspec_sensor)
error("extract_raw_from_sensor_rspec: input is not sensor type");
raw = new_rspec(sens->inf, rspec_raw, sens->nmeas);
raw->mtype = sens->mtype;
raw->mcond = sens->mcond;
raw->state = sens->state;
raw->inttime = sens->inttime;
off = sens->inf->lightrange.off;
for (i = 0; i < raw->nmeas; i++) {
for (j = 0; j < raw->nsamp; j++) {
raw->samp[i][j] = sens->samp[i][off + j];
}
}
return raw;
}
/* Return an interpolated dark reference value from idark */
double ex1_interp_idark_val(rspec_inf *inf, int mix, int six, double inttime) {
double idv;
double w0, w1;
int i, j;
w1 = (inttime - inf->idark[0]->inttime)/(inf->idark[1]->inttime - inf->idark[0]->inttime);
w0 = 1.0 - w1;
idv = w0 * inf->idark[0]->samp[mix][six] + w1 * inf->idark[1]->samp[mix][six];
return idv;
}
/* Return an interpolated dark reference from idark */
rspec *ex1_interp_idark(rspec_inf *inf, double inttime) {
double w0, w1;
int i, j;
rspec *dark;
w1 = (inttime - inf->idark[0]->inttime)/(inf->idark[1]->inttime - inf->idark[0]->inttime);
w0 = 1.0 - w1;
dark = new_rspec_proto(inf->idark[0], 0);
for (i = 0; i < inf->idark[0]->nmeas; i++) {
for (j = 0; j < inf->idark[0]->nsamp; j++)
dark->samp[i][j] = w0 * inf->idark[0]->samp[i][j] + w1 * inf->idark[1]->samp[i][j];
}
return dark;
}
/* Subtract the adaptive black */
void subtract_idark_rspec(rspec *raw) {
rspec_inf *inf = raw->inf;
int i, j;
rspec *dark;
if (raw->state & rspec_dcal)
error("subtract_idark_rspec: already done");
if (raw->stype != inf->idark[0]->stype)
error("subtract_idark_rspect: idark does not match rspec type");
dark = ex1_interp_idark(inf, raw->inttime);
for (i = 0; i < raw->nmeas; i++) {
for (j = 0; j < raw->nsamp; j++) {
raw->samp[i][j] -= dark->samp[0][j];
}
}
raw->state |= rspec_dcal;
}
/* Apply non-linearity */
double linearize_val_rspec(rspec_inf *inf, double ival) {
double oval = ival;
int k;
if (ival >= 0.0) {
for (oval = inf->lin[inf->nlin-1], k = inf->nlin-2; k >= 0; k--) {
oval = oval * ival + inf->lin[k];
}
if (inf->lindiv) /* EX1 divides */
oval = ival/oval;
}
return oval;
}
/* Invert non-linearity. */
/* Since the linearisation is nearly a straight line, */
/* a simple Newton inversion will suffice. */
double inv_linearize_val_rspec(rspec_inf *inf, double targv) {
double oval, ival = targv, del = 100.0;
int i, k;
for (i = 0; i < 200 && fabs(del) > 1e-7; i++) {
for (oval = inf->lin[inf->nlin-1], k = inf->nlin-2; k >= 0; k--)
oval = oval * ival + inf->lin[k];
if (inf->lindiv) /* EX1 divides */
oval = ival/oval;
del = (targv - oval);
ival += 0.99 * del;
}
return ival;
}
/* Correct non-linearity */
void linearize_rspec(rspec *raw) {
rspec_inf *inf = raw->inf;
int i, j;
rspec *dark;
if (raw->state & rspec_lin)
error("linearize_rspec: already done");
if (raw->state & rspec_int)
error("linearize_rspec: can't be integration time adjusted");
if (!(raw->state & rspec_dcal))
error("linearize_rspec: needs black subtract");
if (inf->nlin > 0) {
for (i = 0; i < raw->nmeas; i++) {
for (j = 0; j < raw->nsamp; j++) {
raw->samp[i][j] = linearize_val_rspec(inf, raw->samp[i][j]);
}
}
}
raw->state |= rspec_lin;
}
/* Apply the emsissive calibration */
void emis_calibrate_rspec(rspec *raw) {
rspec_inf *inf = raw->inf;
int i, j;
if (raw->state & rspec_cal)
error("emis_calibrate_rspec: already done");
if (raw->stype != raw->inf->ecaltype)
error("emis_calibrate_rspec: ecaltype does not match rspec type");
for (i = 0; i < raw->nmeas; i++) {
for (j = 0; j < raw->nsamp; j++) {
raw->samp[i][j] *= inf->ecal[j];
}
}
raw->state |= rspec_cal;
}
/* Scale to the integration time */
void inttime_calibrate_rspec(rspec *raw) {
rspec_inf *inf = raw->inf;
int i, j;
if (raw->state & rspec_int)
error("inttime_calibrate_rspec: already done");
for (i = 0; i < raw->nmeas; i++) {
for (j = 0; j < raw->nsamp; j++) {
raw->samp[i][j] /= raw->inttime;
}
}
raw->inttime = 1.0;
raw->state |= rspec_int;
}
/* return a wav rspec from a raw rspec */
/* (This does not make any adjustments to the values) */
rspec *convert_wav_from_raw_rspec(rspec *raw) {
rspec_inf *inf = raw->inf;
rspec *wav;
int cx, sx, i, j, k;
if (raw->stype != rspec_raw)
error("extract_raw_from_sensor_rspec: input is not raw type");
wav = new_rspec(raw->inf, rspec_wav, raw->nmeas);
wav->mtype = raw->mtype;
wav->mcond = raw->mcond;
wav->state = raw->state;
wav->inttime = raw->inttime;
for (i = 0; i < wav->nmeas; i++) { /* For each measurement */
for (cx = j = 0; j < inf->nwav; j++) { /* For each wav sample */
double oval = 0.0;
sx = inf->findex[j]; /* Starting index */
for (k = 0; k < inf->fnocoef[j]; k++, cx++, sx++) /* For each matrix value */
oval += inf->fcoef[cx] * raw->samp[i][sx];
wav->samp[i][j] = oval;
}
}
return wav;
}
/* -------------------------------------------------- */
/* Filter code in i1pro_imp is in:
i1pro_compute_wav_filters() X-Rite way
i1pro_create_hr() Using gausian
*/
/* Resampling kernels. (There are more in i1pro_imp.c) */
/* They aren't expected to be unity area, as they will be */
/* normalized anyway. */
/* wi is the width of the filter */
static double triangle(double wi, double x) {
double y = 0.0;
x = fabs(x/wi);
y = 1.0 - x;
if (y < 0.0)
y = 0.0;
return y;
}
static double gausian(double wi, double x) {
double y = 0.0;
wi = wi/(sqrt(2.0 * log(2.0))); /* Convert width at half max to std. dev. */
x = x/wi;
y = exp(-(x * x)); /* Center at 1.0 */
return y;
}
static double lanczos2(double wi, double x) {
double y = 0.0;
wi *= 1.05; // Improves smoothness. Why ?
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;
}
static double lanczos3(double wi, double x) {
double y = 0.0;
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;
}
static double cubicspline(double wi, double x) {
double y = 0.0;
double xx = x;
double bb, cc;
xx = fabs(1.0 * x/wi);
// bb = cc = 1.0/3.0; /* Mitchell */
bb = 0.5;
cc = 0.5;
if (xx < 1.0) {
y = ( 12.0 - 9.0 * bb - 6.0 * cc) * xx * xx * xx
+ (-18.0 + 12.0 * bb + 6.0 * cc) * xx * xx
+ ( 6.0 - 2.0 * bb);
y /= (6.0 - 2.0 * bb);
} else if (xx < 2.0) {
y = ( -1.0 * bb - 6.0 * cc) * xx * xx * xx
+ ( 6.0 * bb + 30.0 * cc) * xx * xx
+ (-12.0 * bb - 48.0 * cc) * xx
+ ( 8.0 * bb + 24.0 * cc);
y /= (6.0 - 2.0 * bb);
} else {
y = 0.0;
}
return y;
}
/* Create the wavelength resampling filters */
void rspec_make_resample_filters(rspec_inf *inf) {
double twidth = inf->wl_space;
double rawspace; /* Average raw band spacing wl */
double fshmax; /* filter shape max wavelength from center */
double finc; /* Integration step size */
int maxcoeffs; /* Maximum coefficients per filter */
int **coeff_ix; /* [band][coef] Raw index */
double **coeff_we; /* [band][coef] Weighting */
double (*kernel)(double wi, double x) = NULL; /* Filter kernel */
int xcount;
int i, j, k;
if (inf->ktype == rspec_triangle)
kernel = triangle;
else if (inf->ktype == rspec_gausian)
kernel = gausian;
else if (inf->ktype == rspec_lanczos2)
kernel = lanczos2;
else if (inf->ktype == rspec_lanczos3)
kernel = lanczos3;
else if (inf->ktype == rspec_cubicspline)
kernel = cubicspline;
else
error("rspec_make_resample_filters: unknown kernel %d",inf->ktype);
#ifdef NEVER // Check kernel sums to 1.0
{
double x, y;
for (x = 0.0; x < 5.0; x += 0.1) {
y = kernel(1.0, x - 4.0)
+ kernel(1.0, x - 3.0)
+ kernel(1.0, x - 2.0)
+ kernel(1.0, x - 1.0)
+ kernel(1.0, x)
+ kernel(1.0, x + 1.0)
+ kernel(1.0, x + 2.0);
+ kernel(1.0, x + 3.0);
+ kernel(1.0, x + 4.0);
printf("Offset %f sum %f\n",x,y);
}
}
#endif // NEVER
/* Aproximate raw value spacing in nm */
rawspace = (inf->wl_long - inf->wl_short)/inf->rawrange.num;
//printf("~1 rawspace = %f\n",rawspace);
/* Figure the extent of the filter kernel. We assume they */
/* all have a finite extent. */
for (fshmax = 50.0; fshmax >= 0.0; fshmax -= 0.01) {
if (fabs(kernel(twidth, fshmax)) > 1e-6) {
fshmax += 0.01;
break;
}
}
//printf("~1 fshmax = %f\n",fshmax);
if (fshmax <= 0.0)
error("rspec_make_resample_filters: fshmax search failed\n");
a1logd(inf->log, 4,"rspec_make_resample_filters: fshmax = %f\n",fshmax);
/* Figure number of raw samples over kernel extent. */
/* (Allow generous factor for non-linearity) */
maxcoeffs = (int)floor(2.0 * 1.4 * fshmax/rawspace + 3.0);
a1logd(inf->log, 4,"rspec_make_resample_filters: maxcoeffs = %d\n",maxcoeffs);
/* Figure out box integration step size */
#ifdef FAST_HIGH_RES_SETUP
finc = twidth/50.0;
if (rawspace/finc < 10.0)
finc = rawspace/10.0;
#else
finc = twidth/15.0;
if (rawspace/finc < 4.0)
finc = rawspace/4.0;
#endif
a1logd(inf->log, 4,"rspec_make_resample_filters: integration step = %f\n",finc);
if (inf->fnocoef != NULL)
free(inf->fnocoef);
if ((inf->fnocoef = (int *)calloc(inf->nwav, sizeof(int))) == NULL)
error("rspec_make_resample_filters: malloc failure");
/* Space to build filter coeficients */
coeff_ix = imatrix(0, inf->nwav-1, 0, maxcoeffs-1);
coeff_we = dmatrix(0, inf->nwav-1, 0, maxcoeffs-1);
/* For all the usable raw bands */
for (i = inf->rawrange.off+1; i < (inf->rawrange.off+inf->rawrange.num-1); i++) {
double w1, wl, w2;
/* Translate CCD center and boundaries to calibrated wavelength */
wl = rspec_raw2nm(inf, (double)i);
w1 = rspec_raw2nm(inf, (double)i - 0.5);
w2 = rspec_raw2nm(inf, (double)i + 0.5);
// printf("~1 CCD %d, w1 %f, wl %f, w2 %f\n",i,w1,wl,w2);
/* For each output filter */
for (j = 0; j < inf->nwav; j++) {
double cwl, rwl; /* center, relative wavelegth */
double we;
cwl = rspec_wav2nm(inf, (double)j);
rwl = wl - cwl; /* raw relative wavelength to filter */
if (fabs(w1 - cwl) > fshmax && fabs(w2 - cwl) > fshmax)
continue; /* Doesn't fall into this filter */
#ifdef BOX_INTEGRATE
/* Integrate in finc nm increments from filter shape */
/* using triangular integration. */
{
int nn;
double lw, ll;
nn = (int)(fabs(w2 - w1)/finc + 0.5); /* Number to integrate over */
lw = w1; /* start at lower boundary of CCD cell */
ll = kernel(twidth, w1 - cwl);
we = 0.0;
for (k = 0; k < nn; k++) {
double cw, cl;
#if defined(__APPLE__) && defined(__POWERPC__)
gcc_bug_fix(k);
#endif
cw = w1 + (k+1.0)/(nn + 1.0) * fabs(w2 - w1); /* wl to sample */
cl = kernel(twidth, cw - cwl);
we += 0.5 * (cl + ll) * fabs(lw - cw); /* Area under triangle */
ll = cl;
lw = cw;
}
}
#else
we = fabs(w2 - w1) * kernel(twidth, rwl);
#endif
if (inf->fnocoef[j] >= maxcoeffs)
error("rspec_make_resample_filters: run out of high res filter space\n");
coeff_ix[j][inf->fnocoef[j]] = i;
coeff_we[j][inf->fnocoef[j]++] = we;
// printf("~1 filter %d, cwl %f, rwl %f, ix %d, we %f, nocoefs %d\n",j,cwl,rwl,i,we,info->fnocoef[j]);
}
}
/* Convert hires filters into runtime format: */
/* Allocate or reallocate high res filter tables */
if (inf->findex != NULL)
free(inf->findex);
if (inf->fcoef != NULL)
free(inf->fcoef);
if ((inf->findex = (int *)calloc(inf->nraw, sizeof(int))) == NULL)
error("rspec_make_resample_filters: malloc index failed!\n");
/* Count the total number of coefficients */
for (xcount = j = 0; j < inf->nwav; j++) {
inf->findex[j] = coeff_ix[j][0]; /* raw starting index */
xcount += inf->fnocoef[j];
}
//printf("~1 total coefs = %d\n",xcount);
/* Allocate space for them */
if ((inf->fcoef = (double *)calloc(xcount, sizeof(double))) == NULL)
error("rspec_make_resample_filters: malloc index failed!\n");
/* Normalize the weight * nm to 1.0, and pack them into the run-time format */
for (i = j = 0; j < inf->nwav; j++) {
int sx;
double rwi, twe = 0.0;
sx = inf->findex[j]; /* raw starting index */
for (k = 0; k < inf->fnocoef[j]; sx++, k++) {
/* Width of raw band in nm */
rwi = fabs(rspec_raw2nm(inf, (double)sx - 0.5)
- rspec_raw2nm(inf, (double)sx + 0.5));
twe += rwi * coeff_we[j][k];
}
if (twe > 0.0)
twe = 1.0/twe;
else
twe = 1.0;
// printf("Output %d, nocoefs %d, norm weight %f:\n",j,inf->fnocoef[j],twe);
for (k = 0; k < inf->fnocoef[j]; k++, i++) {
inf->fcoef[i] = coeff_we[j][k] * twe;
// printf(" coef %d packed %d from raw %d val %f\n",k,i,inf->findex[j]+k,inf->fcoef[i]);
}
}
free_imatrix(coeff_ix, 0, inf->nwav-1, 0, maxcoeffs-1);
free_dmatrix(coeff_we, 0, inf->nwav-1, 0, maxcoeffs-1);
}
//printf("~1 line %d\n",__LINE__);
/* Plot the wave resampling filters */
void plot_resample_filters(rspec_inf *inf) {
#ifndef SALONEINSTLIB
# ifdef EN_PLOT
double *xx, *ss;
double **yy;
int i, j, k, sx;
//printf("~1 nraw = %d\n",inf->nraw);
xx = dvectorz(0, inf->nraw-1); /* X index */
yy = dmatrixz(0, 5, 0, inf->nraw-1); /* Curves distributed amongst 5 graphs */
/* with 6th holding sum */
for (i = 0; i < inf->nraw; i++)
xx[i] = i;
/* For each output wavelength */
for (i = j = 0; j < inf->nwav; j++) {
sx = inf->findex[j]; /* raw starting index */
//printf("Output %d raw sx %d\n",j,sx);
/* For each matrix value */
for (k = 0; k < inf->fnocoef[j]; k++, sx++, i++) {
yy[5][sx] += 0.5 * inf->fcoef[i];
yy[j % 5][sx] = inf->fcoef[i];
//printf(" filter %d six %d weight = %e\n",k,sx,inf->fcoef[i]);
}
}
printf("Wavelength re-sampling curves:\n");
// do_plot6(xx, yy[0], yy[1], yy[2], yy[3], yy[4], yy[5], 150);
do_plot6(xx, yy[0], yy[1], yy[2], yy[3], yy[4], yy[5], inf->nraw);
free_dvector(xx, 0, inf->nraw-1);
free_dmatrix(yy, 0, 2, 0, inf->nraw-1);
# endif
#endif
}
/* ================================================== */
/* Calibration file support */
/* Open the file. nz wr for write mode, else read */
/* Return nz on error */
int calf_open(calf *x, a1log *log, char *fname, int wr) {
char nmode[10];
char cal_name[200];
char **cal_paths = NULL;
int no_paths = 0;
memset((void *)x, 0, sizeof(calf));
x->log = log;
if (wr)
strcpy(nmode, "w");
else
strcpy(nmode, "r");
#if !defined(O_CREAT) && !defined(_O_CREAT)
# error "Need to #include fcntl.h!"
#endif
#if defined(O_BINARY) || defined(_O_BINARY)
strcat(nmode, "b");
#endif
/* Create the file name */
if (wr)
sprintf(cal_name, "ArgyllCMS/%s", fname);
else
sprintf(cal_name, "ArgyllCMS/%s" SSEPS "color/%s", fname, fname);
if ((no_paths = xdg_bds(NULL, &cal_paths, xdg_cache, xdg_write, xdg_user, xdg_none,
cal_name)) < 1) {
a1logd(x->log,1,"calf_open: xdg_bds returned no paths\n");
return 1;
}
a1logd(x->log,2,"calf_open: %s file '%s'\n",cal_paths[0], wr ? "saving to" : "restoring from");
/* Check the last modification time */
if (!wr) {
struct sys_stat sbuf;
if (sys_stat(cal_paths[0], &sbuf) == 0) {
x->lo_secs = time(NULL) - sbuf.st_mtime;
a1logd(x->log,2,"calf_open:: %d secs from instrument last open\n",x->lo_secs);
} else {
a1logd(x->log,2,"calf_open:: stat on file failed\n");
}
}
if ((wr && create_parent_directories(cal_paths[0]))
|| (x->fp = fopen(cal_paths[0], nmode)) == NULL) {
a1logd(x->log,2,"calf_open: failed to open file for %s\n",wr ? "writing" : "reading");
xdg_free(cal_paths, no_paths);
return 1;
}
xdg_free(cal_paths, no_paths);
a1logd(x->log,2,"calf_open: succeeded\n");
return 0;
}
/* Update the modification time */
/* Return nz on error */
int calf_touch(a1log *log, char *fname) {
char cal_name[200];
char **cal_paths = NULL;
int no_paths = 0;
int rv;
/* Locate the file name */
sprintf(cal_name, "ArgyllCMS/%s" SSEPS "color/%s", fname, fname);
if ((no_paths = xdg_bds(NULL, &cal_paths, xdg_cache, xdg_read, xdg_user, xdg_none,
cal_name)) < 1) {
a1logd(log,2,"calf_touch: xdg_bds failed to locate file'\n");
return 1;
}
a1logd(log,2,"calf_touch: touching file '%s'\n",cal_paths[0]);
if ((rv = sys_utime(cal_paths[0], NULL)) != 0) {
a1logd(log,2,"calf_touch: failed with %d\n",rv);
xdg_free(cal_paths, no_paths);
return 1;
}
xdg_free(cal_paths, no_paths);
return 0;
}
/* Rewind and reset for another read */
void calf_rewind(calf *x) {
x->ef = 0;
x->chsum = 0;
x->nbytes = 0;
rewind(x->fp);
}
/* Close the file and free any memory */
/* return nz on error */
int calf_done(calf *x) {
int rv = 0;
if (x->fp != NULL) {
if (fclose(x->fp)) {
a1logd(x->log,2,"calf_done: closing file failed\n");
rv = 1;
}
}
if (x->buf != NULL)
free(x->buf);
x->buf = NULL;
return rv;
}
static void sizebuf(calf *x, size_t size) {
if (x->bufsz < size)
x->buf = realloc(x->buf, size);
if (x->buf == NULL)
error("calf: sizebuf malloc failed");
}
static void update_chsum(calf *x, unsigned char *p, int nn) {
int i;
for (i = 0; i < nn; i++, p++)
x->chsum = ((x->chsum << 13) | (x->chsum >> (32-13))) + *p;
x->nbytes += nn;
}
/* Write an array of ints to the file. Set the error flag to nz on error */
void calf_wints(calf *x, int *dp, int n) {
if (x->ef)
return;
if (fwrite((void *)dp, sizeof(int), n, x->fp) != n) {
x->ef = 1;
a1logd(x->log,2,"calf_wints: write failed for %d ints at offset %d\n",n,x->nbytes);
} else {
update_chsum(x, (unsigned char *)dp, sizeof(int) * n);
}
}
/* Write an array of doubles to the file. Set the error flag to nz on error */
void calf_wdoubles(calf *x, double *dp, int n) {
if (x->ef)
return;
if (fwrite((void *)dp, sizeof(double), n, x->fp) != n) {
x->ef = 1;
a1logd(x->log,2,"calf_wdoubles: write failed for %d doubles at offset %d\n",n,x->nbytes);
} else {
update_chsum(x, (unsigned char *)dp, sizeof(double) * n);
}
}
/* Write an array of time_t's to the file. Set the error flag to nz on error */
/* (This will cause file checksum fail if different executables on the same */
/* system have different time_t values) */
void calf_wtime_ts(calf *x, time_t *dp, int n) {
if (x->ef)
return;
if (fwrite((void *)dp, sizeof(time_t), n, x->fp) != n) {
x->ef = 1;
a1logd(x->log,2,"calf_wtime_ts: write failed for %d time_ts at offset %d\n",n,x->nbytes);
} else {
update_chsum(x, (unsigned char *)dp, sizeof(time_t) * n);
}
}
/* Write a zero terminated string */
void calf_wstrz(calf *x, char *dp) {
int n;
if (x->ef)
return;
n = strlen(dp) + 1;
calf_wints(x, &n, 1);
if (fwrite((void *)dp, sizeof(char), n, x->fp) != n) {
x->ef = 1;
a1logd(x->log,2,"calf_wstrz: write failed for %d long string at offset %d\n",n,x->nbytes);
} else {
update_chsum(x, (unsigned char *)dp, sizeof(char) * n);
}
}
/* Read an array of ints from the file. Set the error flag to nz on error */
/* Always read (ignore rd flag) */
void calf_rints2(calf *x, int *dp, int n) {
if (x->ef)
return;
if (fread((void *)dp, sizeof(int), n, x->fp) != n) {
x->ef = 1;
a1logd(x->log,2,"calf_rints2: read failed for %d ints at offset %d\n",n,x->nbytes);
} else {
update_chsum(x, (unsigned char *)dp, sizeof(int) * n);
}
}
/* Read an array of ints from the file. Set the error flag to nz on error */
void calf_rints(calf *x, int *dp, int n) {
size_t nbytes = n * sizeof(int);
void *dest = (void *)dp;
if (x->ef)
return;
if (x->rd == 0) { /* Dummy read */
sizebuf(x, nbytes);
dest = x->buf;
}
if (fread(dest, 1, nbytes, x->fp) != nbytes) {
x->ef = 1;
a1logd(x->log,2,"calf_rints: read failed for %d ints at offset %d\n",n,x->nbytes);
} else {
update_chsum(x, dest, nbytes);
}
}
/* Read an array of doubles from the file. Set the error flag to nz on error */
void calf_rdoubles(calf *x, double *dp, int n) {
size_t nbytes = n * sizeof(double);
void *dest = (void *)dp;
if (x->ef)
return;
if (x->rd == 0) { /* Dummy read */
sizebuf(x, nbytes);
dest = x->buf;
}
if (fread(dest, 1, nbytes, x->fp) != nbytes) {
x->ef = 1;
a1logd(x->log,2,"calf_rdoubles: read failed for %d ints at offset %d\n",n,x->nbytes);
} else {
update_chsum(x, dest, nbytes);
}
}
/* Read an array of time_t's from the file. Set the error flag to nz on error */
/* (This will cause file checksum fail if different executables on the same */
/* system have different time_t values) */
void calf_rtime_ts(calf *x, time_t *dp, int n) {
size_t nbytes = n * sizeof(time_t);
void *dest = (void *)dp;
if (x->ef)
return;
if (x->rd == 0) { /* Dummy read */
sizebuf(x, nbytes);
dest = x->buf;
}
if (fread(dest, 1, nbytes, x->fp) != nbytes) {
x->ef = 1;
a1logd(x->log,2,"calf_rtime_ts: read failed for %d ints at offset %d\n",n,x->nbytes);
} else {
update_chsum(x, dest, nbytes);
}
}
/* Read a zero terminated string. */
void calf_rstrz(calf *x, char **dp) {
int n;
size_t nbytes = 0;
char *dest = NULL;
if (x->ef)
return;
calf_rints2(x, &n, 1);
nbytes = sizeof(char) * n;
if (x->ef || n == 0)
return;
if (x->rd != 0) { /* Reading for real */
if (*dp != NULL)
free(*dp);
if ((*dp = dest = malloc(nbytes)) == NULL)
error("calf: calf_rstrz malloc failed");
} else {
sizebuf(x, nbytes);
dest = x->buf;
}
if (fread(dest, 1, nbytes, x->fp) != nbytes) {
x->ef = 1;
a1logd(x->log,2,"calf_rstrz: read failed for %d long string at offset %d\n",n,x->nbytes);
} else {
update_chsum(x, (unsigned char*)dest, nbytes);
}
}
void calf_rstrz2(calf *x, char **dp) {
int rd = x->rd;
x->rd = 1;
calf_rstrz(x, dp);
x->rd = rd;
}
/* ================================================== */
/* Save a rspec to a calibration file */
void calf_wrspec(calf *x, rspec *s) {
int i;
if (x->ef)
return;
calf_wints(x, (int *)&s->stype, 1);
calf_wints(x, (int *)&s->mtype, 1);
calf_wints(x, (int *)&s->mcond, 1);
calf_wints(x, (int *)&s->state, 1);
calf_wdoubles(x, &s->inttime, 1);
calf_wints(x, &s->nmeas, 1);
calf_wints(x, &s->nsamp, 1);
for (i = 0; i < s->nmeas; i++) {
calf_wdoubles(x, s->samp[i], s->nsamp);
}
}
/* Create a rspec from a calibration file */
void calf_rrspec(calf *x, rspec **dp, rspec_inf *inf) {
rspec *s, dumy;
int no, i;
if (x->ef)
return;
if (x->rd != 0) {
if (*dp != NULL)
del_rspec(*dp);
*dp = s = new_rspec(inf, rspec_sensor, 0);
} else {
s = &dumy;
}
calf_rints2(x, (int *)&s->stype, 1);
calf_rints2(x, (int *)&s->mtype, 1);
calf_rints2(x, (int *)&s->mcond, 1);
calf_rints2(x, (int *)&s->state, 1);
calf_rdoubles(x, &s->inttime, 1);
calf_rints2(x, &s->nmeas, 1);
calf_rints2(x, &s->nsamp, 1);
/* Sanity check. */
no = rspec_typesize(inf, s->stype);
if (no != s->nsamp) {
a1logd(inf->log, 4,"calf_rrspec: unexpected nsamp %d (expect %d)\n",s->nsamp,no);
x->ef = 1;
return;
}
if (x->rd != 0) {
s->samp = dmatrix(0, s->nmeas-1, 0, s->nsamp-1);
for (i = 0; i < s->nmeas; i++) {
calf_rdoubles(x, s->samp[i], s->nsamp);
}
} else {
for (i = 0; i < s->nmeas; i++) {
calf_rdoubles(x, NULL, s->nsamp);
}
}
}
Reference in New Issue View Git Blame Copy Permalink