1
2 /*
3 * Argyll Color Correction System
4 *
5 * Author: Graeme W. Gill
6 * Date: 20015
7 *
8 * Copyright 2006 - 2015 Graeme W. Gill
9 * All rights reserved.
10 *
11 * This material is licenced under the GNU GENERAL PUBLIC LICENSE Version 2 or later :-
12 * see the License2.txt file for licencing details.
13 *
14 * Derived from i1pro_imp.c & munki_imp.c
15 */
16
17 /*
18 * A library for processing raw spectrometer values.
19 *
20 * Currently this is setup for the EX1 spectrometer,
21 * but the longer term plan is to expand the functionality
22 * so that it becomes more generic, and can replace a lot
23 * of common code in i1pro_imp.c & munki_imp.c.
24 */
25
26 #include <stdio.h>
27 #include <stdlib.h>
28 #include <ctype.h>
29 #include <string.h>
30 #include <time.h>
31 #include <fcntl.h>
32 #if defined(UNIX)
33 # include <utime.h>
34 #else
35 # include <sys/utime.h>
36 #endif
37 #include <sys/stat.h>
38 #include <stdarg.h>
39 #ifndef SALONEINSTLIB
40 #include "copyright.h"
41 #include "aconfig.h"
42 #include "numlib.h"
43 #else /* !SALONEINSTLIB */
44 #include "sa_config.h"
45 #include "numsup.h"
46 #endif /* !SALONEINSTLIB */
47 #ifndef SALONEINSTLIB
48 # include "plot.h"
49 #endif
50 #include "xspect.h"
51 #include "insttypes.h"
52 #include "conv.h"
53 #include "icoms.h"
54 #include "inst.h"
55 #include "rspec.h"
56
57 #define BOX_INTEGRATE /* [und] Integrate raw samples as if they were +/-0.5 boxes */
58 /* (This improves coeficient consistency a bit ?) */
59
60 /* -------------------------------------------------- */
61 #if defined(__APPLE__) && defined(__POWERPC__)
62
63 /* Workaround for a ppc gcc 3.3 optimiser bug... */
gcc_bug_fix(int i)64 static int gcc_bug_fix(int i) {
65 static int nn;
66 nn += i;
67 return nn;
68 }
69 #endif /* APPLE */
70
71 /* -------------------------------------------------- */
72 /* Setup code */
73
74 /* Fit a wavelength polynomial to a set of mapping points */
75 // ~~~~9999
76
77 /* Completely clear an rspec_inf. */
clear_rspec_inf(rspec_inf * inf)78 void clear_rspec_inf(rspec_inf *inf) {
79 memset(inf, 0, sizeof(rspec_inf));
80 }
81
82 /* Completely free contesnt of rspec_inf. */
free_rspec_inf(rspec_inf * inf)83 void free_rspec_inf(rspec_inf *inf) {
84
85 if (inf != NULL) {
86 if (inf->straylight != NULL) {
87 error("rspec_inf: help - don't know how to free straylight!");
88 }
89
90 if (inf->wlcal)
91 free(inf->wlcal);
92
93 if (inf->findex != NULL)
94 free(inf->findex);
95 if (inf->fnocoef != NULL)
96 free(inf->fnocoef);
97 if (inf->fcoef != NULL)
98 free(inf->fcoef);
99
100 if (inf->lin != NULL)
101 free(inf->lin);
102
103 if (inf->idark[0] != NULL)
104 del_rspec(inf->idark[0]);
105 if (inf->idark[1] != NULL)
106 del_rspec(inf->idark[1]);
107
108 if (inf->ecal != NULL)
109 free(inf->ecal);
110
111 clear_rspec_inf(inf); /* In case it gets reused */
112 }
113 }
114
115
116 /* return the number of samples for the given spectral type */
rspec_typesize(rspec_inf * inf,rspec_type ty)117 int rspec_typesize(rspec_inf *inf, rspec_type ty) {
118 int no;
119 if (ty == rspec_sensor)
120 no = inf->nsen;
121 else if (ty == rspec_raw)
122 no = inf->nraw;
123 else if (ty == rspec_wav)
124 no = inf->nwav;
125 else
126 error("rspec_typesize type %d unknown",ty);
127 return no;
128 }
129
130 /* Compute the valid raw range from the calibration information */
rspec_comp_raw_range_from_ecal(rspec_inf * inf)131 void rspec_comp_raw_range_from_ecal(rspec_inf *inf) {
132 int i;
133
134 if (inf->ecaltype != rspec_raw)
135 error("rspec_comp_raw_range_from_ecal: ecaltype not raw");
136
137 for (i = 0; i < inf->nraw; i++) {
138 if (inf->ecal[i] != 0.0) {
139 inf->rawrange.off = i;
140 break;
141 }
142 }
143 if (i >= inf->nraw)
144 error("rspec_comp_raw_range_from_ecal: ecal is zero");
145
146 for (i = inf->rawrange.off; i < inf->nraw; i++) {
147 if (inf->ecal[i] == 0.0) {
148 break;
149 }
150 }
151 inf->rawrange.num = i - inf->rawrange.off;
152 }
153
154 /* Convert a raw index to nm using polynomial */
rspec_raw2nm(rspec_inf * inf,double rix)155 double rspec_raw2nm(rspec_inf *inf, double rix) {
156 int k;
157 double wl;
158
159 if (inf->nwlcal == 0)
160 error("rspec_raw2nm: nwlcal == 0");
161
162 /* Compute polinomial */
163 for (wl = inf->wlcal[inf->nwlcal-1], k = inf->nwlcal-2; k >= 0; k--)
164 wl = wl * rix + inf->wlcal[k];
165
166 return wl;
167 }
168
169 /* Convert a cooked index to nm */
rspec_wav2nm(rspec_inf * inf,double ix)170 double rspec_wav2nm(rspec_inf *inf, double ix) {
171 return inf->wl_short + ix * inf->wl_space;
172 }
173
174 /* -------------------------------------------------- */
175
176 /* Create a new rspec from scratch. */
177 /* Don't allocate samp if nmeas == 0 */
178 /* This always succeeds (i.e. application bombs if malloc fails) */
new_rspec(rspec_inf * inf,rspec_type ty,int nmeas)179 rspec *new_rspec(rspec_inf *inf, rspec_type ty, int nmeas) {
180 rspec *p;
181 int no;
182
183 if ((p = (rspec *)calloc(1, sizeof(rspec))) == NULL) {
184 error("Malloc failure in rspec()");
185 }
186 p->inf = inf;
187 p->stype = ty;
188
189 p->nmeas = nmeas;
190 p->nsamp = rspec_typesize(inf, p->stype);
191 if (nmeas > 0)
192 p->samp = dmatrix(0, p->nmeas-1, 0, p->nsamp-1);
193
194 return p;
195 }
196
197 /* Create a new rspec based on an existing prototype */
198 /* If nmeas == 0, create space for the same number or measurements */
new_rspec_proto(rspec * rs,int nmeas)199 rspec *new_rspec_proto(rspec *rs, int nmeas) {
200 rspec *p;
201
202 if ((p = (rspec *)calloc(1, sizeof(rspec))) == NULL) {
203 error("Malloc failure in rspec()");
204 }
205 p->inf = rs->inf;
206 p->stype = rs->stype;
207 p->mtype = rs->mtype;
208 p->state = rs->state;
209 p->inttime = rs->inttime;
210
211 if (nmeas == 0)
212 p->nmeas = rs->nmeas;
213 else
214 p->nmeas = nmeas;
215 p->nsamp = rs->nsamp;
216 p->samp = dmatrix(0, p->nmeas-1, 0, p->nsamp-1);
217
218 return p;
219 }
220
221 /* Create a new rspec by cloning an existing one */
new_rspec_clone(rspec * rs)222 rspec *new_rspec_clone(rspec *rs) {
223 rspec *p;
224 int i, j;
225
226 if ((p = (rspec *)calloc(1, sizeof(rspec))) == NULL) {
227 error("Malloc failure in rspec()");
228 }
229 p->inf = rs->inf;
230 p->stype = rs->stype;
231 p->mtype = rs->mtype;
232 p->state = rs->state;
233 p->inttime = rs->inttime;
234
235 p->nmeas = rs->nmeas;
236 p->nsamp = rs->nsamp;
237 p->samp = dmatrix(0, p->nmeas-1, 0, p->nsamp-1);
238
239 for (i = 0; i < p->nmeas; i++) {
240 for (j = 0; j < p->nsamp; j++) {
241 p->samp[i][j] = rs->samp[i][j];
242 }
243 }
244
245 return p;
246 }
247
248 /* Free a rspec */
del_rspec(rspec * p)249 void del_rspec(rspec *p) {
250 if (p != NULL) {
251 if (p->samp != NULL)
252 free_dmatrix(p->samp, 0, p->nmeas-1, 0, p->nsamp-1);
253 free(p);
254 }
255 }
256
257 /* Plot the first rspec */
plot_rspec1(rspec * p)258 void plot_rspec1(rspec *p) {
259 #ifndef SALONEINSTLIB
260 int i, no;
261 double xx[RSPEC_MAXSAMP];
262 double yy[RSPEC_MAXSAMP];
263
264 no = rspec_typesize(p->inf, p->stype);
265
266 for (i = 0; i < no; i++) {
267 if (p->stype == rspec_wav)
268 xx[i] = rspec_wav2nm(p->inf, (double)i);
269 else
270 xx[i] = (double)i;
271 yy[i] = p->samp[0][i];
272 }
273 do_plot(xx, yy, NULL, NULL, no);
274 #endif
275 }
276
277 /* Plot the first rspec of 2 */
plot_rspec2(rspec * p1,rspec * p2)278 void plot_rspec2(rspec *p1, rspec *p2) {
279 #ifndef SALONEINSTLIB
280 int i, no;
281 double xx[RSPEC_MAXSAMP];
282 double y1[RSPEC_MAXSAMP];
283 double y2[RSPEC_MAXSAMP];
284
285 // Should check p1 & p2 are compatible ??
286
287 no = rspec_typesize(p1->inf, p1->stype);
288
289 for (i = 0; i < no; i++) {
290 if (p1->stype == rspec_wav)
291 xx[i] = rspec_wav2nm(p1->inf, (double)i);
292 else
293 xx[i] = (double)i;
294 y1[i] = p1->samp[0][i];
295 y2[i] = p2->samp[0][i];
296 }
297 do_plot(xx, y1, y2, NULL, no);
298 #endif
299 }
300
plot_ecal(rspec_inf * inf)301 void plot_ecal(rspec_inf *inf) {
302 #ifndef SALONEINSTLIB
303 int i, no;
304 double xx[RSPEC_MAXSAMP];
305 double yy[RSPEC_MAXSAMP];
306
307 no = rspec_typesize(inf, inf->ecaltype);
308
309 for (i = 0; i < no; i++) {
310 if (inf->ecaltype == rspec_wav)
311 xx[i] = rspec_wav2nm(inf, (double)i);
312 else
313 xx[i] = (double)i;
314 yy[i] = inf->ecal[i];
315 }
316 do_plot(xx, yy, NULL, NULL, no);
317 #endif
318 }
319
320
321 /* -------------------------------------------------- */
322
323 /* Return the largest value */
324 /* Optionally return the measurement and sample idex of that sample */
largest_val_rspec(int * pmix,int * psix,rspec * raw)325 double largest_val_rspec(int *pmix, int *psix, rspec *raw) {
326 double mx = -1e38;
327 int mi = -1, mj = -1;
328 int i, j;
329
330 if (raw->nmeas <= 0)
331 error("largest_val_rspec: raw has zero measurements");
332
333 for (i = 0; i < raw->nmeas; i++) {
334 for (j = 0; j < raw->nsamp; j++) {
335 if (raw->samp[i][j] > mx) {
336 mx = raw->samp[i][j];
337 mi = i;
338 mj = j;
339 }
340 }
341 }
342 if (pmix != NULL)
343 *pmix = mi;
344 if (psix != NULL)
345 *psix = mj;
346
347 return mx;
348 }
349
350 /* return a raw rspec from a sensor rspec */
351 /* (This does not make any adjustments to the values) */
extract_raw_from_sensor_rspec(rspec * sens)352 rspec *extract_raw_from_sensor_rspec(rspec *sens) {
353 rspec *raw;
354 int off, i, j;
355
356 if (sens->stype != rspec_sensor)
357 error("extract_raw_from_sensor_rspec: input is not sensor type");
358
359 raw = new_rspec(sens->inf, rspec_raw, sens->nmeas);
360
361 raw->mtype = sens->mtype;
362 raw->state = sens->state;
363 raw->inttime = sens->inttime;
364
365 off = sens->inf->lightrange.off;
366 for (i = 0; i < raw->nmeas; i++) {
367 for (j = 0; j < raw->nsamp; j++) {
368 raw->samp[i][j] = sens->samp[i][off + j];
369 }
370 }
371
372 return raw;
373 }
374
375 /* Return an interpolated dark reference value from idark */
ex1_interp_idark_val(rspec_inf * inf,int mix,int six,double inttime)376 double ex1_interp_idark_val(rspec_inf *inf, int mix, int six, double inttime) {
377 double idv;
378 double w0, w1;
379 int i, j;
380
381 w1 = (inttime - inf->idark[0]->inttime)/(inf->idark[1]->inttime - inf->idark[0]->inttime);
382 w0 = 1.0 - w1;
383
384 idv = w0 * inf->idark[0]->samp[mix][six] + w1 * inf->idark[1]->samp[mix][six];
385
386 return idv;
387 }
388
389 /* Return an interpolated dark reference from idark */
ex1_interp_idark(rspec_inf * inf,double inttime)390 rspec *ex1_interp_idark(rspec_inf *inf, double inttime) {
391 double w0, w1;
392 int i, j;
393 rspec *dark;
394
395 w1 = (inttime - inf->idark[0]->inttime)/(inf->idark[1]->inttime - inf->idark[0]->inttime);
396 w0 = 1.0 - w1;
397
398 dark = new_rspec_proto(inf->idark[0], 0);
399
400 for (i = 0; i < inf->idark[0]->nmeas; i++) {
401 for (j = 0; j < inf->idark[0]->nsamp; j++)
402 dark->samp[i][j] = w0 * inf->idark[0]->samp[i][j] + w1 * inf->idark[1]->samp[i][j];
403 }
404
405 return dark;
406 }
407
408 /* Subtract the adaptive black */
subtract_idark_rspec(rspec * raw)409 void subtract_idark_rspec(rspec *raw) {
410 rspec_inf *inf = raw->inf;
411 int i, j;
412 rspec *dark;
413
414 if (raw->state & rspec_dcal)
415 error("subtract_idark_rspec: already done");
416
417 if (raw->stype != inf->idark[0]->stype)
418 error("subtract_idark_rspect: idark does not match rspec type");
419
420 dark = ex1_interp_idark(inf, raw->inttime);
421
422 for (i = 0; i < raw->nmeas; i++) {
423 for (j = 0; j < raw->nsamp; j++) {
424 raw->samp[i][j] -= dark->samp[0][j];
425 }
426 }
427
428 raw->state |= rspec_dcal;
429 }
430
431 /* Apply non-linearity */
linearize_val_rspec(rspec_inf * inf,double ival)432 double linearize_val_rspec(rspec_inf *inf, double ival) {
433 double oval = ival;
434 int k;
435
436 if (ival >= 0.0) {
437 for (oval = inf->lin[inf->nlin-1], k = inf->nlin-2; k >= 0; k--) {
438 oval = oval * ival + inf->lin[k];
439 }
440
441 if (inf->lindiv) /* EX1 divides */
442 oval = ival/oval;
443 }
444 return oval;
445 }
446
447 /* Invert non-linearity. */
448 /* Since the linearisation is nearly a straight line, */
449 /* a simple Newton inversion will suffice. */
inv_linearize_val_rspec(rspec_inf * inf,double targv)450 double inv_linearize_val_rspec(rspec_inf *inf, double targv) {
451 double oval, ival = targv, del = 100.0;
452 int i, k;
453
454 for (i = 0; i < 200 && fabs(del) > 1e-7; i++) {
455 for (oval = inf->lin[inf->nlin-1], k = inf->nlin-2; k >= 0; k--)
456 oval = oval * ival + inf->lin[k];
457 if (inf->lindiv) /* EX1 divides */
458 oval = ival/oval;
459
460 del = (targv - oval);
461 ival += 0.99 * del;
462 }
463 return ival;
464 }
465
466 /* Correct non-linearity */
linearize_rspec(rspec * raw)467 void linearize_rspec(rspec *raw) {
468 rspec_inf *inf = raw->inf;
469 int i, j;
470 rspec *dark;
471
472 if (raw->state & rspec_lin)
473 error("linearize_rspec: already done");
474
475 if (raw->state & rspec_int)
476 error("linearize_rspec: can't be integration time adjusted");
477
478 if (!(raw->state & rspec_dcal))
479 error("linearize_rspec: needs black subtract");
480
481 if (inf->nlin > 0) {
482 for (i = 0; i < raw->nmeas; i++) {
483 for (j = 0; j < raw->nsamp; j++) {
484 raw->samp[i][j] = linearize_val_rspec(inf, raw->samp[i][j]);
485 }
486 }
487 }
488 raw->state |= rspec_lin;
489 }
490
491 /* Apply the emsissive calibration */
emis_calibrate_rspec(rspec * raw)492 void emis_calibrate_rspec(rspec *raw) {
493 rspec_inf *inf = raw->inf;
494 int i, j;
495
496 if (raw->state & rspec_cal)
497 error("emis_calibrate_rspec: already done");
498
499 if (raw->stype != raw->inf->ecaltype)
500 error("emis_calibrate_rspec: ecaltype does not match rspec type");
501
502 for (i = 0; i < raw->nmeas; i++) {
503 for (j = 0; j < raw->nsamp; j++) {
504 raw->samp[i][j] *= inf->ecal[j];
505 }
506 }
507 raw->state |= rspec_cal;
508 }
509
510 /* Scale to the integration time */
inttime_calibrate_rspec(rspec * raw)511 void inttime_calibrate_rspec(rspec *raw) {
512 rspec_inf *inf = raw->inf;
513 int i, j;
514
515 if (raw->state & rspec_int)
516 error("inttime_calibrate_rspec: already done");
517
518 for (i = 0; i < raw->nmeas; i++) {
519 for (j = 0; j < raw->nsamp; j++) {
520 raw->samp[i][j] /= raw->inttime;
521 }
522 }
523
524 raw->inttime = 1.0;
525
526 raw->state |= rspec_int;
527 }
528
529 /* return a wav rspec from a raw rspec */
530 /* (This does not make any adjustments to the values) */
convert_wav_from_raw_rspec(rspec * raw)531 rspec *convert_wav_from_raw_rspec(rspec *raw) {
532 rspec_inf *inf = raw->inf;
533 rspec *wav;
534 int cx, sx, i, j, k;
535
536 if (raw->stype != rspec_raw)
537 error("extract_raw_from_sensor_rspec: input is not raw type");
538
539 wav = new_rspec(raw->inf, rspec_wav, raw->nmeas);
540
541 wav->mtype = raw->mtype;
542 wav->state = raw->state;
543 wav->inttime = raw->inttime;
544
545 for (i = 0; i < wav->nmeas; i++) { /* For each measurement */
546 for (cx = j = 0; j < inf->nwav; j++) { /* For each wav sample */
547 double oval = 0.0;
548
549 sx = inf->findex[j]; /* Starting index */
550 for (k = 0; k < inf->fnocoef[j]; k++, cx++, sx++) /* For each matrix value */
551 oval += inf->fcoef[cx] * raw->samp[i][sx];
552 wav->samp[i][j] = oval;
553 }
554 }
555
556 return wav;
557 }
558
559 /* -------------------------------------------------- */
560
561
562 /* Filter code in i1pro_imp is in:
563
564 i1pro_compute_wav_filters() X-Rite way
565 i1pro_create_hr() Using gausian
566
567 */
568
569 /* Resampling kernels. (There are more in i1pro_imp.c) */
570
571 /* They aren't expected to be unity area, as they will be */
572 /* normalized anyway. */
573 /* wi is the width of the filter */
574
triangle(double wi,double x)575 static double triangle(double wi, double x) {
576 double y = 0.0;
577
578 x = fabs(x/wi);
579 y = 1.0 - x;
580 if (y < 0.0)
581 y = 0.0;
582
583 return y;
584 }
585
gausian(double wi,double x)586 static double gausian(double wi, double x) {
587 double y = 0.0;
588
589 wi = wi/(sqrt(2.0 * log(2.0))); /* Convert width at half max to std. dev. */
590 x = x/wi;
591 y = exp(-(x * x)); /* Center at 1.0 */
592
593 return y;
594 }
595
lanczos2(double wi,double x)596 static double lanczos2(double wi, double x) {
597 double y = 0.0;
598
599 wi *= 1.05; // Improves smoothness. Why ?
600
601 x = fabs(1.0 * x/wi);
602 if (x >= 2.0)
603 return 0.0;
604 if (x < 1e-6)
605 return 1.0;
606 y = sin(DBL_PI * x)/(DBL_PI * x) * sin(DBL_PI * x/2.0)/(DBL_PI * x/2.0);
607
608 return y;
609 }
610
lanczos3(double wi,double x)611 static double lanczos3(double wi, double x) {
612 double y = 0.0;
613
614 x = fabs(1.0 * x/wi);
615 if (x >= 3.0)
616 return 0.0;
617 if (x < 1e-6)
618 return 1.0;
619 y = sin(DBL_PI * x)/(DBL_PI * x) * sin(DBL_PI * x/3.0)/(DBL_PI * x/3.0);
620
621 return y;
622 }
623
cubicspline(double wi,double x)624 static double cubicspline(double wi, double x) {
625 double y = 0.0;
626 double xx = x;
627 double bb, cc;
628
629 xx = fabs(1.0 * x/wi);
630
631 // bb = cc = 1.0/3.0; /* Mitchell */
632 bb = 0.5;
633 cc = 0.5;
634
635 if (xx < 1.0) {
636 y = ( 12.0 - 9.0 * bb - 6.0 * cc) * xx * xx * xx
637 + (-18.0 + 12.0 * bb + 6.0 * cc) * xx * xx
638 + ( 6.0 - 2.0 * bb);
639 y /= (6.0 - 2.0 * bb);
640 } else if (xx < 2.0) {
641 y = ( -1.0 * bb - 6.0 * cc) * xx * xx * xx
642 + ( 6.0 * bb + 30.0 * cc) * xx * xx
643 + (-12.0 * bb - 48.0 * cc) * xx
644 + ( 8.0 * bb + 24.0 * cc);
645 y /= (6.0 - 2.0 * bb);
646 } else {
647 y = 0.0;
648 }
649
650 return y;
651 }
652
653 /* Create the wavelength resampling filters */
rspec_make_resample_filters(rspec_inf * inf)654 void rspec_make_resample_filters(rspec_inf *inf) {
655 double twidth = inf->wl_space;
656 double rawspace; /* Average raw band spacing wl */
657 double fshmax; /* filter shape max wavelength from center */
658 double finc; /* Integration step size */
659 int maxcoeffs; /* Maximum coefficients per filter */
660 int **coeff_ix; /* [band][coef] Raw index */
661 double **coeff_we; /* [band][coef] Weighting */
662 double (*kernel)(double wi, double x) = NULL; /* Filter kernel */
663 int xcount;
664 int i, j, k;
665
666 if (inf->ktype == rspec_triangle)
667 kernel = triangle;
668 else if (inf->ktype == rspec_gausian)
669 kernel = gausian;
670 else if (inf->ktype == rspec_lanczos2)
671 kernel = lanczos2;
672 else if (inf->ktype == rspec_lanczos3)
673 kernel = lanczos3;
674 else if (inf->ktype == rspec_cubicspline)
675 kernel = cubicspline;
676 else
677 error("rspec_make_resample_filters: unknown kernel %d",inf->ktype);
678
679 #ifdef NEVER // Check kernel sums to 1.0
680 {
681 double x, y;
682
683 for (x = 0.0; x < 5.0; x += 0.1) {
684 y = kernel(1.0, x - 4.0)
685 + kernel(1.0, x - 3.0)
686 + kernel(1.0, x - 2.0)
687 + kernel(1.0, x - 1.0)
688 + kernel(1.0, x)
689 + kernel(1.0, x + 1.0)
690 + kernel(1.0, x + 2.0);
691 + kernel(1.0, x + 3.0);
692 + kernel(1.0, x + 4.0);
693
694 printf("Offset %f sum %f\n",x,y);
695 }
696
697 }
698 #endif // NEVER
699
700 /* Aproximate raw value spacing in nm */
701 rawspace = (inf->wl_long - inf->wl_short)/inf->rawrange.num;
702 //printf("~1 rawspace = %f\n",rawspace);
703
704 /* Figure the extent of the filter kernel. We assume they */
705 /* all have a finite extent. */
706 for (fshmax = 50.0; fshmax >= 0.0; fshmax -= 0.01) {
707 if (fabs(kernel(twidth, fshmax)) > 1e-6) {
708 fshmax += 0.01;
709 break;
710 }
711 }
712 //printf("~1 fshmax = %f\n",fshmax);
713
714 if (fshmax <= 0.0)
715 error("rspec_make_resample_filters: fshmax search failed\n");
716
717 a1logd(inf->log, 4,"rspec_make_resample_filters: fshmax = %f\n",fshmax);
718
719 /* Figure number of raw samples over kernel extent. */
720 /* (Allow generous factor for non-linearity) */
721 maxcoeffs = (int)floor(2.0 * 1.4 * fshmax/rawspace + 3.0);
722
723 a1logd(inf->log, 4,"rspec_make_resample_filters: maxcoeffs = %d\n",maxcoeffs);
724
725 /* Figure out box integration step size */
726 #ifdef FAST_HIGH_RES_SETUP
727 finc = twidth/50.0;
728 if (rawspace/finc < 10.0)
729 finc = rawspace/10.0;
730 #else
731 finc = twidth/15.0;
732 if (rawspace/finc < 4.0)
733 finc = rawspace/4.0;
734 #endif
735
736 a1logd(inf->log, 4,"rspec_make_resample_filters: integration step = %f\n",finc);
737
738 if (inf->fnocoef != NULL)
739 free(inf->fnocoef);
740 if ((inf->fnocoef = (int *)calloc(inf->nwav, sizeof(int))) == NULL)
741 error("rspec_make_resample_filters: malloc failure");
742
743 /* Space to build filter coeficients */
744 coeff_ix = imatrix(0, inf->nwav-1, 0, maxcoeffs-1);
745 coeff_we = dmatrix(0, inf->nwav-1, 0, maxcoeffs-1);
746
747 /* For all the usable raw bands */
748 for (i = inf->rawrange.off+1; i < (inf->rawrange.off+inf->rawrange.num-1); i++) {
749 double w1, wl, w2;
750
751 /* Translate CCD center and boundaries to calibrated wavelength */
752 wl = rspec_raw2nm(inf, (double)i);
753 w1 = rspec_raw2nm(inf, (double)i - 0.5);
754 w2 = rspec_raw2nm(inf, (double)i + 0.5);
755
756 // printf("~1 CCD %d, w1 %f, wl %f, w2 %f\n",i,w1,wl,w2);
757
758 /* For each output filter */
759 for (j = 0; j < inf->nwav; j++) {
760 double cwl, rwl; /* center, relative wavelegth */
761 double we;
762
763 cwl = rspec_wav2nm(inf, (double)j);
764 rwl = wl - cwl; /* raw relative wavelength to filter */
765
766 if (fabs(w1 - cwl) > fshmax && fabs(w2 - cwl) > fshmax)
767 continue; /* Doesn't fall into this filter */
768
769 #ifdef BOX_INTEGRATE
770 /* Integrate in finc nm increments from filter shape */
771 /* using triangular integration. */
772 {
773 int nn;
774 double lw, ll;
775
776 nn = (int)(fabs(w2 - w1)/finc + 0.5); /* Number to integrate over */
777
778 lw = w1; /* start at lower boundary of CCD cell */
779 ll = kernel(twidth, w1 - cwl);
780 we = 0.0;
781 for (k = 0; k < nn; k++) {
782 double cw, cl;
783
784 #if defined(__APPLE__) && defined(__POWERPC__)
785 gcc_bug_fix(k);
786 #endif
787 cw = w1 + (k+1.0)/(nn + 1.0) * fabs(w2 - w1); /* wl to sample */
788 cl = kernel(twidth, cw - cwl);
789 we += 0.5 * (cl + ll) * fabs(lw - cw); /* Area under triangle */
790 ll = cl;
791 lw = cw;
792 }
793 }
794 #else
795 we = fabs(w2 - w1) * kernel(twidth, rwl);
796 #endif
797
798 if (inf->fnocoef[j] >= maxcoeffs)
799 error("rspec_make_resample_filters: run out of high res filter space\n");
800
801 coeff_ix[j][inf->fnocoef[j]] = i;
802 coeff_we[j][inf->fnocoef[j]++] = we;
803 // printf("~1 filter %d, cwl %f, rwl %f, ix %d, we %f, nocoefs %d\n",j,cwl,rwl,i,we,info->fnocoef[j]);
804 }
805 }
806
807 /* Convert hires filters into runtime format: */
808
809 /* Allocate or reallocate high res filter tables */
810 if (inf->findex != NULL)
811 free(inf->findex);
812 if (inf->fcoef != NULL)
813 free(inf->fcoef);
814
815 if ((inf->findex = (int *)calloc(inf->nraw, sizeof(int))) == NULL)
816 error("rspec_make_resample_filters: malloc index failed!\n");
817
818 /* Count the total number of coefficients */
819 for (xcount = j = 0; j < inf->nwav; j++) {
820 inf->findex[j] = coeff_ix[j][0]; /* raw starting index */
821 xcount += inf->fnocoef[j];
822 }
823
824 //printf("~1 total coefs = %d\n",xcount);
825
826 /* Allocate space for them */
827 if ((inf->fcoef = (double *)calloc(xcount, sizeof(double))) == NULL)
828 error("rspec_make_resample_filters: malloc index failed!\n");
829
830 /* Normalize the weight * nm to 1.0, and pack them into the run-time format */
831 for (i = j = 0; j < inf->nwav; j++) {
832 int sx;
833 double rwi, twe = 0.0;
834
835 sx = inf->findex[j]; /* raw starting index */
836
837 for (k = 0; k < inf->fnocoef[j]; sx++, k++) {
838 /* Width of raw band in nm */
839 rwi = fabs(rspec_raw2nm(inf, (double)sx - 0.5)
840 - rspec_raw2nm(inf, (double)sx + 0.5));
841 twe += rwi * coeff_we[j][k];
842 }
843
844 if (twe > 0.0)
845 twe = 1.0/twe;
846 else
847 twe = 1.0;
848
849 // printf("Output %d, nocoefs %d, norm weight %f:\n",j,inf->fnocoef[j],twe);
850
851 for (k = 0; k < inf->fnocoef[j]; k++, i++) {
852 inf->fcoef[i] = coeff_we[j][k] * twe;
853 // printf(" coef %d packed %d from raw %d val %f\n",k,i,inf->findex[j]+k,inf->fcoef[i]);
854 }
855 }
856
857 free_imatrix(coeff_ix, 0, inf->nwav-1, 0, maxcoeffs-1);
858 free_dmatrix(coeff_we, 0, inf->nwav-1, 0, maxcoeffs-1);
859 }
860
861 //printf("~1 line %d\n",__LINE__);
862
863 /* Plot the wave resampling filters */
plot_resample_filters(rspec_inf * inf)864 void plot_resample_filters(rspec_inf *inf) {
865 #ifndef SALONEINSTLIB
866 double *xx, *ss;
867 double **yy;
868 int i, j, k, sx;
869
870 //printf("~1 nraw = %d\n",inf->nraw);
871
872 xx = dvectorz(0, inf->nraw-1); /* X index */
873 yy = dmatrixz(0, 5, 0, inf->nraw-1); /* Curves distributed amongst 5 graphs */
874 /* with 6th holding sum */
875 for (i = 0; i < inf->nraw; i++)
876 xx[i] = i;
877
878 /* For each output wavelength */
879 for (i = j = 0; j < inf->nwav; j++) {
880
881 sx = inf->findex[j]; /* raw starting index */
882 //printf("Output %d raw sx %d\n",j,sx);
883
884 /* For each matrix value */
885 for (k = 0; k < inf->fnocoef[j]; k++, sx++, i++) {
886 yy[5][sx] += 0.5 * inf->fcoef[i];
887 yy[j % 5][sx] = inf->fcoef[i];
888 //printf(" filter %d six %d weight = %e\n",k,sx,inf->fcoef[i]);
889 }
890 }
891
892 printf("Wavelength re-sampling curves:\n");
893 // do_plot6(xx, yy[0], yy[1], yy[2], yy[3], yy[4], yy[5], 150);
894 do_plot6(xx, yy[0], yy[1], yy[2], yy[3], yy[4], yy[5], inf->nraw);
895 free_dvector(xx, 0, inf->nraw-1);
896 free_dmatrix(yy, 0, 2, 0, inf->nraw-1);
897 #endif
898 }
899
900 /* ================================================== */
901 /* Calibration file support */
902
903 /* Open the file. nz wr for write mode, else read */
904 /* Return nz on error */
calf_open(calf * x,a1log * log,char * fname,int wr)905 int calf_open(calf *x, a1log *log, char *fname, int wr) {
906 char nmode[10];
907 char cal_name[200];
908 char **cal_paths = NULL;
909 int no_paths = 0;
910
911 memset((void *)x, 0, sizeof(calf));
912 x->log = log;
913
914 if (wr)
915 strcpy(nmode, "w");
916 else
917 strcpy(nmode, "r");
918
919 #if !defined(O_CREAT) && !defined(_O_CREAT)
920 # error "Need to #include fcntl.h!"
921 #endif
922 #if defined(O_BINARY) || defined(_O_BINARY)
923 strcat(nmode, "b");
924 #endif
925
926 /* Create the file name */
927 if (wr)
928 sprintf(cal_name, "ArgyllCMS/%s", fname);
929 else
930 sprintf(cal_name, "ArgyllCMS/%s" SSEPS "color/%s", fname, fname);
931 if ((no_paths = xdg_bds(NULL, &cal_paths, xdg_cache, xdg_write, xdg_user, cal_name)) < 1) {
932 a1logd(x->log,1,"calf_open: xdg_bds returned no paths\n");
933 return 1;
934 }
935
936 a1logd(x->log,2,"calf_open: %s file '%s'\n",cal_paths[0], wr ? "saving to" : "restoring from");
937
938 /* Check the last modification time */
939 if (!wr) {
940 struct sys_stat sbuf;
941
942 if (sys_stat(cal_paths[0], &sbuf) == 0) {
943 x->lo_secs = time(NULL) - sbuf.st_mtime;
944 a1logd(x->log,2,"calf_open:: %d secs from instrument last open\n",x->lo_secs);
945 } else {
946 a1logd(x->log,2,"calf_open:: stat on file failed\n");
947 }
948 }
949
950 if ((wr && create_parent_directories(cal_paths[0]))
951 || (x->fp = fopen(cal_paths[0], nmode)) == NULL) {
952 a1logd(x->log,2,"calf_open: failed to open file for %s\n",wr ? "writing" : "reading");
953 xdg_free(cal_paths, no_paths);
954 return 1;
955 }
956 xdg_free(cal_paths, no_paths);
957
958 a1logd(x->log,2,"calf_open: suceeded\n");
959
960 return 0;
961 }
962
963 /* Update the modification time */
964 /* Return nz on error */
calf_touch(a1log * log,char * fname)965 int calf_touch(a1log *log, char *fname) {
966 char cal_name[200];
967 char **cal_paths = NULL;
968 int no_paths = 0;
969 int rv;
970
971 /* Locate the file name */
972 sprintf(cal_name, "ArgyllCMS/%s" SSEPS "color/%s", fname, fname);
973
974 if ((no_paths = xdg_bds(NULL, &cal_paths, xdg_cache, xdg_read, xdg_user, cal_name)) < 1) {
975 a1logd(log,2,"calf_touch: xdg_bds failed to locate file'\n");
976 return 1;
977 }
978
979 a1logd(log,2,"calf_touch: touching file '%s'\n",cal_paths[0]);
980
981 if ((rv = sys_utime(cal_paths[0], NULL)) != 0) {
982 a1logd(log,2,"calf_touch: failed with %d\n",rv);
983 xdg_free(cal_paths, no_paths);
984 return 1;
985 }
986 xdg_free(cal_paths, no_paths);
987
988 return 0;
989 }
990
991 /* Rewind and reset for another read */
calf_rewind(calf * x)992 void calf_rewind(calf *x) {
993 x->ef = 0;
994 x->chsum = 0;
995 x->nbytes = 0;
996 rewind(x->fp);
997 }
998
999 /* Close the file and free any memory */
1000 /* return nz on error */
calf_done(calf * x)1001 int calf_done(calf *x) {
1002 int rv = 0;
1003
1004 if (x->fp != NULL) {
1005 if (fclose(x->fp)) {
1006 a1logd(x->log,2,"calf_done: closing file failed\n");
1007 rv = 1;
1008 }
1009 }
1010
1011 if (x->buf != NULL)
1012 free(x->buf);
1013 x->buf = NULL;
1014 return rv;
1015 }
1016
sizebuf(calf * x,size_t size)1017 static void sizebuf(calf *x, size_t size) {
1018 if (x->bufsz < size)
1019 x->buf = realloc(x->buf, size);
1020 if (x->buf == NULL)
1021 error("calf: sizebuf malloc failed");
1022 }
1023
update_chsum(calf * x,unsigned char * p,int nn)1024 static void update_chsum(calf *x, unsigned char *p, int nn) {
1025 int i;
1026 for (i = 0; i < nn; i++, p++)
1027 x->chsum = ((x->chsum << 13) | (x->chsum >> (32-13))) + *p;
1028 x->nbytes += nn;
1029 }
1030
1031 /* Write an array of ints to the file. Set the error flag to nz on error */
calf_wints(calf * x,int * dp,int n)1032 void calf_wints(calf *x, int *dp, int n) {
1033 if (x->ef)
1034 return;
1035
1036 if (fwrite((void *)dp, sizeof(int), n, x->fp) != n) {
1037 x->ef = 1;
1038 a1logd(x->log,2,"calf_wints: write failed for %d ints at offset %d\n",n,x->nbytes);
1039 } else {
1040 update_chsum(x, (unsigned char *)dp, sizeof(int) * n);
1041 }
1042 }
1043
1044 /* Write an array of doubles to the file. Set the error flag to nz on error */
calf_wdoubles(calf * x,double * dp,int n)1045 void calf_wdoubles(calf *x, double *dp, int n) {
1046 if (x->ef)
1047 return;
1048
1049 if (fwrite((void *)dp, sizeof(double), n, x->fp) != n) {
1050 x->ef = 1;
1051 a1logd(x->log,2,"calf_wdoubles: write failed for %d doubles at offset %d\n",n,x->nbytes);
1052 } else {
1053 update_chsum(x, (unsigned char *)dp, sizeof(double) * n);
1054 }
1055 }
1056
1057 /* Write an array of time_t's to the file. Set the error flag to nz on error */
1058 /* (This will cause file checksum fail if different executables on the same */
1059 /* system have different time_t values) */
calf_wtime_ts(calf * x,time_t * dp,int n)1060 void calf_wtime_ts(calf *x, time_t *dp, int n) {
1061 if (x->ef)
1062 return;
1063
1064 if (fwrite((void *)dp, sizeof(time_t), n, x->fp) != n) {
1065 x->ef = 1;
1066 a1logd(x->log,2,"calf_wtime_ts: write failed for %d time_ts at offset %d\n",n,x->nbytes);
1067 } else {
1068 update_chsum(x, (unsigned char *)dp, sizeof(time_t) * n);
1069 }
1070 }
1071
1072 /* Write a zero terminated string */
calf_wstrz(calf * x,char * dp)1073 void calf_wstrz(calf *x, char *dp) {
1074 int n;
1075
1076 if (x->ef)
1077 return;
1078
1079 n = strlen(dp) + 1;
1080
1081 calf_wints(x, &n, 1);
1082
1083 if (fwrite((void *)dp, sizeof(char), n, x->fp) != n) {
1084 x->ef = 1;
1085 a1logd(x->log,2,"calf_wstrz: write failed for %d long string at offset %d\n",n,x->nbytes);
1086 } else {
1087 update_chsum(x, (unsigned char *)dp, sizeof(char) * n);
1088 }
1089 }
1090
1091
1092 /* Read an array of ints from the file. Set the error flag to nz on error */
1093 /* Always read (ignore rd flag) */
calf_rints2(calf * x,int * dp,int n)1094 void calf_rints2(calf *x, int *dp, int n) {
1095 if (x->ef)
1096 return;
1097
1098 if (fread((void *)dp, sizeof(int), n, x->fp) != n) {
1099 x->ef = 1;
1100 a1logd(x->log,2,"calf_rints2: read failed for %d ints at offset %d\n",n,x->nbytes);
1101 } else {
1102 update_chsum(x, (unsigned char *)dp, sizeof(int) * n);
1103 }
1104 }
1105
1106
1107 /* Read an array of ints from the file. Set the error flag to nz on error */
calf_rints(calf * x,int * dp,int n)1108 void calf_rints(calf *x, int *dp, int n) {
1109 size_t nbytes = n * sizeof(int);
1110 void *dest = (void *)dp;
1111
1112 if (x->ef)
1113 return;
1114
1115 if (x->rd == 0) { /* Dummy read */
1116 sizebuf(x, nbytes);
1117 dest = x->buf;
1118 }
1119
1120 if (fread(dest, 1, nbytes, x->fp) != nbytes) {
1121 x->ef = 1;
1122 a1logd(x->log,2,"calf_rints: read failed for %d ints at offset %d\n",n,x->nbytes);
1123 } else {
1124 update_chsum(x, dest, nbytes);
1125 }
1126 }
1127
1128 /* Read an array of doubles from the file. Set the error flag to nz on error */
calf_rdoubles(calf * x,double * dp,int n)1129 void calf_rdoubles(calf *x, double *dp, int n) {
1130 size_t nbytes = n * sizeof(double);
1131 void *dest = (void *)dp;
1132
1133 if (x->ef)
1134 return;
1135
1136 if (x->rd == 0) { /* Dummy read */
1137 sizebuf(x, nbytes);
1138 dest = x->buf;
1139 }
1140
1141 if (fread(dest, 1, nbytes, x->fp) != nbytes) {
1142 x->ef = 1;
1143 a1logd(x->log,2,"calf_rdoubles: read failed for %d ints at offset %d\n",n,x->nbytes);
1144 } else {
1145 update_chsum(x, dest, nbytes);
1146 }
1147 }
1148
1149 /* Read an array of time_t's from the file. Set the error flag to nz on error */
1150 /* (This will cause file checksum fail if different executables on the same */
1151 /* system have different time_t values) */
calf_rtime_ts(calf * x,time_t * dp,int n)1152 void calf_rtime_ts(calf *x, time_t *dp, int n) {
1153 size_t nbytes = n * sizeof(time_t);
1154 void *dest = (void *)dp;
1155
1156 if (x->ef)
1157 return;
1158
1159 if (x->rd == 0) { /* Dummy read */
1160 sizebuf(x, nbytes);
1161 dest = x->buf;
1162 }
1163
1164 if (fread(dest, 1, nbytes, x->fp) != nbytes) {
1165 x->ef = 1;
1166 a1logd(x->log,2,"calf_rtime_ts: read failed for %d ints at offset %d\n",n,x->nbytes);
1167 } else {
1168 update_chsum(x, dest, nbytes);
1169 }
1170 }
1171
1172 /* Read a zero terminated string. */
calf_rstrz(calf * x,char ** dp)1173 void calf_rstrz(calf *x, char **dp) {
1174 int n;
1175 size_t nbytes = 0;
1176 char *dest = NULL;
1177
1178 if (x->ef)
1179 return;
1180
1181 calf_rints2(x, &n, 1);
1182 nbytes = sizeof(char) * n;
1183
1184 if (x->ef || n == 0)
1185 return;
1186
1187 if (x->rd != 0) { /* Reading for real */
1188 if (*dp != NULL)
1189 free(*dp);
1190 if ((*dp = dest = malloc(nbytes)) == NULL)
1191 error("calf: calf_rstrz malloc failed");
1192 } else {
1193 sizebuf(x, nbytes);
1194 dest = x->buf;
1195 }
1196 if (fread(dest, 1, nbytes, x->fp) != nbytes) {
1197 x->ef = 1;
1198 a1logd(x->log,2,"calf_rstrz: read failed for %d long string at offset %d\n",n,x->nbytes);
1199 } else {
1200 update_chsum(x, (unsigned char*)dest, nbytes);
1201 }
1202 }
1203
calf_rstrz2(calf * x,char ** dp)1204 void calf_rstrz2(calf *x, char **dp) {
1205 int rd = x->rd;
1206 x->rd = 1;
1207 calf_rstrz(x, dp);
1208 x->rd = rd;
1209 }
1210
1211 /* ================================================== */
1212
1213 /* Save a rspec to a calibration file */
calf_wrspec(calf * x,rspec * s)1214 void calf_wrspec(calf *x, rspec *s) {
1215 int i;
1216
1217 if (x->ef)
1218 return;
1219
1220 calf_wints(x, (int *)&s->stype, 1);
1221 calf_wints(x, (int *)&s->mtype, 1);
1222 calf_wints(x, (int *)&s->state, 1);
1223 calf_wdoubles(x, &s->inttime, 1);
1224 calf_wints(x, &s->nmeas, 1);
1225 calf_wints(x, &s->nsamp, 1);
1226
1227 for (i = 0; i < s->nmeas; i++) {
1228 calf_wdoubles(x, s->samp[i], s->nsamp);
1229 }
1230 }
1231
1232 /* Create a rspec from a calibration file */
calf_rrspec(calf * x,rspec ** dp,rspec_inf * inf)1233 void calf_rrspec(calf *x, rspec **dp, rspec_inf *inf) {
1234 rspec *s, dumy;
1235 int no, i;
1236
1237 if (x->ef)
1238 return;
1239
1240 if (x->rd != 0) {
1241 if (*dp != NULL)
1242 del_rspec(*dp);
1243 *dp = s = new_rspec(inf, rspec_sensor, 0);
1244 } else {
1245 s = &dumy;
1246 }
1247
1248 calf_rints2(x, (int *)&s->stype, 1);
1249 calf_rints2(x, (int *)&s->mtype, 1);
1250 calf_rints2(x, (int *)&s->state, 1);
1251 calf_rdoubles(x, &s->inttime, 1);
1252 calf_rints2(x, &s->nmeas, 1);
1253 calf_rints2(x, &s->nsamp, 1);
1254
1255 /* Sanity check. */
1256 no = rspec_typesize(inf, s->stype);
1257 if (no != s->nsamp) {
1258 a1logd(inf->log, 4,"calf_rrspec: unexpected nsamp %d (expect %d)\n",s->nsamp,no);
1259 x->ef = 1;
1260 return;
1261 }
1262
1263 if (x->rd != 0) {
1264 s->samp = dmatrix(0, s->nmeas-1, 0, s->nsamp-1);
1265 for (i = 0; i < s->nmeas; i++) {
1266 calf_rdoubles(x, s->samp[i], s->nsamp);
1267 }
1268 } else {
1269 for (i = 0; i < s->nmeas; i++) {
1270 calf_rdoubles(x, NULL, s->nsamp);
1271 }
1272 }
1273 }
1274
1275