summaryrefslogtreecommitdiff
path: root/spectro/munki_imp.c
diff options
context:
space:
mode:
Diffstat (limited to 'spectro/munki_imp.c')
-rw-r--r--spectro/munki_imp.c1263
1 files changed, 906 insertions, 357 deletions
diff --git a/spectro/munki_imp.c b/spectro/munki_imp.c
index c76e9bc..9f431cc 100644
--- a/spectro/munki_imp.c
+++ b/spectro/munki_imp.c
@@ -5,7 +5,7 @@
* Author: Graeme W. Gill
* Date: 12/1/2009
*
- * Copyright 2006 - 2013, Graeme W. Gill
+ * Copyright 2006 - 2014, Graeme W. Gill
* All rights reserved.
*
* This material is licenced under the GNU GENERAL PUBLIC LICENSE Version 2 or later :-
@@ -67,14 +67,22 @@
#include "sort.h"
/* Configuration */
-#undef USE_HIGH_GAIN_MODE /* [Und] Make use of high gain mode */
+#undef USE_HIGH_GAIN_MODE /* [Und] Make use of high gain mode in emissive measurements */
#define USE_THREAD /* [Def] Need to use thread, or there are 1.5 second internal */
/* instrument delays ! */
-#define ENABLE_NONVCAL /* [Def] Enable saving calibration state between program runs in a file */
+#define ENABLE_NONVCAL /* [Def] Enable saving calibration state between program runs to a file */
#define ENABLE_NONLINCOR /* [Def] Enable non-linear correction */
/* NOTE :- high gain scaling will be stuffed if disabled! */
#define ENABLE_LEDTEMPC /* [Def] Enable LED temperature compensation */
-#define ENABLE_SPOS_CHECK /* [Def] Chech the sensor position is reasonable for measurement */
+#define ENABLE_BKDRIFTC /* [Def] Enable Emis. Black drift compensation using sheilded cell values */
+#define HEURISTIC_BKDRIFTC /* [Def] Use heusristic black drift correction version */
+#undef ENABLE_REFSTRAYC /* [Und] Enable Reflective stray light compensation */
+#define REFSTRAYC_FACTOR 0.000086 /* [0.00043] Compensation factor */
+#undef ENABLE_REFLEDINTER /* [Und] Enable Reflective LED interference correction */
+
+#define ENABLE_SPOS_CHECK /* [Def] Check the sensor position is reasonable for measurement */
+#define FILTER_SPOS_EVENTS /* [Def] Use a thread to filter SPOS event changes */
+#define FILTER_TIME 500 /* [500] Filter time in msec */
#define DCALTOUT (1 * 60 * 60) /* [1 Hrs] Dark Calibration timeout in seconds */
#define WCALTOUT (24 * 60 * 60) /* [24 Hrs] White Calibration timeout in seconds */
#define MAXSCANTIME 20.0 /* [20 Sec] Maximum scan time in seconds */
@@ -83,24 +91,27 @@
#define SINGLE_READ /* [Def] Use a single USB read for scan to eliminate latency issues. */
#define HIGH_RES /* [Def] Enable high resolution spectral mode code. Disable */
/* to break dependency on rspl library. */
+# undef FAST_HIGH_RES_SETUP /* Slightly better accuracy ? */
/* Debug [Und] */
#undef DEBUG /* Turn on extra messages & plots */
#undef PLOT_DEBUG /* Use plot to show readings & processing */
#undef PLOT_REFRESH /* Plot refresh rate measurement info */
+#undef PLOT_UPDELAY /* Plot data used to determine display update delay */
#undef RAWR_DEBUG /* Print out raw reading processing values */
#undef DUMP_SCANV /* Dump scan readings to a file "mkdump.txt" */
#undef DUMP_DARKM /* Append raw dark readings to file "mkddump.txt" */
+#undef DUMP_BKLED /* Save REFSTRAYC & REFLEDNOISE comp plot to "refbk1.txt" & "refbk2.txt" */
#undef APPEND_MEAN_EMMIS_VAL /* Append averaged uncalibrated reading to file "mkdump.txt" */
#undef IGNORE_WHITE_INCONS /* Ignore define reference reading inconsistency */
#undef TEST_DARK_INTERP /* Test out the dark interpolation (need DEBUG for plot) */
#undef PLOT_RCALCURVE /* Plot the reflection reference curve */
#undef PLOT_ECALCURVES /* Plot the emission reference curves */
-#undef PLOT_TEMPCOMP /* Plot before and after temp. compensation */
-#undef PATREC_DEBUG /* Print & Plot patch/flash recognition information */
+#undef PLOT_TEMPCOMP /* Plot before and after LED temp. compensation */
+#undef PLOT_PATREC /* Print & Plot patch/flash recognition information */
#undef HIGH_RES_DEBUG
#undef HIGH_RES_PLOT
-#undef HIGH_RES_PLOT_STRAYL /* Plot strat light upsample */
+#undef HIGH_RES_PLOT_STRAYL /* Plot stray light upsample */
#define DISP_INTT 0.7 /* Seconds per reading in display spot mode */
@@ -116,24 +127,22 @@
#define ADARKINT_MAX 2.0 /* Max cal time for adaptive dark cal with high gain mode */
#define ADARKINT_MAX2 4.0 /* Max cal time for adaptive dark for no high gain */
-#define SCAN_OP_LEV 0.10 /* Degree of optimimum sensor value to aim for */
- /* Make it scan as fast as possible */
-
#define RDEAD_TIME 0.004432 /* Fudge figure to make reflecting intn. time scale linearly */
#define EMIS_SCALE_FACTOR 1.0 /* Emission mode scale factor */
-#define AMB_SCALE_FACTOR (1.0/3.141592654) /* Ambient mode scale factor - convert */
- /* from Lux to Lux/PI */
- /* These factors get the same behaviour as the GMB drivers. */
+#define AMB_SCALE_FACTOR 1.0 /* Ambient mode scale factor for Lux */
#define NSEN_MAX 140 /* Maximum nsen/raw value we can cope with */
+/* Wavelength to start duplicating values below, because it is too noisy */
+#define WL_REF_MIN 420.0
+#define WL_EMIS_MIN 400.0
+
/* High res mode settings */
-#define HIGHRES_SHORT 360 /* Wavelength to calculate */
-#define HIGHRES_LONG 740
+#define HIGHRES_SHORT 380 /* Wavelength to calculate. Too noisy to try expanding range. */
+#define HIGHRES_LONG 730
#define HIGHRES_WIDTH (10.0/3.0) /* (The 3.3333 spacing and lanczos2 seems a good combination) */
-#define HIGHRES_REF_MIN 410.0 /* Too much stray light below this in reflective mode */
-#define HIGHRES_TRANS_MIN 380.0 /* Too much stray light below this in reflective mode */
+
#include "munki.h"
#include "munki_imp.h"
@@ -177,7 +186,7 @@ static void dump_bytes(a1log *log, char *pfx, unsigned char *buf, int base, int
bp += sprintf(bp,".");
}
bp += sprintf(bp,"\n");
- a1logd(log,0,oline);
+ a1logd(log,0,"%s",oline);
bp = oline;
}
}
@@ -186,7 +195,7 @@ static void dump_bytes(a1log *log, char *pfx, unsigned char *buf, int base, int
/* ============================================================ */
/* Debugging plot support */
-#if defined(DEBUG) || defined(PLOT_DEBUG) || defined(PATREC_DEBUG) || defined(HIGH_RES_PLOT) || defined(HIGH_RES_PLOT_STRAYL)
+#if defined(DEBUG) || defined(PLOT_DEBUG) || defined(PLOT_PATREC) || defined(HIGH_RES_PLOT) || defined(HIGH_RES_PLOT_STRAYL)
# include <plot.h>
@@ -322,6 +331,11 @@ void del_munkiimp(munki *p) {
munkiimp *m = (munkiimp *)p->m;
munki_state *s;
+#ifdef FILTER_SPOS_EVENTS
+ if (m->spos_th != NULL)
+ m->spos_th_term = 1;
+#endif
+
if (m->th != NULL) { /* Terminate switch monitor thread by simulating an event */
m->th_term = 1; /* Tell thread to exit on error */
munki_simulate_event(p, mk_eve_spos_change, 0);
@@ -331,9 +345,20 @@ void del_munkiimp(munki *p) {
a1logd(p->log,3,"Munki switch thread termination failed\n");
}
m->th->del(m->th);
- usb_uninit_cancel(&m->cancelt); /* Don't need cancel token now */
+ usb_uninit_cancel(&m->sw_cancel); /* Don't need cancel token now */
}
+#ifdef FILTER_SPOS_EVENTS
+ if (m->spos_th != NULL) {
+ for (i = 0; m->spos_th_termed == 0 && i < 5; i++)
+ msec_sleep(50); /* Wait for thread to terminate */
+ if (i >= 5) {
+ a1logd(p->log,3,"Munki spos thread termination failed\n");
+ }
+ m->spos_th->del(m->spos_th);
+ }
+#endif
+
/* Free any per mode data */
for (i = 0; i < mk_no_modes; i++) {
s = &m->ms[i];
@@ -568,11 +593,17 @@ munki_code munki_imp_init(munki *p) {
#ifdef USE_THREAD
/* Setup the switch monitoring thread */
- usb_init_cancel(&m->cancelt); /* Get cancel token ready */
+ usb_init_cancel(&m->sw_cancel); /* Get cancel token ready */
if ((m->th = new_athread(munki_switch_thread, (void *)p)) == NULL)
return MUNKI_INT_THREADFAILED;
#endif
+#ifdef FILTER_SPOS_EVENTS
+ /* Setup the sensor position filter thread */
+ if ((m->spos_th = new_athread(munki_spos_thread, (void *)p)) == NULL)
+ return MUNKI_INT_THREADFAILED;
+#endif
+
/* Set up the current state of each mode */
{
int i, j;
@@ -589,6 +620,11 @@ munki_code munki_imp_init(munki *p) {
s->targmaxitime = 2.0; /* Maximum integration time to aim for */
s->targoscale2 = 0.15; /* Proportion of targoscale to meed targmaxitime */
+#ifdef USE_HIGH_GAIN_MODE
+ s->auto_gain = 1; /* No high gain by default */
+#else
+ s->auto_gain = 0; /* No high gain by default */
+#endif
s->gainmode = 0; /* Normal gain mode */
s->inttime = 0.5; /* Initial integration time */
@@ -608,18 +644,14 @@ munki_code munki_imp_init(munki *p) {
s->idark_valid = 0; /* Interpolatable Dark cal invalid */
s->idark_data = dmatrixz(0, 3, -1, m->nraw-1);
- s->min_wl = 0.0; /* Default minimum to report */
-
s->dark_int_time = DISP_INTT; /* 0.7 */
s->dark_int_time2 = DISP_INTT2; /* 0.3 */
s->dark_int_time3 = DISP_INTT3; /* 0.1 */
s->idark_int_time[0] = s->idark_int_time[2] = m->min_int_time;
-#ifdef USE_HIGH_GAIN_MODE
- s->idark_int_time[1] = s->idark_int_time[3] = ADARKINT_MAX; /* 2.0 */
-#else
- s->idark_int_time[1] = s->idark_int_time[3] = ADARKINT_MAX2; /* 4.0 */
-#endif
+ s->idark_int_time[1] = ADARKINT_MAX; /* 2.0 */
+ s->idark_int_time[3] = ADARKINT_MAX2; /* 4.0 */
+
s->want_calib = 1; /* By default want an initial calibration */
s->want_dcalib = 1;
}
@@ -629,7 +661,8 @@ munki_code munki_imp_init(munki *p) {
s = &m->ms[i];
switch(i) {
case mk_refl_spot:
- s->targoscale = 1.0; /* Optimised sensor scaling to full */
+ s->auto_gain = 0; /* Don't automatically set gain */
+ s->targoscale = 1.0; /* Optimised sensor scaling to full */
s->reflective = 1;
s->adaptive = 0;
s->inttime = s->targoscale * m->cal_int_time;
@@ -642,17 +675,19 @@ munki_code munki_imp_init(munki *p) {
s->dreadtime = 0.5; /* same as reading */
s->wreadtime = 0.5;
s->maxscantime = 0.0;
- s->min_wl = HIGHRES_REF_MIN; /* Not enogh illumination to go below this */
break;
case mk_refl_scan:
- s->targoscale = SCAN_OP_LEV; /* Maximize scan rate */
+ s->auto_gain = 0; /* Don't automatically set gain */
+ s->targoscale = 1.0; /* Maximize level */
s->reflective = 1;
s->scan = 1;
s->adaptive = 0;
- s->inttime = (s->targoscale * m->cal_int_time - RDEAD_TIME) + RDEAD_TIME;
- if (s->inttime < m->min_int_time)
- s->inttime = m->min_int_time;
+// s->inttime = (s->targoscale * m->cal_int_time - RDEAD_TIME) + RDEAD_TIME;
+// if (s->inttime < m->min_int_time)
+// s->inttime = m->min_int_time;
+// s->dark_int_time = s->inttime;
+ s->inttime = s->targoscale * m->cal_int_time;
s->dark_int_time = s->inttime;
s->dpretime = 0.20; /* Pre-measure time */
@@ -662,7 +697,6 @@ munki_code munki_imp_init(munki *p) {
s->dreadtime = 0.10;
s->wreadtime = 0.10;
s->maxscantime = MAXSCANTIME;
- s->min_wl = HIGHRES_REF_MIN; /* Not enogh illumination to go below this */
break;
case mk_emiss_spot_na: /* Emissive spot not adaptive */
@@ -767,14 +801,13 @@ munki_code munki_imp_init(munki *p) {
s->dreadtime = 0.0;
s->wreadtime = 1.0;
s->maxscantime = 0.0;
- s->min_wl = HIGHRES_TRANS_MIN; /* Too much stray light below this ? */
break;
case mk_trans_scan:
- s->targoscale = 0.90; /* Allow extra 10% margine */
+ // ~~99 should we use high gain mode ??
+ s->targoscale = 0.10; /* Scan as fast as possible */
s->trans = 1;
s->scan = 1;
- s->targoscale = SCAN_OP_LEV;
s->inttime = s->targoscale * m->cal_int_time;
if (s->inttime < m->min_int_time)
s->inttime = m->min_int_time;
@@ -788,7 +821,6 @@ munki_code munki_imp_init(munki *p) {
s->dreadtime = 0.00;
s->wreadtime = 0.10;
s->maxscantime = MAXSCANTIME;
- s->min_wl = HIGHRES_TRANS_MIN; /* Too much stray light below this ? */
break;
}
}
@@ -837,12 +869,12 @@ char *munki_imp_get_serial_no(munki *p) {
/* Set the measurement mode. It may need calibrating */
munki_code munki_imp_set_mode(
munki *p,
- mk_mode mmode,
- int spec_en /* nz to enable reporting spectral */
+ mk_mode mmode, /* Operating mode */
+ inst_mode mode /* Full mode mask for options */
) {
munkiimp *m = (munkiimp *)p->m;
- a1logd(p->log,3,"munki_imp_set_mode called with %d\n",mmode);
+ a1logd(p->log,2,"munki_imp_set_mode called with mode no. %d and mask 0x%x\n",mmode,m);
switch(mmode) {
case mk_refl_spot:
case mk_refl_scan:
@@ -856,12 +888,21 @@ munki_code munki_imp_set_mode(
case mk_trans_spot:
case mk_trans_scan:
m->mmode = mmode;
- m->spec_en = spec_en ? 1 : 0;
- return MUNKI_OK;
- default:
break;
+ default:
+ return MUNKI_INT_ILLEGALMODE;
}
- return MUNKI_INT_ILLEGALMODE;
+ m->spec_en = (mode & inst_mode_spectral) != 0;
+
+ if ((mode & inst_mode_highres) != 0) {
+ munki_code rv;
+ if ((rv = munki_set_highres(p)) != MUNKI_OK)
+ return rv;
+ } else {
+ munki_set_stdres(p); /* Ignore any error */
+ }
+
+ return MUNKI_OK;
}
/* Return needed and available inst_cal_type's */
@@ -875,7 +916,7 @@ munki_code munki_imp_get_n_a_cals(munki *p, inst_cal_type *pn_cals, inst_cal_typ
a1logd(p->log,3,"munki_imp_get_n_a_cals: checking mode %d\n",m->mmode);
/* Timout calibrations that are too old */
- a1logd(p->log,4,"curtime = %u, iddate = %u\n",curtime,cs->iddate);
+ a1logd(p->log,4,"curtime %u, iddate %u, ddate %u, cfdate %u\n",curtime,cs->iddate,cs->ddate,cs->cfdate);
if ((curtime - cs->iddate) > DCALTOUT) {
a1logd(p->log,3,"Invalidating adaptive dark cal as %d secs from last cal\n",curtime - cs->iddate);
cs->idark_valid = 0;
@@ -972,7 +1013,7 @@ munki_code munki_imp_calibrate(
else if (*calt == inst_calt_available)
*calt = available & inst_calt_n_dfrble_mask;
- a1logd(p->log,4,"munki_imp_calibrate: doing calt 0x%x\n",calt);
+ a1logd(p->log,4,"munki_imp_calibrate: doing calt 0x%x\n",*calt);
if ((*calt & inst_calt_n_dfrble_mask) == 0) /* Nothing todo */
return MUNKI_OK;
@@ -989,6 +1030,16 @@ munki_code munki_imp_calibrate(
}
a1logd(p->log,4,"munki sensor position = 0x%x\n",spos);
+#ifdef NEVER
+ /* We can set the *calc to the actual conditions, in which case */
+ /* the calibration will commence immediately. */
+ if (!m->nosposcheck && spos == mk_spos_calib) {
+ *calc = inst_calc_man_cal_smode;
+ a1logd(p->log,4,"munki set calc to cal conditions\n",spos);
+
+ }
+#endif
+
/* Make sure that the instrument configuration matches the */
/* conditions */
if (*calc == inst_calc_man_cal_smode) {
@@ -1001,13 +1052,6 @@ munki_code munki_imp_calibrate(
}
}
- /* If the instrument is in the calibration position, */
- /* we know what the conditions are. */
- if (!m->nosposcheck && spos == mk_spos_calib) {
- *calc = inst_calc_man_cal_smode;
- a1logd(p->log,4,"munki set calc to cal conditions\n",spos);
- }
-
a1logd(p->log,4,"munki_imp_calibrate has right conditions\n");
if (*calt & inst_calt_ap_flag) {
@@ -1025,9 +1069,9 @@ munki_code munki_imp_calibrate(
/* Sanity check scan mode settings, in case something strange */
/* has been restored from the persistence file. */
- if (s->scan && s->inttime > (2.1 * m->min_int_time)) {
- s->inttime = m->min_int_time; /* Maximize scan rate */
- }
+// if (s->scan && s->inttime > (2.1 * m->min_int_time)) {
+// s->inttime = m->min_int_time; /* Maximize scan rate */
+// }
/* We are now either in inst_calc_man_cal_smode, */
/* inst_calc_man_trans_white, inst_calc_disp_white or inst_calc_proj_white */
@@ -1053,7 +1097,7 @@ munki_code munki_imp_calibrate(
nummeas = munki_comp_nummeas(p, s->dcaltime, s->inttime);
- a1logd(p->log,3,"\nDoing initial display black calibration with dcaltime %f, int_time %f, nummeas %d, gainmode %d\n", s->dcaltime, s->inttime, nummeas, s->gainmode);
+ a1logd(p->log,3,"\nDoing initial black calibration with dcaltime %f, int_time %f, nummeas %d, gainmode %d\n", s->dcaltime, s->inttime, nummeas, s->gainmode);
stm = msec_time();
if ((ev = munki_dark_measure(p, s->dark_data, nummeas, &s->inttime, s->gainmode))
!= MUNKI_OK) {
@@ -1207,25 +1251,25 @@ if (ss->dark_int_time2 != s->dark_int_time2
return ev;
}
-#ifdef USE_HIGH_GAIN_MODE
- s->idark_int_time[2] = m->min_int_time; /* 0.01 */
- nummeas = munki_comp_nummeas(p, s->dcaltime, s->idark_int_time[2]);
- a1logd(p->log,3,"Doing adaptive interpolated black calibration, dcaltime %f, idark_int_time[2] %f, nummeas %d, gainmode %d\n", s->dcaltime, s->idark_int_time[2], nummeas, 1);
- if ((ev = munki_dark_measure(p, s->idark_data[2], nummeas, &s->idark_int_time[2], 1))
- != MUNKI_OK) {
- m->mmode = mmode; /* Restore actual mode */
- return ev;
- }
-
- s->idark_int_time[3] = ADARKINT_MAX; /* 2.0 */
- a1logd(p->log,3,"Doing adaptive interpolated black calibration, dcaltime %f, idark_int_time[3] %f, nummeas %d, gainmode %d\n", s->dcaltime, s->idark_int_time[3], nummeas, 1);
- nummeas = munki_comp_nummeas(p, s->dcaltime, s->idark_int_time[3]);
- if ((ev = munki_dark_measure(p, s->idark_data[3], nummeas, &s->idark_int_time[3], 1))
- != MUNKI_OK) {
- m->mmode = mmode; /* Restore actual mode */
- return ev;
+ if (s->auto_gain) { /* If high gain is permitted */
+ s->idark_int_time[2] = m->min_int_time; /* 0.01 */
+ nummeas = munki_comp_nummeas(p, s->dcaltime, s->idark_int_time[2]);
+ a1logd(p->log,3,"Doing adaptive interpolated black calibration, dcaltime %f, idark_int_time[2] %f, nummeas %d, gainmode %d\n", s->dcaltime, s->idark_int_time[2], nummeas, 1);
+ if ((ev = munki_dark_measure(p, s->idark_data[2], nummeas, &s->idark_int_time[2], 1))
+ != MUNKI_OK) {
+ m->mmode = mmode; /* Restore actual mode */
+ return ev;
+ }
+
+ s->idark_int_time[3] = ADARKINT_MAX; /* 2.0 */
+ a1logd(p->log,3,"Doing adaptive interpolated black calibration, dcaltime %f, idark_int_time[3] %f, nummeas %d, gainmode %d\n", s->dcaltime, s->idark_int_time[3], nummeas, 1);
+ nummeas = munki_comp_nummeas(p, s->dcaltime, s->idark_int_time[3]);
+ if ((ev = munki_dark_measure(p, s->idark_data[3], nummeas, &s->idark_int_time[3], 1))
+ != MUNKI_OK) {
+ m->mmode = mmode; /* Restore actual mode */
+ return ev;
+ }
}
-#endif /* USE_HIGH_GAIN_MODE */
munki_prepare_idark(p);
@@ -1263,12 +1307,7 @@ if (ss->dark_int_time != s->dark_int_time
|| ss->dark_gain_mode != s->dark_gain_mode)
a1logd(p->log,1,"copying cal to mode with different cal/gain mode: %d -> %d\n",s->mode, ss->mode);
#endif
-#ifdef USE_HIGH_GAIN_MODE
- for (j = 0; j < 4; j++)
-#else
- for (j = 0; j < 2; j++)
-#endif
- {
+ for (j = 0; j < (s->auto_gain ? 4 : 2); j++) {
ss->idark_int_time[j] = s->idark_int_time[j];
#ifndef NEVER // ~~99
if (ss->idark_int_time[j] != s->idark_int_time[j])
@@ -1318,30 +1357,30 @@ if (ss->idark_int_time[j] != s->idark_int_time[j])
ddumpdarkm = 0;
#endif
-#ifdef USE_HIGH_GAIN_MODE
- // fprintf(stderr,"High gain offsets, base:\n");
- // plot_raw(s->idark_data[2]);
- // fprintf(stderr,"High gain offsets, multiplier:\n");
- // plot_raw(s->idark_data[3]);
+ if (s->auto_gain) {
+// fprintf(stderr,"High gain offsets, base:\n");
+// plot_raw(s->idark_data[2]);
+// fprintf(stderr,"High gain offsets, multiplier:\n");
+// plot_raw(s->idark_data[3]);
- for (tinttime = m->min_int_time; ; tinttime *= 2.0) {
- if (tinttime >= m->max_int_time)
- tinttime = m->max_int_time;
-
- nummeas = munki_comp_nummeas(p, s->dcaltime, tinttime);
- if ((ev = munki_dark_measure(p, ref, nummeas, &tinttime, 1)) != MUNKI_OK) {
- m->mmode = mmode; /* Restore actual mode */
- return ev;
- }
- munki_interp_dark(p, interp, tinttime, 1);
+ for (tinttime = m->min_int_time; ; tinttime *= 2.0) {
+ if (tinttime >= m->max_int_time)
+ tinttime = m->max_int_time;
+
+ nummeas = munki_comp_nummeas(p, s->dcaltime, tinttime);
+ if ((ev = munki_dark_measure(p, ref, nummeas, &tinttime, 1)) != MUNKI_OK) {
+ m->mmode = mmode; /* Restore actual mode */
+ return ev;
+ }
+ munki_interp_dark(p, interp, tinttime, 1);
#ifdef DEBUG
- printf("High gain, int time %f:\n",tinttime);
- plot_raw2(ref, interp);
+ printf("High gain, int time %f:\n",tinttime);
+ plot_raw2(ref, interp);
#endif
- if ((tinttime * 1.1) > m->max_int_time)
- break;
- }
-#endif /* USE_HIGH_GAIN_MODE */
+ if ((tinttime * 1.1) > m->max_int_time)
+ break;
+ }
+ }
}
#endif /* TEST_DARK_INTERP */
@@ -1370,27 +1409,24 @@ if (ss->idark_int_time[j] != s->idark_int_time[j])
return ev;
}
-#ifdef USE_HIGH_GAIN_MODE
- s->idark_int_time[2] = s->inttime;
- nummeas = munki_comp_nummeas(p, s->dcaltime, s->idark_int_time[2]);
- a1logd(p->log,3,"Doing adaptive scan black calibration, dcaltime %f, idark_int_time[2] %f, nummeas %d, gainmode %d\n", s->dcaltime, s->idark_int_time[2], nummeas, s->gainmode);
- if ((ev = munki_dark_measure(p, s->idark_data[2], nummeas, &s->idark_int_time[2], 1))
- != MUNKI_OK) {
- m->mmode = mmode; /* Restore actual mode */
- return ev;
+ if (s->auto_gain) {
+ s->idark_int_time[2] = s->inttime;
+ nummeas = munki_comp_nummeas(p, s->dcaltime, s->idark_int_time[2]);
+ a1logd(p->log,3,"Doing adaptive scan black calibration, dcaltime %f, idark_int_time[2] %f, nummeas %d, gainmode %d\n", s->dcaltime, s->idark_int_time[2], nummeas, s->gainmode);
+ if ((ev = munki_dark_measure(p, s->idark_data[2], nummeas, &s->idark_int_time[2], 1))
+ != MUNKI_OK) {
+ m->mmode = mmode; /* Restore actual mode */
+ return ev;
+ }
}
-#endif
s->idark_valid = 1;
s->iddate = cdate;
-#ifdef USE_HIGH_GAIN_MODE
- if (s->gainmode) {
+ if (s->auto_gain && s->gainmode) {
for (j = -1; j < m->nraw; j++)
s->dark_data[j] = s->idark_data[2][j];
- } else
-#endif
- {
+ } else {
for (j = -1; j < m->nraw; j++)
s->dark_data[j] = s->idark_data[0][j];
}
@@ -1418,12 +1454,8 @@ if (ss->idark_int_time[j] != s->idark_int_time[j])
ss->iddate = s->iddate;
ss->dark_int_time = s->dark_int_time;
ss->dark_gain_mode = s->dark_gain_mode;
-#ifdef USE_HIGH_GAIN_MODE
- for (j = 0; j < 4; j += 2)
-#else
- for (j = 0; j < 2; j += 2)
-#endif
- {
+
+ for (j = 0; j < (s->auto_gain ? 4 : 2); j += 2) {
ss->idark_int_time[j] = s->idark_int_time[j];
for (k = -1; k < m->nraw; k++)
ss->idark_data[j][k] = s->idark_data[j][k];
@@ -1724,6 +1756,154 @@ int icoms2munki_err(int se) {
return MUNKI_OK;
}
+/* - - - - - - - - - - - - - - - - */
+/* Measure a display update delay. It is assumed that a */
+/* white to black change has been made to the displayed color, */
+/* and this will measure the time it took for the update to */
+/* be noticed by the instrument, up to 0.6 seconds. */
+/* inst_misread will be returned on failure to find a transition to black. */
+#define NDMXTIME 0.7 /* Maximum time to take */
+#define NDSAMPS 500 /* Debug samples */
+
+typedef struct {
+ double sec;
+ double rgb[3];
+ double tot;
+} i1rgbdsamp;
+
+munki_code munki_imp_meas_delay(
+munki *p,
+int *msecdelay) { /* Return the number of msec */
+ munki_code ev = MUNKI_OK;
+ munkiimp *m = (munkiimp *)p->m;
+ munki_state *s = &m->ms[m->mmode];
+ int i, j, k, mm;
+ double **multimeas; /* Spectral measurements */
+ int nummeas;
+ double rgbw[3] = { 610.0, 520.0, 460.0 };
+ double ucalf = 1.0; /* usec_time calibration factor */
+ double inttime;
+ i1rgbdsamp *samp;
+ double stot, etot, del, thr;
+ double etime;
+ int isdeb;
+
+ /* Read the samples */
+ inttime = m->min_int_time;
+ nummeas = (int)(NDMXTIME/inttime + 0.5);
+ multimeas = dmatrix(0, nummeas-1, -1, m->nwav-1);
+ if ((samp = (i1rgbdsamp *)calloc(sizeof(i1rgbdsamp), nummeas)) == NULL) {
+ a1logd(p->log, 1, "munki_meas_delay: malloc failed\n");
+ return MUNKI_INT_MALLOC;
+ }
+
+//printf("~1 %d samples at %f int\n",nummeas,inttime);
+ if ((ev = munki_read_patches_all(p, multimeas, nummeas, &inttime, 0)) != inst_ok) {
+ free_dmatrix(multimeas, 0, nummeas-1, 0, m->nwav-1);
+ free(samp);
+ return ev;
+ }
+
+ /* Convert the samples to RGB */
+ for (i = 0; i < nummeas; i++) {
+ samp[i].sec = i * inttime;
+ samp[i].rgb[0] = samp[i].rgb[1] = samp[i].rgb[2] = 0.0;
+ for (j = 0; j < m->nwav; j++) {
+ double wl = XSPECT_WL(m->wl_short, m->wl_long, m->nwav, j);
+
+//printf("~1 multimeas %d %d = %f\n",i, j, multimeas[i][j]);
+ for (k = 0; k < 3; k++) {
+ double tt = (double)(wl - rgbw[k]);
+ tt = (50.0 - fabs(tt))/50.0;
+ if (tt < 0.0)
+ tt = 0.0;
+ samp[i].rgb[k] += sqrt(tt) * multimeas[i][j];
+ }
+ }
+ samp[i].tot = samp[i].rgb[0] + samp[i].rgb[1] + samp[i].rgb[2];
+ }
+ free_dmatrix(multimeas, 0, nummeas-1, 0, m->nwav-1);
+
+ a1logd(p->log, 3, "munki_measure_refresh: Read %d samples for refresh calibration\n",nummeas);
+
+#ifdef PLOT_UPDELAY
+ /* Plot the raw sensor values */
+ {
+ double xx[NDSAMPS];
+ double y1[NDSAMPS];
+ double y2[NDSAMPS];
+ double y3[NDSAMPS];
+ double y4[NDSAMPS];
+
+ for (i = 0; i < nummeas && i < NDSAMPS; i++) {
+ xx[i] = samp[i].sec;
+ y1[i] = samp[i].rgb[0];
+ y2[i] = samp[i].rgb[1];
+ y3[i] = samp[i].rgb[2];
+ y4[i] = samp[i].tot;
+//printf("%d: %f -> %f\n",i,samp[i].sec, samp[i].tot);
+ }
+ printf("Display update delay measure sensor values and time (sec)\n");
+ do_plot6(xx, y1, y2, y3, y4, NULL, NULL, nummeas);
+ }
+#endif
+
+ /* Over the first 100msec, locate the maximum value */
+ etime = samp[nummeas-1].sec;
+ stot = -1e9;
+ for (i = 0; i < nummeas; i++) {
+ if (samp[i].tot > stot)
+ stot = samp[i].tot;
+ if (samp[i].sec > 0.1)
+ break;
+ }
+
+ /* Over the last 100msec, locate the maximum value */
+ etime = samp[nummeas-1].sec;
+ etot = -1e9;
+ for (i = nummeas-1; i >= 0; i--) {
+ if (samp[i].tot > etot)
+ etot = samp[i].tot;
+ if ((etime - samp[i].sec) > 0.1)
+ break;
+ }
+
+ del = stot - etot;
+ thr = etot + 0.30 * del; /* 30% of transition threshold */
+
+#ifdef PLOT_UPDELAY
+ a1logd(p->log, 0, "munki_meas_delay: start tot %f end tot %f del %f, thr %f\n", stot, etot, del, thr);
+#endif
+
+ /* Check that there has been a transition */
+ if (del < 5.0) {
+ free(samp);
+ a1logd(p->log, 1, "munki_meas_delay: can't detect change from white to black\n");
+ return MUNKI_RD_NOTRANS_FOUND;
+ }
+
+ /* Locate the time at which the values are above the end values */
+ for (i = nummeas-1; i >= 0; i--) {
+ if (samp[i].tot > thr)
+ break;
+ }
+ if (i < 0) /* Assume the update was so fast that we missed it */
+ i = 0;
+
+ a1logd(p->log, 2, "munki_meas_delay: stoped at sample %d time %f\n",i,samp[i].sec);
+
+ *msecdelay = (int)(samp[i].sec * 1000.0 + 0.5);
+
+#ifdef PLOT_UPDELAY
+ a1logd(p->log, 0, "munki_meas_delay: returning %d msec\n",*msecdelay);
+#endif
+ free(samp);
+
+ return MUNKI_OK;
+}
+#undef NDSAMPS
+#undef NDMXTIME
+
/* - - - - - - - - - - - - - - - - */
/* Measure a patch or strip or flash in the current mode. */
@@ -2276,7 +2456,7 @@ munki_code munki_imp_meas_refrate(
tt = (40.0 - fabs(tt))/40.0;
if (tt < 0.0)
tt = 0.0;
- samp[i].rgb[k] += tt * multimeas[i][j];
+ samp[i].rgb[k] += sqrt(tt) * multimeas[i][j];
}
}
}
@@ -3461,7 +3641,8 @@ munki_code munki_dark_measure_2(
free_dmatrix(multimes, 0, nummeas-1, -1, m->nraw-1);
#ifdef PLOT_DEBUG
- printf("Average absolute sensor readings, average = %f, darkthresh %f\n",sensavg,darkthresh);
+ a1logd(p->log,4,"dark_measure_2: Avg abs. sensor readings = %f, sh/darkthresh %f\n",sensavg,darkthresh);
+ printf("sens data:\n");
plot_raw(sens);
#endif
@@ -3681,6 +3862,7 @@ munki_code munki_whitemeasure(
a1logd(p->log,3,"Average absolute sensor readings, avg %f, max %f, darkth %f satth %f\n",
sensavg,maxval,darkthresh,trackmax[2]);
#ifdef PLOT_DEBUG
+ printf("absraw whitemeas:\n");
plot_raw(absraw);
#endif
}
@@ -3720,7 +3902,7 @@ munki_code munki_compute_wav_whitemeas(
munki_absraw_to_abswav1(p, 1, &abswav1, &absraw);
#ifdef PLOT_DEBUG
- printf("Converted to wavelengths std res:\n");
+ printf("White meas converted to wavelengths std res:\n");
plot_wav1(m, abswav1);
#endif
}
@@ -3810,7 +3992,7 @@ munki_code munki_ledtemp_whitemeasure(
*reftemp = 0.5 * (ledtemp[0] + ledtemp[nummeas-1]);
#ifdef PLOT_DEBUG
- printf("Dark data:\n");
+ printf("Ledtemp Dark data:\n");
plot_raw(s->dark_data);
#endif
@@ -3961,8 +4143,7 @@ munki_code munki_ledtemp_comp(
return MUNKI_OK;
}
-/* Take a measurement reading using the current mode, part 1 */
-/* Converts to completely processed output readings. */
+/* Trigger measure and gather raw readings using the current mode. */
munki_code munki_read_patches_1(
munki *p,
int ninvmeas, /* Number of extra invalid measurements at start */
@@ -3972,7 +4153,7 @@ munki_code munki_read_patches_1(
int gainmode, /* Gain mode to use, 0 = normal, 1 = high */
int *nmeasuered, /* Number actually measured (excluding ninvmeas) */
unsigned char *buf, /* Raw USB reading buffer */
- unsigned int bsize
+ unsigned int bsize /* Raw USB readings buffer size in bytes */
) {
munki_code ev = MUNKI_OK;
munkiimp *m = (munkiimp *)p->m;
@@ -4002,8 +4183,8 @@ munki_code munki_read_patches_1(
return ev;
}
-/* Take a measurement reading using the current mode, part 2 */
-/* Converts to completely processed output readings. */
+/* Given a buffer full of raw USB values, process them into */
+/* completely processed spectral output patch readings. */
munki_code munki_read_patches_2(
munki *p,
double *duration, /* Return flash duration in seconds */
@@ -4014,7 +4195,7 @@ munki_code munki_read_patches_2(
int ninvmeas, /* Number of extra invalid measurements at start */
int nummeas, /* Number of actual measurements */
unsigned char *buf, /* Raw USB reading buffer */
- unsigned int bsize
+ unsigned int bsize /* Raw USB reading buffer size */
) {
munki_code ev = MUNKI_OK;
munkiimp *m = (munkiimp *)p->m;
@@ -4049,32 +4230,6 @@ munki_code munki_read_patches_2(
munki_sub_raw_to_absraw(p, nummeas, inttime, gainmode, multimes, s->dark_data,
&darkthresh, 1, NULL);
-#ifdef PLOT_TEMPCOMP
- /* Plot the raw spectra, 6 at a time */
- if (s->reflective) {
- int i, j, k;
- double *indx;
- double **mod;
- indx = dvectorz(0, nummeas-1);
- mod = dmatrix(0, 5, 0, nummeas-1);
- for (i = 0; i < nummeas; i++)
- indx[i] = (double)i;
-
-// for (j = 0; (j+5) < m->nraw; j += 6) {
- for (j = 50; (j+5) < 56; j += 6) {
- for (i = 0; i < nummeas; i++) {
- for (k = j; k < (j + 6); k++) {
- mod[k-j][i] = multimes[i][k];
- }
- }
-
- printf("Before temp comp, bands %d - %d\n",j, j+5);
- do_plot6(indx, mod[0], mod[1], mod[2], mod[3], mod[4], mod[5], nummeas);
- }
- free_dvector(indx, 0, nummeas-1);
- free_dmatrix(mod, 0, 5, 0, nummeas-1);
- }
-#endif
#ifdef DUMP_SCANV
/* Dump raw scan readings to a file "mkdump.txt" */
{
@@ -4097,8 +4252,38 @@ munki_code munki_read_patches_2(
#endif
#ifdef ENABLE_LEDTEMPC
+
+
/* Do LED temperature compensation of absraw values */
if (s->reflective) {
+
+#ifdef PLOT_TEMPCOMP
+ /* Plot the raw spectra, 6 at a time */
+ {
+ int i, j, k;
+ double *indx;
+ double **mod;
+ indx = dvectorz(0, nummeas-1);
+ mod = dmatrix(0, 5, 0, nummeas-1);
+ for (i = 0; i < nummeas; i++)
+ indx[i] = (double)i;
+
+// for (j = 0; (j+5) < m->nraw; j += 6) {
+ for (j = 50; (j+5) < 56; j += 6) {
+ for (i = 0; i < nummeas; i++) {
+ for (k = j; k < (j + 6); k++) {
+ mod[k-j][i] = multimes[i][k];
+ }
+ }
+
+ printf("Before temp comp, bands %d - %d\n",j, j+5);
+ do_plot6(indx, mod[0], mod[1], mod[2], mod[3], mod[4], mod[5], nummeas);
+ }
+ free_dvector(indx, 0, nummeas-1);
+ free_dmatrix(mod, 0, 5, 0, nummeas-1);
+ }
+#endif
+ /* Do the LED temperature compensation */
if ((ev = munki_ledtemp_comp(p, multimes, ledtemp, nummeas, s->reftemp, s->iwhite_data))
!= MUNKI_OK) {
free_dvector(ledtemp, 0, nummeas-1);
@@ -4177,7 +4362,7 @@ munki_code munki_read_patches_2(
/* Recognise the required number of ref/trans patch locations, */
/* and average the measurements within each patch. */
if ((ev = munki_extract_patches_multimeas(p, &rv, absraw, numpatches, multimes,
- nummeas, inttime)) != MUNKI_OK) {
+ nummeas, inttime)) != MUNKI_OK) {
free_dvector(ledtemp, 0, nummeas-1);
free_dmatrix(multimes, 0, nummeas-1, -1, m->nraw-1);
free_dmatrix(absraw, 0, numpatches-1, -1, m->nraw-1);
@@ -4195,6 +4380,52 @@ munki_code munki_read_patches_2(
return MUNKI_RD_READINCONS;
}
+#ifdef ENABLE_REFSTRAYC /* Enable Reflective stray light compensation */
+ if (s->reflective) {
+ int i, j;
+# if defined(PLOT_DEBUG) || defined(DUMP_BKLED)
+ double xx[140];
+ double yy[3][140];
+# endif
+
+ double fact = REFSTRAYC_FACTOR; /* Slightly conservative */
+
+ for (i = 0; i < numpatches; i++) {
+
+ for (j = 0; j < m->nraw; j++) {
+# if defined(PLOT_DEBUG) || defined(DUMP_BKLED)
+ yy[0][j] = absraw[i][j];
+# endif
+ absraw[i][j] -= fact * s->white_data[j];
+# if defined(PLOT_DEBUG) || defined(DUMP_BKLED)
+ xx[j] = j;
+ yy[1][j] = fact * s->white_data[j];
+ yy[2][j] = absraw[i][j];
+# endif
+ }
+# ifdef PLOT_DEBUG
+ printf("Before/After subtracting stray ref. light %d:\n",i);
+ do_plot6(xx, yy[0], yy[1], yy[2], NULL, NULL, NULL, m->nraw);
+# endif
+# ifdef DUMP_BKLED /* Save REFSTRAYC & REFLEDNOISE comp plot to "refbk1.txt" & "refbk2.txt" */
+ {
+ xspect sp[3];
+ for (i = 0; i < 3; i++) {
+ sp[i].spec_n = 128;
+ sp[i].spec_wl_short = 0.0;
+ sp[i].spec_wl_long = 127.0;
+ sp[i].norm = 1.0;
+ for (j = 0; j < 128; j++)
+ sp[i].spec[j] = yy[i][j];
+ }
+ write_nxspect("refbk2.txt", sp, 3, 0);
+# pragma message("######### munki DUMP_BKLED enabled! ########")
+ }
+# endif /* DUMP_BKLED */
+ }
+ }
+#endif /* ENABLE_REFSTRAYC */
+
/* Convert an absraw array from raw wavelengths to output wavelenths */
munki_absraw_to_abswav(p, numpatches, specrd, absraw);
free_dmatrix(absraw, 0, numpatches-1, -1, m->nraw-1);
@@ -4243,9 +4474,11 @@ munki_code munki_read_patches_2(
return ev;
}
-/* Take a measurement reading using the current mode, part 2a */
-/* Converts to completely processed output readings, */
+/* Given a buffer full of raw USB values, process them into */
+/* completely processed spectral output readings, */
/* but don't average together or extract patches or flash. */
+/* This is used for delay & refresh rate measurement. */
+/* !! Doesn't do LED temperature compensation for reflective !! */
/* (! Note that we aren't currently detecting saturation here!) */
munki_code munki_read_patches_2a(
munki *p,
@@ -4309,7 +4542,8 @@ munki_code munki_read_patches_2a(
/* Take a measurement reading using the current mode (combined parts 1 & 2a) */
/* Converts to completely processed output readings, without averaging or extracting */
-/* sample patches. */
+/* sample patches, for emissive measurement mode. */
+/* This is used for delay & refresh rate measurement. */
munki_code munki_read_patches_all(
munki *p,
double **specrd, /* Return array [numpatches][nwav] of spectral reading values */
@@ -4535,7 +4769,7 @@ munki_code munki_sens_to_raw(
int ninvalid, /* Number of initial invalid readings to skip */
int nummeas, /* Number of readings measured */
double satthresh, /* Saturation threshold to trigger error in raw units (if > 0.0) */
- double *pdarkthresh /* Return a dark threshold value */
+ double *pdarkthresh /* Return a dark threshold value = sheilded cell values */
) {
munkiimp *m = (munkiimp *)p->m;
int i, j, k;
@@ -4551,11 +4785,6 @@ munki_code munki_sens_to_raw(
return MUNKI_INT_ASSERT;
}
- if (ninvalid >= nummeas) {
- a1logw(p->log,"ninvalid %d is >= nummeas %d!\n",ninvalid,nummeas);
- return MUNKI_INT_ASSERT;
- }
-
if (ninvalid > 0)
a1logd(p->log, 4, "munki_sens_to_raw: Skipping %d invalid readings\n",ninvalid);
@@ -4624,12 +4853,13 @@ void munki_sub_raw_to_absraw(
double inttime, /* Integration time used */
int gainmode, /* Gain mode, 0 = normal, 1 = high */
double **absraw, /* Source/Desination array [-1 nraw] */
- double *sub, /* Value to subtract [-1 nraw] */
+ double *sub, /* Value to subtract [-1 nraw] (ie. cal dark data) */
double *trackmax, /* absraw values that should be offset the same as max */
int ntrackmax, /* Number of trackmax values */
double *maxv /* If not NULL, return the maximum value */
) {
munkiimp *m = (munkiimp *)p->m;
+ munki_state *s = &m->ms[m->mmode];
int npoly; /* Number of linearisation coefficients */
double *polys; /* the coeficients */
double scale; /* Absolute scale value */
@@ -4639,6 +4869,11 @@ void munki_sub_raw_to_absraw(
double rawmax, maxval = -1e38;
double maxzero = 0.0;
int i, j, k;
+
+ /* Heusristic correction for LED interference bump for 0.018 secs int_time */
+ int pos[] = { 0, 20, 56, 62, 75, 127 };
+// double off[] = { 0.7, 0.0, 0.6, -0.9, -1.2, -0.7 };
+ double off[] = { 0.7, 0.0, 0.6, -0.9, -0.8, -0.5 };
if (gainmode) { /* High gain */
npoly = m->nlin1; /* Encodes gain too */
@@ -4681,18 +4916,129 @@ void munki_sub_raw_to_absraw(
a1logd(p->log,4,"Black shielded value = %f, Reading shielded value = %f\n",sub[-1], avgscell);
- /* Compute the adjusted black */
- for (j = 0; j < m->nraw; j++) {
-#ifdef NEVER
- /* simple additive correction */
-# pragma message("######### munki Simple shielded cell temperature correction! ########")
- asub[j] = sub[j] + avgscell - sub[-1];
+ /* Compute the adjusted black for each band */
+ if (s->reflective) {
+
+ /* It seems that having the LED on shifts the shielded cell values */
+ /* by about 2.5, and this stuffs up the reflective measurement. */
+ /* This seems to be from the LED PWM driver, which perhaps */
+ /* is synchronous to the sensor clock, and so switches */
+ /* at a certain point in the transfer of data from the sensor. */
+ /* The result is a step up from 0-60, and then down from 61-128. */
+ /* Perhaps altering the LED PWM setting and seeing if this point */
+ /* shifts would be a way of confirming this ? */
+ /* There is also some stray light reflected into the sensor */
+ /* from the LED, but due to the LED step, the sensor reading is */
+ /* less than the dark data at some wavelengths. */
+ /* The details of the LED step seem to be integration time dependent, */
+ /* but decresing the scanning rate therebye increasing integration */
+ /* time and light level reduces the impact of this error. */
+
+ /* Since we do an on the fly black measurement before each */
+ /* reflective measurement, ignoring the shielded cell values */
+ /* shouldn't affect accuracy so much. */
+
+#ifdef ENABLE_REFLEDINTER
+ /* A heuristic to correct for the LED noise. */
+ /* This is only valid for int_time of 0.0182 secs, */
+ /* and it's not clear how well it works across different */
+ /* temperatures or examples of the ColorMunki. */
+ /* in another revision ?? */
+ for (j = 0; j < m->nraw; j++) {
+ int ix;
+ double bl, val;
+
+ for (ix = 0; ; ix++) {
+ if (j >= pos[ix] && j <= pos[ix+1])
+ break;
+ }
+ bl = (j - pos[ix])/((double)pos[ix+1] - pos[ix]);
+ val = (1.0 - bl) * off[ix] + bl * off[ix+1];
+ asub[j] = sub[j] + val;
+ }
+#else
+ for (j = 0; j < m->nraw; j++)
+ asub[j] = sub[j]; /* Just use the calibration dark data */
+#endif
+
+ } else {
+ /* No LED on operation - use sheilded cell values */
+ for (j = 0; j < m->nraw; j++) {
+#ifdef ENABLE_BKDRIFTC
+
+# ifdef HEURISTIC_BKDRIFTC
+ /* heuristic scaled correction */
+ asub[j] = zero - (zero - sub[j]) * (zero - avgscell)/(zero - sub[-1]);
+# else
+ /* simple additive correction */
+# pragma message("######### munki Simple shielded cell temperature correction! ########")
+ asub[j] = sub[j] + (avgscell - sub[-1]);
+# endif
#else
- /* heuristic scaled correction */
- asub[j] = zero - (zero - sub[j]) * (zero - avgscell)/(zero - sub[-1]);
+# pragma message("######### munki No shielded cell temperature correction! ########")
+ asub[j] = sub[j]; /* Just use the calibration dark data */
#endif
+ }
}
-
+
+#if defined(PLOT_DEBUG) || defined(DUMP_BKLED)
+ {
+ double xx[130];
+ double yy[3][130];
+
+ for (j = -1; j < m->nraw+1; j++)
+ yy[0][j+1] = 0.0;
+
+ for (i = 0; i < nummeas; i++) {
+ for (j = -1; j < m->nraw; j++)
+ yy[0][j+1] += absraw[i][j];
+ }
+ for (j = -1; j < m->nraw; j++)
+ yy[0][j+1] /= (double)nummeas;
+
+ for (j = -1; j < m->nraw; j++)
+ yy[1][j+1]= sub[j];
+
+ /* Show what ENABLE_REFLEDINTER would do */
+ for (j = 0; j < m->nraw; j++) {
+ int ix;
+ double bl, val;
+
+ for (ix = 0; ; ix++) {
+ if (j >= pos[ix] && j <= pos[ix+1])
+ break;
+ }
+ bl = (j - pos[ix])/((double)pos[ix+1] - pos[ix]);
+ val = (1.0 - bl) * off[ix] + bl * off[ix+1];
+ yy[2][j+1] = yy[0][j+1] - val;
+ }
+ yy[2][0] = yy[0][0];
+
+ for (j = -1; j < m->nraw; j++)
+ xx[j+1] = (double)j;
+
+ xx[0]= -10.0;
+
+# ifdef PLOT_DEBUG
+ printf("sub_raw_to_absraw %d samp avg - dark ref:\n",nummeas);
+ do_plot(xx, yy[0], yy[1], yy[2], 129);
+# endif
+# ifdef DUMP_BKLED
+ {
+ xspect sp[3];
+ for (i = 0; i < 3; i++) {
+ sp[i].spec_n = 128;
+ sp[i].spec_wl_short = 0.0;
+ sp[i].spec_wl_long = 127.0;
+ sp[i].norm = 1.0;
+ for (j = 0; j < 128; j++)
+ sp[i].spec[j] = yy[i][j+1];
+ }
+ write_nxspect("refbk1.txt", sp, 3, 0);
+ }
+# endif /* DUMP_BKLED */
+ }
+#endif /* PLOT_DEBUG || DUMP_BKLED */
/* For each measurement */
for (i = 0; i < nummeas; i++) {
@@ -4701,6 +5047,7 @@ void munki_sub_raw_to_absraw(
for (j = 0; j < m->nraw; j++) {
rval = absraw[i][j];
+
sval = rval - asub[j]; /* Make zero based */
#ifdef ENABLE_NONLINCOR
@@ -4825,7 +5172,7 @@ int munki_average_multimeas(
#define BH 105 /* End */
#define BW 5 /* Width */
-/* Record of possible patch */
+/* Record of possible patch within a reading buffer */
typedef struct {
int ss; /* Start sample index */
int no; /* Number of samples */
@@ -4866,7 +5213,7 @@ munki_code munki_extract_patches_multimeas(
double white_avg; /* Average of (aproximate) white data */
int rv = 0;
double patch_cons_thr = PATCH_CONS_THR * m->scan_toll_ratio;
-#ifdef PATREC_DEBUG
+#ifdef PLOT_PATREC
double **plot;
#endif
@@ -4888,7 +5235,7 @@ munki_code munki_extract_patches_multimeas(
maxval[j] = 1.0;
}
-#ifdef PATREC_DEBUG
+#ifdef PLOT_PATREC
/* Plot out 6 lots of 6 values each */
plot = dmatrixz(0, 6, 0, nummeas-1);
// for (j = 1; j < (m->nraw-6); j += 6) /* Plot all the bands */
@@ -4906,6 +5253,16 @@ munki_code munki_extract_patches_multimeas(
}
#endif
+#ifdef NEVER
+ /* Plot the shielded cell value */
+ for (i = 0; i < nummeas; i++) {
+ plot[0][i] = multimeas[i][-1];
+ plot[6][i] = (double)i;
+ }
+ printf("Sheilded values\n");
+ do_plot6(plot[6], plot[0], NULL, NULL, NULL, NULL, NULL, nummeas);
+#endif
+
sslope = dmatrixz(0, nummeas-1, -1, m->nraw-1);
slope = dvectorz(0, nummeas-1);
fslope = dvectorz(0, nummeas-1);
@@ -5144,7 +5501,7 @@ munki_code munki_extract_patches_multimeas(
}
#endif
-#ifdef PATREC_DEBUG
+#ifdef PLOT_PATREC
printf("Slope filter output\n");
for (i = 0; i < nummeas; i++) {
int jj;
@@ -5196,7 +5553,8 @@ munki_code munki_extract_patches_multimeas(
}
a1logd(p->log,7,"Number of patches = %d\n",npat);
- /* We don't count the first and last patches, as we assume they are white leader */
+ /* We don't count the first and last patches, as we assume they are white leader. */
+ /* (They are marked !use in list anyway) */
if (npat < (tnpatch + 2)) {
free_ivector(sizepop, 0, nummeas-1);
free_dvector(slope, 0, nummeas-1);
@@ -5214,7 +5572,7 @@ munki_code munki_extract_patches_multimeas(
}
avglegth /= (double)npat;
-#ifdef PATREC_DEBUG
+#ifdef PLOT_PATREC
for (i = 0; i < npat; i++) {
printf("Raw patch %d, start %d, length %d\n",i, pat[i].ss, pat[i].no);
}
@@ -5227,7 +5585,7 @@ munki_code munki_extract_patches_multimeas(
/* Locate the median potential patch width */
for (j = 0, i = 0; i < nummeas; i++) {
j += sizepop[i];
- if (j >= ((npat-2)/2))
+ if (j > ((npat-2)/2))
break;
}
median = (double)i;
@@ -5308,7 +5666,7 @@ munki_code munki_extract_patches_multimeas(
return MUNKI_RD_NOTENOUGHPATCHES;
}
-#ifdef PATREC_DEBUG
+#ifdef PLOT_PATREC
printf("Got %d patches out of potentional %d:\n",tnpatch, npat);
printf("Average patch legth %f\n",avglegth);
for (i = 1; i < (npat-1); i++) {
@@ -5325,15 +5683,15 @@ munki_code munki_extract_patches_multimeas(
if (pat[k].use == 0)
continue;
-
nno = (pat[k].no * 3)/4;
+// nno = (pat[k].no * 2)/3; // experimental tighter trim
trim = (pat[k].no - nno)/2;
pat[k].ss += trim;
pat[k].no = nno;
}
-#ifdef PATREC_DEBUG
+#ifdef PLOT_PATREC
printf("After trimming got:\n");
for (i = 1; i < (npat-1); i++) {
if (pat[i].use == 0)
@@ -5351,7 +5709,8 @@ munki_code munki_extract_patches_multimeas(
slope[i] = 0.0;
}
- printf("Trimmed output:\n");
+ printf("Trimmed output - averaged bands:\n");
+ /* Plot box averaged bands */
for (i = 0; i < nummeas; i++) {
int jj;
for (jj = 0, j = BL; jj < 6 && j < BH; jj++, j += ((BH-BL)/6)) {
@@ -5364,8 +5723,27 @@ munki_code munki_extract_patches_multimeas(
for (i = 0; i < nummeas; i++)
plot[6][i] = (double)i;
do_plot6(plot[6], slope, plot[0], plot[1], plot[2], plot[3], plot[4], nummeas);
+
+#ifdef NEVER
+ /* Plot all the bands */
+ printf("Trimmed output - all bands:\n");
+ for (j = 0; j < (m->nraw-5); j += 5) {
+ for (k = 0; k < 5; k ++) {
+ for (i = 0; i < nummeas; i++) {
+ plot[k][i] = multimeas[i][j+k]/maxval[j+k];
+ }
+ }
+ for (i = 0; i < nummeas; i++)
+ plot[6][i] = (double)i;
+ printf("Bands %d - %d\n",j,j+5);
+ do_plot6(plot[6], slope, plot[0], plot[1], plot[2], plot[3], plot[4], nummeas);
+ }
+#endif
+
#endif
+ /* Now compute averaged patch values */
+
/* Compute average of (aproximate) white */
white_avg = 0.0;
for (j = 1; j < (m->nraw-1); j++)
@@ -5429,14 +5807,14 @@ munki_code munki_extract_patches_multimeas(
if (flags != NULL)
*flags = rv;
-#ifdef PATREC_DEBUG
+#ifdef PLOT_PATREC
free_dmatrix(plot, 0, 6, 0, nummeas-1);
#endif
free_dvector(slope, 0, nummeas-1);
free_ivector(sizepop, 0, nummeas-1);
free_dvector(maxval, -1, m->nraw-1);
- free(pat);
+ free(pat); /* Otherwise caller will have to do it */
a1logd(p->log,3,"munki_extract_patches_multimeas done, sat = %s, inconsist = %s\n",
rv & 2 ? "true" : "false", rv & 1 ? "true" : "false");
@@ -5475,7 +5853,7 @@ munki_code munki_extract_patches_flash(
double *aavg; /* ambient average [-1 nraw] */
double finttime; /* Flash integration time */
int rv = 0;
-#ifdef PATREC_DEBUG
+#ifdef PLOT_PATREC
double **plot;
#endif
@@ -5511,7 +5889,7 @@ munki_code munki_extract_patches_flash(
thresh = (3.0 * mean + maxval)/4.0;
a1logd(p->log,7,"munki_extract_patches_flash band %d minval %f maxval %f, mean = %f, thresh = %f\n",maxband,minval,maxval,mean, thresh);
-#ifdef PATREC_DEBUG
+#ifdef PLOT_PATREC
/* Plot out 6 lots of 6 values each */
plot = dmatrixz(0, 6, 0, nummeas-1);
for (j = maxband -3; j>= 0 && j < (m->nraw-6); j += 100) /* Do one set around max */
@@ -5529,7 +5907,7 @@ munki_code munki_extract_patches_flash(
free_dmatrix(plot,0,6,0,nummeas-1);
#endif
-#ifdef PATREC_DEBUG
+#ifdef PLOT_PATREC
/* Plot just the pulses */
{
int start, end;
@@ -5841,8 +6219,11 @@ void munki_scale_specrd(
#ifdef HIGH_RES
/* High res congiguration */
-#undef EXISTING_SHAPE /* Else generate filter shape */
-#undef USE_GAUSSIAN /* Use gaussian filter shape, else lanczos2 */
+#undef EXISTING_SHAPE /* [und] Else generate filter shape */
+#define USE_GAUSSIAN /* [def] Use gaussian filter shape, else lanczos2 */
+
+#define DO_CCDNORM /* [def] Normalise CCD values to original */
+#define DO_CCDNORMAVG /* [unde] Normalise averages rather than per CCD bin */
#ifdef NEVER
/* Plot the matrix coefficients */
@@ -5940,7 +6321,8 @@ static double lanczos2(double wi, double x) {
/* gausian */
wi = wi/(2.0 * sqrt(2.0 * log(2.0))); /* Convert width at half max to std. dev. */
x = x/(sqrt(2.0) * wi);
- y = 1.0/(wi * sqrt(2.0 * DBL_PI)) * exp(-(x * x));
+// y = 1.0/(wi * sqrt(2.0 * DBL_PI)) * exp(-(x * x)); /* Unity area */
+ y = exp(-(x * x)); /* Center at 1.0 */
#else
/* lanczos2 */
@@ -5964,6 +6346,10 @@ static int gcc_bug_fix(int i) {
}
#endif /* APPLE */
+#ifdef SALONEINSTLIB
+# define ONEDSTRAYLIGHTUS
+#endif
+
/* Create high resolution mode references, */
/* Create Reflective if ref nz, else create Emissive */
/* We expect this to be called twice, once for each. */
@@ -5982,6 +6368,8 @@ munki_code munki_create_hr(munki *p, int ref) {
int *mtx_index1, **pmtx_index2, *mtx_index2;
int *mtx_nocoef1, **pmtx_nocoef2, *mtx_nocoef2;
double *mtx_coef1, **pmtx_coef2, *mtx_coef2;
+
+ double min_wl = ref ? WL_REF_MIN : WL_EMIS_MIN;
/* Start with nominal values. May alter these if filters are not unique */
nwav1 = m->nwav1;
@@ -5993,7 +6381,8 @@ munki_code munki_create_hr(munki *p, int ref) {
mtx_index1 = m->rmtx_index1;
mtx_nocoef1 = m->rmtx_nocoef1;
mtx_coef1 = m->rmtx_coef1;
- mtx_index2 = mtx_nocoef2 = NULL;
+ mtx_index2 = NULL;
+ mtx_nocoef2 = NULL;
mtx_coef2 = NULL;
pmtx_index2 = &m->rmtx_index2;
pmtx_nocoef2 = &m->rmtx_nocoef2;
@@ -6002,7 +6391,8 @@ munki_code munki_create_hr(munki *p, int ref) {
mtx_index1 = m->emtx_index1;
mtx_nocoef1 = m->emtx_nocoef1;
mtx_coef1 = m->emtx_coef1;
- mtx_index2 = mtx_nocoef2 = NULL;
+ mtx_index2 = NULL;
+ mtx_nocoef2 = NULL;
mtx_coef2 = NULL;
pmtx_index2 = &m->emtx_index2;
pmtx_nocoef2 = &m->emtx_nocoef2;
@@ -6020,7 +6410,7 @@ munki_code munki_create_hr(munki *p, int ref) {
/* For each matrix value */
sx = mtx_index1[j]; /* Starting index */
- if (jj == 0 && (j+1) < m->nwav1 && sx == mtx_index1[j+1]) { /* Same index */
+ if (j < (m->nwav1-1) && sx == mtx_index1[j+1]) { /* Skip duplicates + last */
// printf("~1 skipping %d\n",j);
wl_short1 += wl_step1;
nwav1--;
@@ -6076,6 +6466,7 @@ munki_code munki_create_hr(munki *p, int ref) {
/* Compute the crossover points between each filter */
for (i = 0; i < (nwav1-1); i++) {
double den, y1, y2, y3, y4, yn, xn; /* Location of intersection */
+ double eps = 1e-6; /* Numerical tollerance */
double besty = -1e6;
/* between filter i and i+1, we want to find the two */
@@ -6085,13 +6476,8 @@ munki_code munki_create_hr(munki *p, int ref) {
for (j = 0; j < (mtx_nocoef1[i]-1); j++) {
for (k = 0; k < (mtx_nocoef1[i+1]-1); k++) {
if (coeff[i][j].ix == coeff[i+1][k].ix
- && coeff[i][j+1].ix == coeff[i+1][k+1].ix
- && coeff[i][j].we > 0.0 && coeff[i][j+1].we > 0.0
- && coeff[i+1][k].we > 0.0 && coeff[i+1][k+1].we > 0.0
- && (( coeff[i][j].we >= coeff[i+1][k].we
- && coeff[i][j+1].we <= coeff[i+1][k+1].we)
- || ( coeff[i][j].we <= coeff[i+1][k].we
- && coeff[i][j+1].we >= coeff[i+1][k+1].we))) {
+ && coeff[i][j+1].ix == coeff[i+1][k+1].ix) {
+
// a1logd(p->log,3,"got it at %d, %d: %d = %d, %d = %d\n",j,k, coeff[i][j].ix, coeff[i+1][k].ix, coeff[i][j+1].ix, coeff[i+1][k+1].ix);
/* Compute the intersection of the two line segments */
@@ -6099,17 +6485,22 @@ munki_code munki_create_hr(munki *p, int ref) {
y2 = coeff[i][j+1].we;
y3 = coeff[i+1][k].we;
y4 = coeff[i+1][k+1].we;
+// a1logd(p->log,3,"y1 %f, y2 %f, y3 %f, y4 %f\n",y1, y2, y3, y4);
den = -y4 + y3 + y2 - y1;
+ if (fabs(den) < eps)
+ continue;
yn = (y2 * y3 - y1 * y4)/den;
xn = (y3 - y1)/den;
+ if (xn < -eps || xn > (1.0 + eps))
+ continue;
// a1logd(p->log,3,"den = %f, yn = %f, xn = %f\n",den,yn,xn);
-// a1logd(p->log,3,"Intersection %d: wav %f, raw %f, wei %f\n",i+1,xp[i+1].wav,xp[i+1].raw,xp[i+1].wei);
if (yn > besty) {
xp[i+1].wav = XSPECT_WL(wl_short1, wl_long1, nwav1, i + 0.5);
xp[i+1].raw = (1.0 - xn) * coeff[i][j].ix + xn * coeff[i][j+1].ix;
xp[i+1].wei = yn;
besty = yn;
-// a1logd(p->log,3,"Found new best y\n");
+// a1logd(p->log,3,"Intersection %d: wav %f, raw %f, wei %f\n",i+1,xp[i+1].wav,xp[i+1].raw,xp[i+1].wei);
+// a1logd(p->log,3,"Found new best y %f\n",yn);
}
// a1logd(p->log,3,"\n");
}
@@ -6151,7 +6542,7 @@ munki_code munki_create_hr(munki *p, int ref) {
}
}
if (j >= mtx_nocoef1[0]) { /* Assert */
- a1logw(p->log,"munki: failed to end crossover\n");
+ a1logw(p->log,"munki: failed to find end crossover\n");
return MUNKI_INT_ASSERT;
}
den = -y4 + y3 + y2 - y1;
@@ -6189,7 +6580,7 @@ munki_code munki_create_hr(munki *p, int ref) {
}
}
if (j >= mtx_nocoef1[nwav1-1]) { /* Assert */
- a1logw(p->log, "munki: failed to end crossover\n");
+ a1logw(p->log, "munki: failed to find end crossover\n");
return MUNKI_INT_ASSERT;
}
den = -y4 + y3 + y2 - y1;
@@ -6397,14 +6788,15 @@ munki_code munki_create_hr(munki *p, int ref) {
{
double fshmax; /* filter shape max wavelength from center */
+#define MXNOWL 500 /* Max hires bands */
#define MXNOFC 64
- munki_fc coeff2[500][MXNOFC]; /* New filter cooefficients */
+ munki_fc coeff2[MXNOWL][MXNOFC]; /* New filter cooefficients */
double twidth;
/* Construct a set of filters that uses more CCD values */
twidth = HIGHRES_WIDTH;
- if (m->nwav2 > 500) { /* Assert */
+ if (m->nwav2 > MXNOWL) { /* Assert */
a1logw(p->log,"High res filter has too many bands\n");
return MUNKI_INT_ASSERT;
}
@@ -6474,7 +6866,7 @@ munki_code munki_create_hr(munki *p, int ref) {
raw2wav->interp(raw2wav, &pp);
w2 = pp.v[0];
-// a1logd(p->log,1,"CCD %d, wl %f\n",i,wl);
+ a1logd(p->log,1,"CCD %d, wl %f - %f\n",i,w1,w2);
/* For each filter */
for (j = 0; j < m->nwav2; j++) {
@@ -6482,6 +6874,10 @@ munki_code munki_create_hr(munki *p, int ref) {
double we;
cwl = m->wl_short2 + (double)j * (m->wl_long2 - m->wl_short2)/(m->nwav2-1.0);
+
+ if (cwl < min_wl) /* Duplicate below this wl */
+ cwl = min_wl;
+
rwl = wl - cwl; /* relative wavelgth to filter */
if (fabs(w1 - cwl) > fshmax && fabs(w2 - cwl) > fshmax)
@@ -6492,8 +6888,12 @@ munki_code munki_create_hr(munki *p, int ref) {
{
int nn;
double lw, ll;
-
- nn = (int)(fabs(w2 - w1)/0.05 + 0.5);
+#ifdef FAST_HIGH_RES_SETUP
+# define FINC 0.2
+#else
+# define FINC 0.05
+#endif
+ nn = (int)(fabs(w2 - w1)/FINC + 0.5);
lw = w1;
#ifdef EXISTING_SHAPE
@@ -6538,7 +6938,21 @@ munki_code munki_create_hr(munki *p, int ref) {
coeff2[j][mtx_nocoef2[j]].ix = i;
coeff2[j][mtx_nocoef2[j]++].we = we;
-// a1logd(p->log,1,"filter %d, cwl %f, rwl %f, ix %d, we %f\n",j,cwl,rwl,i,we);
+ a1logd(p->log,1,"filter %d, cwl %f, rwl %f, ix %d, we %f\n",j,cwl,rwl,i,we);
+ }
+ }
+
+ /* Dump the filter coefficients */
+ if (p->log->debug >= 1) {
+
+ /* For each output wavelength */
+ for (j = 0; j < m->nwav2; j++) {
+
+ a1logd(p->log,1,"filter %d, cwl %f\n",j,XSPECT_WL(m->wl_short2, m->wl_long2, m->nwav2, j));
+ /* For each matrix value */
+ for (k = 0; k < mtx_nocoef2[j]; k++) {
+ a1logd(p->log,1," CCD %d, we %f\n",coeff2[j][k].ix,coeff2[j][k].we);
+ }
}
}
@@ -6575,89 +6989,230 @@ munki_code munki_create_hr(munki *p, int ref) {
}
#endif /* HIGH_RES_PLOT */
+ /* Convert into runtime format */
+ {
+ int xcount;
+
+ if ((*pmtx_index2 = mtx_index2 = (int *)calloc(m->nwav2, sizeof(int))) == NULL) {
+ a1logd(p->log,1,"munki: malloc mtx_index2 failed!\n");
+ return MUNKI_INT_MALLOC;
+ }
+
+ xcount = 0;
+ for (j = 0; j < m->nwav2; j++) {
+ mtx_index2[j] = coeff2[j][0].ix;
+ xcount += mtx_nocoef2[j];
+ }
+
+ if ((*pmtx_coef2 = mtx_coef2 = (double *)calloc(xcount, sizeof(double))) == NULL) {
+ a1logd(p->log,1,"munki: malloc mtx_coef2 failed!\n");
+ return MUNKI_INT_MALLOC;
+ }
+
+ for (i = j = 0; j < m->nwav2; j++)
+ for (k = 0; k < mtx_nocoef2[j]; k++, i++)
+ mtx_coef2[i] = coeff2[j][k].we;
+ }
+
/* Normalise the filters area in CCD space, while maintaining the */
/* total contribution of each CCD at the target too. */
+ /* Hmm. This will wreck super-sample. We should fix it */
+#ifdef DO_CCDNORM /* Normalise CCD values to original */
{
- int ii;
- double tot = 0.0;
- double ccdweight[NSEN_MAX], avgw; /* Weighting determined by cell widths */
- double ccdsum[NSEN_MAX];
+ double x[4], y[4];
+ double avg[2], max[2];
+ double ccdsum[2][128]; /* Target weight/actual for each CCD */
+ double dth[2];
+
+ avg[0] = avg[1] = 0.0;
+ max[0] = max[1] = 0.0;
+ for (j = 0; j < 128; j++) {
+ ccdsum[0][j] = 0.0;
+ ccdsum[1][j] = 0.0;
+ }
- /* Normalize the overall filter weightings */
- for (j = 0; j < m->nwav2; j++)
- for (k = 0; k < mtx_nocoef2[j]; k++)
- tot += coeff2[j][k].we;
- tot /= (double)m->nwav2;
- for (j = 0; j < m->nwav2; j++)
- for (k = 0; k < mtx_nocoef2[j]; k++)
- coeff2[j][k].we /= tot;
+ /* Compute the weighting of each CCD value in the normal output */
+ for (cx = j = 0; j < m->nwav1; j++) { /* For each wavelength */
- /* Determine the relative weights for each CCD */
- avgw = 0.0;
- for (i = 0; i < m->nraw; i++) {
- co pp;
+ /* For each matrix value */
+ sx = mtx_index1[j]; /* Starting index */
+ if (j < (m->nwav1-2) && sx == mtx_index1[j+1]) {
+ cx += mtx_nocoef1[j];
+ continue; /* Skip all duplicate filters */
+ }
+ for (k = 0; k < mtx_nocoef1[j]; k++, cx++, sx++) {
+ ccdsum[0][sx] += mtx_coef1[cx];
+//printf("~1 Norm CCD [%d] %f += [%d] %f\n",sx,ccdsum[0][sx],cx, mtx_coef1[cx]);
+ }
+ }
- pp.p[0] = i - 0.5;
- raw2wav->interp(raw2wav, &pp);
- ccdweight[i] = pp.v[0];
+ /* Compute the weighting of each CCD value in the hires output */
+ for (cx = j = 0; j < m->nwav2; j++) { /* For each wavelength */
- pp.p[0] = i + 0.5;
- raw2wav->interp(raw2wav, &pp);
- ccdweight[i] = fabs(ccdweight[i] - pp.v[0]);
- avgw += ccdweight[i];
+ /* For each matrix value */
+ sx = mtx_index2[j]; /* Starting index */
+ if (j < (m->nwav2-2) && sx == mtx_index2[j+1]) {
+ cx += mtx_nocoef2[j];
+ continue; /* Skip all duplicate filters */
+ }
+ for (k = 0; k < mtx_nocoef2[j]; k++, cx++, sx++) {
+ ccdsum[1][sx] += mtx_coef2[cx];
+//printf("~1 HiRes CCD [%d] %f += [%d] %f\n",sx,ccdsum[1][sx],cx, mtx_coef2[cx]);
+ }
}
- avgw /= 126.0;
- // ~~this isn't right because not all CCD's get used!!
-#ifdef NEVER
- /* Itterate */
- for (ii = 0; ; ii++) {
-
- /* Normalize the individual filter weightings */
- for (j = 0; j < m->nwav2; j++) {
- double err;
- tot = 0.0;
- for (k = 0; k < mtx_nocoef2[j]; k++)
- tot += coeff2[j][k].we;
- err = 1.0 - tot;
+#ifdef HIGH_RES_PLOT
+ /* Plot target CCD values */
+ {
+ double xx[128], y1[128], y2[128];
+
+ for (i = 0; i < 128; i++) {
+ xx[i] = i;
+ y1[i] = ccdsum[0][i];
+ y2[i] = ccdsum[1][i];
+ }
- for (k = 0; k < mtx_nocoef2[j]; k++)
- coeff2[j][k].we += err/mtx_nocoef2[j];
-// for (k = 0; k < mtx_nocoef2[j]; k++)
-// coeff2[j][k].we *= 1.0/tot;
+ printf("Raw target and actual CCD weight sums:\n");
+ do_plot(xx, y1, y2, NULL, 128);
+ }
+#endif
+
+ /* Figure valid range and extrapolate to edges */
+ dth[0] = 0.0; /* ref */
+ dth[1] = 0.007; /* hires */
+
+ for (k = 0; k < 2; k++) {
+
+ for (i = 0; i < 128; i++) {
+ if (ccdsum[k][i] > max[k])
+ max[k] = ccdsum[k][i];
}
- /* Check CCD sums */
- for (i = 0; i < nraw; i++)
- ccdsum[i] = 0.0;
+//printf("~1 max[%d] = %f\n",k, max[k]);
+ /* Figure out the valid range */
+ for (i = 64; i >= 0; i--) {
+ if (ccdsum[k][i] > (0.8 * max[k])) {
+ x[0] = (double)i;
+ } else {
+ break;
+ }
+ }
+ for (i = 64; i < 128; i++) {
+ if (ccdsum[k][i] > (0.8 * max[k])) {
+ x[3] = (double)i;
+ } else {
+ break;
+ }
+ }
+ /* Space off the last couple of entries */
+ x[0] += 2.0;
+ x[3] -= 6.0;
+ x[1] = floor((2 * x[0] + x[3])/3.0);
+ x[2] = floor((x[0] + 2 * x[3])/3.0);
+
+ for (i = 0; i < 4; i++)
+ y[i] = ccdsum[k][(int)x[i]];
+
+//printf("~1 extrap nodes %f, %f, %f, %f\n",x[0],x[1],x[2],x[3]);
+//printf("~1 extrap value %f, %f, %f, %f\n",y[0],y[1],y[2],y[3]);
- for (j = 0; j < m->nwav2; j++) {
- for (k = 0; k < mtx_nocoef2[j]; k++)
- ccdsum[coeff2[j][k].ix] += coeff2[j][k].we;
+ for (i = 0; i < 128; i++) {
+ double xw, yw;
+
+ xw = (double)i;
+
+ /* Compute interpolated value using Lagrange: */
+ yw = y[0] * (xw-x[1]) * (xw-x[2]) * (xw-x[3])
+ /((x[0]-x[1]) * (x[0]-x[2]) * (x[0]-x[3]))
+ + y[1] * (xw-x[0]) * (xw-x[2]) * (xw-x[3])
+ /((x[1]-x[0]) * (x[1]-x[2]) * (x[1]-x[3]))
+ + y[2] * (xw-x[0]) * (xw-x[1]) * (xw-x[3])
+ /((x[2]-x[0]) * (x[2]-x[1]) * (x[2]-x[3]))
+ + y[3] * (xw-x[0]) * (xw-x[1]) * (xw-x[2])
+ /((x[3]-x[0]) * (x[3]-x[1]) * (x[3]-x[2]));
+
+ if ((xw < x[0] || xw > x[3])
+ && fabs(ccdsum[k][i] - yw)/yw > dth[k]) {
+ ccdsum[k][i] = yw;
+ }
+ avg[k] += ccdsum[k][i];
}
+ avg[k] /= 128.0;
+ }
- if (ii >= 6)
- break;
+#ifdef HIGH_RES_PLOT
+ /* Plot target CCD values */
+ {
+ double xx[129], y1[129], y2[129];
+
+ for (i = 0; i < 128; i++) {
+ xx[i] = i;
+ y1[i] = ccdsum[0][i]/avg[0];
+ y2[i] = ccdsum[1][i]/avg[1];
+ }
+ xx[i] = i;
+ y1[i] = 0.0;
+ y2[i] = 0.0;
+
+ printf("Extrap. target and actual CCD weight sums:\n");
+ do_plot(xx, y1, y2, NULL, 129);
+ }
+#endif
+
+#ifdef DO_CCDNORMAVG /* Just correct by average */
+ for (cx = j = 0; j < m->nwav2; j++) { /* For each wavelength */
- /* Readjust CCD sum */
- for (i = 0; i < nraw; i++) {
- ccdsum[i] = ccdtsum[i]/ccdsum[i]; /* Amount to adjust */
+ /* For each matrix value */
+ sx = mtx_index2[j]; /* Starting index */
+ for (k = 0; k < mtx_nocoef2[j]; k++, cx++, sx++) {
+ mtx_coef2[cx] *= 10.0/twidth * avg[0]/avg[1];
}
+ }
+
+#else /* Correct by CCD bin */
+
+ /* Correct the weighting of each CCD value in the hires output */
+ for (i = 0; i < 128; i++) {
+ ccdsum[1][i] = 10.0/twidth * ccdsum[0][i]/ccdsum[1][i]; /* Correction factor */
+ }
+ for (cx = j = 0; j < m->nwav2; j++) { /* For each wavelength */
- for (j = 0; j < m->nwav2; j++) {
- for (k = 0; k < mtx_nocoef2[j]; k++)
- coeff2[j][k].we *= ccdsum[coeff2[j][k].ix];
+ /* For each matrix value */
+ sx = mtx_index2[j]; /* Starting index */
+ for (k = 0; k < mtx_nocoef2[j]; k++, cx++, sx++) {
+ mtx_coef2[cx] *= ccdsum[1][sx];
}
}
-#endif /* NEVER */
+#endif
}
+#endif /* DO_CCDNORM */
#ifdef HIGH_RES_PLOT
- /* Plot resampled curves */
{
+ static munki_fc coeff2[MXNOWL][MXNOFC];
double *xx, *ss;
double **yy;
+ /* Convert the native filter cooeficient representation to */
+ /* a 2D array we can randomly index. */
+ for (cx = j = 0; j < m->nwav2; j++) { /* For each output wavelength */
+ if (j >= MXNOWL) { /* Assert */
+ a1loge(p->log,1,"munki: number of hires output wavelenths is > %d\n",MXNOWL);
+ return MUNKI_INT_ASSERT;
+ }
+
+ /* For each matrix value */
+ sx = mtx_index2[j]; /* Starting index */
+ for (k = 0; k < mtx_nocoef2[j]; k++, cx++, sx++) {
+ if (k >= MXNOFC) { /* Assert */
+ a1loge(p->log,1,"munki: number of hires filter coeefs is > %d\n",MXNOFC);
+ return MUNKI_INT_ASSERT;
+ }
+ coeff2[j][k].ix = sx;
+ coeff2[j][k].we = mtx_coef2[cx];
+ }
+ }
+
xx = dvectorz(-1, m->nraw-1); /* X index */
yy = dmatrixz(0, 5, -1, m->nraw-1); /* Curves distributed amongst 5 graphs */
@@ -6675,36 +7230,14 @@ munki_code munki_create_hr(munki *p, int ref) {
}
}
- printf("Normalized Hi-Res wavelength sampling curves:\n");
+ printf("Normalized Hi-Res wavelength sampling curves: %s\n",ref ? "refl" : "emis");
do_plot6(xx, yy[0], yy[1], yy[2], yy[3], yy[4], yy[5], m->nraw);
free_dvector(xx, -1, m->nraw-1);
free_dmatrix(yy, 0, 2, -1, m->nraw-1);
}
#endif /* HIGH_RES_PLOT */
- /* Convert into runtime format */
- {
- int xcount;
-
- if ((*pmtx_index2 = mtx_index2 = (int *)calloc(m->nwav2, sizeof(int))) == NULL) {
- a1logd(p->log,1,"munki: malloc mtx_index2 failed!\n");
- return MUNKI_INT_MALLOC;
- }
-
- xcount = 0;
- for (j = 0; j < m->nwav2; j++) {
- mtx_index2[j] = coeff2[j][0].ix;
- xcount += mtx_nocoef2[j];
- }
-
- if ((*pmtx_coef2 = mtx_coef2 = (double *)calloc(xcount, sizeof(double))) == NULL) {
- a1logd(p->log,1,"munki: malloc mtx_coef2 failed!\n");
- return MUNKI_INT_MALLOC;
- }
-
- for (i = j = 0; j < m->nwav2; j++)
- for (k = 0; k < mtx_nocoef2[j]; k++, i++)
- mtx_coef2[i] = coeff2[j][k].we;
- }
+#undef MXNOWL
+#undef MXNOFC
/* Basic capability is initialised */
m->hr_inited++;
@@ -6713,9 +7246,9 @@ munki_code munki_create_hr(munki *p, int ref) {
/* If both reflective and emissive samplings have been created, */
/* deal with upsampling the references and calibrations */
if (m->hr_inited == 2) {
-#ifdef SALONEINSTLIB
+#ifdef ONEDSTRAYLIGHTUS
double **slp; /* 2D Array of stray light values */
-#endif /* !SALONEINSTLIB */
+#endif /* !ONEDSTRAYLIGHTUS */
rspl *trspl; /* Upsample rspl */
cow sd[40 * 40]; /* Scattered data points of existing references */
datai glow, ghigh;
@@ -6765,6 +7298,7 @@ munki_code munki_create_hr(munki *p, int ref) {
vhigh[0] = sd[i].v[0];
}
+#ifdef NEVER
/* Add some corrections at short wavelengths */
/* (The combination of the diffraction grating and */
/* LED light source doesn't give us much to work with here.) */
@@ -6788,13 +7322,14 @@ munki_code munki_create_hr(munki *p, int ref) {
sd[i].v[0] = 0.2 * sd[0].v[0];
sd[i++].w = 1.0;
}
+#endif
glow[0] = m->wl_short2;
ghigh[0] = m->wl_long2;
gres[0] = m->nwav2;
avgdev[0] = 0.0;
- trspl->fit_rspl_w(trspl, 0, sd, i, glow, ghigh, gres, vlow, vhigh, 0.5, avgdev, NULL);
+ trspl->fit_rspl_w(trspl, 0, sd, i, glow, ghigh, gres, vlow, vhigh, 5.0, avgdev, NULL);
if ((*ref2 = (double *)calloc(m->nwav2, sizeof(double))) == NULL) {
raw2wav->del(raw2wav);
@@ -6806,6 +7341,8 @@ munki_code munki_create_hr(munki *p, int ref) {
/* Create upsampled version */
for (i = 0; i < m->nwav2; i++) {
pp.p[0] = XSPECT_WL(m->wl_short2, m->wl_long2, m->nwav2, i);
+ if (pp.p[0] < min_wl) /* Duplicate below this wl */
+ pp.p[0] = min_wl;
trspl->interp(trspl, &pp);
if (pp.v[0] < 0.0)
pp.v[0] = 0.0;
@@ -6813,6 +7350,7 @@ munki_code munki_create_hr(munki *p, int ref) {
}
+#ifdef NEVER
/* Add some corrections at short wavelengths */
if (ii == 0) {
/* 376.67 - 470 */
@@ -6839,6 +7377,7 @@ munki_code munki_create_hr(munki *p, int ref) {
}
}
+#endif
#ifdef HIGH_RES_PLOT
/* Plot original and upsampled reference */
@@ -6885,7 +7424,7 @@ munki_code munki_create_hr(munki *p, int ref) {
}
trspl->del(trspl);
-#ifdef SALONEINSTLIB
+#ifdef ONEDSTRAYLIGHTUS
/* Then the 2D stray light using linear interpolation */
slp = dmatrix(0, m->nwav1-1, 0, m->nwav1-1);
@@ -6929,7 +7468,7 @@ munki_code munki_create_hr(munki *p, int ref) {
}
}
}
-#else /* !SALONEINSTLIB */
+#else /* !ONEDSTRAYLIGHTUS */
/* Then setup 2D stray light using rspl */
if ((trspl = new_rspl(RSPL_NOFLAGS, 2, 1)) == NULL) {
a1logd(p->log,3,"munki: creating rspl for high res conversion failed\n");
@@ -6961,7 +7500,7 @@ munki_code munki_create_hr(munki *p, int ref) {
avgdev[1] = 0.0;
trspl->fit_rspl_w(trspl, 0, sd, m->nwav1 * m->nwav1, glow, ghigh, gres, NULL, NULL, 0.5, avgdev, NULL);
-#endif /* !SALONEINSTLIB */
+#endif /* !ONEDSTRAYLIGHTUS */
m->straylight2 = dmatrixz(0, m->nwav2-1, 0, m->nwav2-1);
@@ -6970,8 +7509,12 @@ munki_code munki_create_hr(munki *p, int ref) {
for (j = 0; j < m->nwav2; j++) { /* Input wavelength */
double p0, p1;
p0 = XSPECT_WL(m->wl_short2, m->wl_long2, m->nwav2, i);
+ if (p0 < min_wl) /* Duplicate below this wl */
+ p0 = min_wl;
p1 = XSPECT_WL(m->wl_short2, m->wl_long2, m->nwav2, j);
-#ifdef SALONEINSTLIB
+ if (p1 < min_wl) /* Duplicate below this wl */
+ p1 = min_wl;
+#ifdef ONEDSTRAYLIGHTUS
/* Do linear interp with clipping at ends */
{
int x0, x1, y0, y1;
@@ -7002,11 +7545,11 @@ munki_code munki_create_hr(munki *p, int ref) {
+ w1 * v0 * slp[x1][y0]
+ w1 * v1 * slp[x1][y1];
}
-#else /* !SALONEINSTLIB */
+#else /* !ONEDSTRAYLIGHTUS */
pp.p[0] = p0;
pp.p[1] = p1;
trspl->interp(trspl, &pp);
-#endif /* !SALONEINSTLIB */
+#endif /* !ONEDSTRAYLIGHTUS */
m->straylight2[i][j] = pp.v[0] * HIGHRES_WIDTH/10.0;
if (m->straylight2[i][j] > 0.0)
m->straylight2[i][j] = 0.0;
@@ -7077,11 +7620,11 @@ munki_code munki_create_hr(munki *p, int ref) {
}
#endif /* HIGH_RES_PLOT */
-#ifdef SALONEINSTLIB
+#ifdef ONEDSTRAYLIGHTUS
free_dmatrix(slp, 0, m->nwav1-1, 0, m->nwav1-1);
-#else /* !SALONEINSTLIB */
+#else /* !ONEDSTRAYLIGHTUS */
trspl->del(trspl);
-#endif /* !SALONEINSTLIB */
+#endif /* !ONEDSTRAYLIGHTUS */
/* - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Allocate space for per mode calibration reference */
@@ -7278,22 +7821,9 @@ munki_code munki_conv2XYZ(
if (conv == NULL)
return MUNKI_INT_CIECONVFAIL;
- /* Don't report any wavelengths below the minimum for this mode */
- if ((s->min_wl-1e-3) > wl_short) {
- double wl;
- for (j = 0; j < m->nwav; j++) {
- wl = XSPECT_WL(m->wl_short, m->wl_long, m->nwav, j);
- if (wl >= (s->min_wl-1e-3))
- break;
- }
- six = j;
- wl_short = wl;
- nwl -= six;
- }
-
- a1logd(p->log,3,"munki_conv2XYZ got wl_short %f, wl_long %f, nwav %d, min_wl %f\n"
+ a1logd(p->log,3,"munki_conv2XYZ got wl_short %f, wl_long %f, nwav %d\n"
" after skip got wl_short %f, nwl = %d\n",
- m->wl_short, m->wl_long, m->nwav, s->min_wl, wl_short, nwl);
+ m->wl_short, m->wl_long, m->nwav, wl_short, nwl);
for (sms = 0.0, i = 1; i < 21; i++)
sms += opt_adj_weights[i];
@@ -7493,15 +8023,15 @@ munki_code munki_optimise_sensor(
new_int_time *= s->targoscale2;
a1logd(p->log,3,"Using compromse sensor target\n");
}
-#ifdef USE_HIGH_GAIN_MODE
- /* Hmm. It may not be a good idea to use high gain mode if it compromises */
+ /* Hmm. It seems not be a good idea to use high gain mode if it compromises */
/* the longer integration time which reduces noise. */
- if (new_int_time > m->max_int_time && permithg) {
- new_int_time /= m->highgain;
- new_gain_mode = 1;
- a1logd(p->log,3,"Switching to high gain mode\n");
+ if (s->auto_gain) {
+ if (new_int_time > m->max_int_time && permithg) {
+ new_int_time /= m->highgain;
+ new_gain_mode = 1;
+ a1logd(p->log,3,"Switching to high gain mode\n");
+ }
}
-#endif
}
a1logd(p->log,3,"after low light adjust, inttime %f, gain mode %d\n",new_int_time,new_gain_mode);
@@ -7616,10 +8146,8 @@ munki_interp_dark(
return MUNKI_INT_NOTCALIBRATED;
i = 0;
-#ifdef USE_HIGH_GAIN_MODE
- if (gainmode)
+ if (s->auto_gain && gainmode)
i = 2;
-#endif
for (j = -1; j < m->nraw; j++) {
double tt;
@@ -7654,7 +8182,7 @@ inst_opt_type munki_get_trig(munki *p) {
}
/* Switch thread handler */
-int munki_switch_thread(void *pp) {
+static int munki_switch_thread(void *pp) {
int nfailed = 0;
munki *p = (munki *)pp;
munkiimp *m = (munkiimp *)p->m;
@@ -7683,15 +8211,51 @@ int munki_switch_thread(void *pp) {
p->eventcallback(p->event_cntx, inst_event_switch);
}
} else if (ecode == mk_eve_spos_change) {
- if (p->eventcallback != NULL) {
- p->eventcallback(p->event_cntx, inst_event_mconf);
+#ifdef FILTER_SPOS_EVENTS
+ /* Signal change to filer thread */
+ m->spos_msec = msec_time();
+ m->spos_change++;
+#else
+ if (m->eventcallback != NULL) {
+ m->eventcallback(p->event_cntx, inst_event_mconf);
}
+#endif
}
}
a1logd(p->log,3,"Switch thread returning\n");
return rv;
}
+#ifdef FILTER_SPOS_EVENTS
+static int munki_spos_thread(void *pp) {
+ munki *p = (munki *)pp;
+ munkiimp *m = (munkiimp *)p->m;
+ int change = m->spos_change; /* Current count */
+
+ a1logd(p->log,3,"spos thread started\n");
+
+ for (;;) {
+
+ if (m->spos_th_term) {
+ m->spos_th_termed = 1;
+ break;
+ }
+
+ /* Do callback if change has persisted for 1 second */
+ if (change != m->spos_change
+ && (msec_time() - m->spos_msec) >= FILTER_TIME) {
+ change = m->spos_change;
+ if (p->eventcallback != NULL) {
+ p->eventcallback(p->event_cntx, inst_event_mconf);
+ }
+ }
+ msec_sleep(100);
+ }
+ return 0;
+}
+#endif
+
+
/* ============================================================ */
/* Low level commands */
@@ -8246,7 +8810,7 @@ munki_code munki_simulate_event(munki *p, mk_eve ecode, int timestamp) {
msec_sleep(50);
if (m->th_termed == 0) {
a1logd(p->log,1,"munki_simulate_event: terminate switch thread failed, canceling I/O\n");
- p->icom->usb_cancel_io(p->icom, &m->cancelt);
+ p->icom->usb_cancel_io(p->icom, &m->sw_cancel);
}
return rv;
@@ -8323,7 +8887,7 @@ munki_code munki_waitfor_switch_th(munki *p, mk_eve *ecode, int *timest, double
a1logd(p->log,2,"munki_waitfor_switch_th: Read 8 bytes from switch hit port\n");
/* Now read 8 bytes */
- se = p->icom->usb_read(p->icom, &m->cancelt, 0x83, buf, 8, &rwbytes, top);
+ se = p->icom->usb_read(p->icom, &m->sw_cancel, 0x83, buf, 8, &rwbytes, top);
if (se & ICOM_TO) {
a1logd(p->log,1,"munki_waitfor_switch_th: read 0x%x bytes, timed out\n",rwbytes);
@@ -8561,21 +9125,6 @@ munki_code munki_parse_eeprom(munki *p, unsigned char *buf, unsigned int len) {
if ((m->rmtx_coef1 = d->get_32_doubles(d, NULL, 184, 36 * 16)) == NULL)
return MUNKI_DATA_RANGE;
-#ifdef NEVER
-// ~~~ hack !!!
-m->rmtx_index1[4] = m->rmtx_index1[5] + 3;
-m->rmtx_index1[3] = m->rmtx_index1[4] + 3;
-m->rmtx_index1[2] = m->rmtx_index1[3] + 3;
-m->rmtx_index1[1] = m->rmtx_index1[2] + 3;
-m->rmtx_index1[0] = m->rmtx_index1[1] + 3;
-
-for(i = 0; i < 5; i++) {
- for (j = 0; j < 16; j++) {
- m->rmtx_coef1[i * 16 + j] = m->rmtx_coef1[5 * 16 + j];
- }
-}
-#endif
-
if (p->log->debug >= 7) {
a1logd(p->log,7,"Reflectance matrix:\n");
for(i = 0; i < 36; i++) {
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-25 10:51:09 +0000