diff options
Diffstat (limited to 'spectro/munki_imp.c')
-rw-r--r-- | spectro/munki_imp.c | 1263 |
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++) { |