Quantification of dual-state 5-ALA-induced PpIX fluorescence: methodology and validation in tissue-mimicking phantoms.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches

[1] § Materials & methods › Optical properties extraction ↔ src/iad_main.c, lines 52–156 · score 0.67 · reduced scattering coefficient, scattering anisotropy, absorption coefficient, IAD, optical properties, transmittance
[2] § Materials & methods › Optical properties extraction ↔ src/iad_main.c, lines 52–156 · score 0.52 · dual beam spectrophotometer, Optical properties, transmittance, wavelength, sphere, port

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

C · 3,065 lines · 118 KB · MIT · 2 matches

/* Autogenerated v4-0-0 from https://github.com/scottprahl/iad */
#define _CRT_SECURE_NO_WARNINGS
#define _CRT_NONSTDC_NO_WARNINGS
#define NO_SLIDES 0
#define ONE_SLIDE_ON_TOP 1
#define TWO_IDENTICAL_SLIDES 2
#define ONE_SLIDE_ON_BOTTOM 3
#define ONE_SLIDE_NEAR_SPHERE 4
#define ONE_SLIDE_NOT_NEAR_SPHERE 5
#define MR_IS_ONLY_RD 1
#define MT_IS_ONLY_TD 2
#define NO_UNSCATTERED_LIGHT 3
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
#include <time.h>
#include <math.h>
#include <ctype.h>
#include <errno.h>
#include "ad_globl.h"
#include "ad_prime.h"
#include "iad_type.h"
#include "iad_pub.h"
#include "iad_io.h"
#include "iad_calc.h"
#include "iad_util.h"
#include "version.h"
#include "mc_lost.h"
#include "ad_frsnl.h"
static void print_version(int verbosity)
{
if (verbosity == 0) {
fprintf(stdout, "%s", VersionShort);
}
else {
fprintf(stdout, "iad %s\n", Version);
fprintf(stdout, "Copyright 1993-2026 Scott Prahl, [email hidden]\n");
fprintf(stdout, " (see Applied Optics, 32:559-568, 1993)\n\n");
fprintf(stdout, "This is free software; see the source for copying conditions.\n");
fprintf(stdout, "There is no warranty; not even for MERCHANTABILITY or FITNESS.\n");
fprintf(stdout, "FOR A PARTICULAR PURPOSE.\n");
}
}
static void print_usage(void)
{
fprintf(stdout, "iad %s\n\n", Version);
fprintf(stdout, "iad finds optical properties from measurements\n\n");
fprintf(stdout, "Usage: iad [options] input\n\n");
fprintf(stdout, "Options:\n");
fprintf(stdout, " -1 '# # # # #' reflection sphere parameters \n");
fprintf(stdout, " 'd_sphere d d_sample_port d_entrance_port d_detector_port r_wall'\n");
fprintf(stdout, " -2 '# # # # #' transmission sphere parameters \n");
fprintf(stdout, " 'd_sphere d d_sample_port d_third_port d_detector_port r_wall'\n");
fprintf(stdout, " -a # use this albedo \n");
fprintf(stdout, " -A # use this absorption coefficient \n");
fprintf(stdout, " -b # use this optical thickness \n");
fprintf(stdout, " -B # beam diameter \n");
fprintf(stdout, " -c # fraction of unscattered refl in MR\n");
fprintf(stdout, " -C # fraction of unscattered trans in MT\n");
fprintf(stdout, " -d # thickness of sample \n");
fprintf(stdout, " -D # thickness of slide \n");
fprintf(stdout, " -e # error tolerance (default 0.0001) \n");
fprintf(stdout, " -E # optical depth (=mua*D) for slides\n");
fprintf(stdout, " -f # allow a fraction 0.0-1.0 of light to hit sphere wall first\n");
fprintf(stdout, " -F # constrain scattering coefficient \n");
fprintf(stdout, " # = constant: use constant scattering coefficient \n");
fprintf(stdout, " # = 'P lambda0 mus0 gamma' then mus=mus0*(lambda/lambda0)^gamma\n");
fprintf(stdout, " -g # scattering anisotropy (default 0) \n");
fprintf(stdout, " -G # type of boundary '0', '2', 't', 'b', 'n', 'f' \n");
fprintf(stdout, " '0' or '2' --- number of slides\n");
fprintf(stdout, " 't' (top) or 'b' (bottom) \
--- one slide that is hit by light first\n");
fprintf(stdout, " 'n' (near) or 'f' (far) \
--- one slide position relative to sphere\n");
fprintf(stdout, " -h display help\n");
fprintf(stdout, " -H # # = 0, no baffles for R or T spheres\n");
fprintf(stdout, " # = 1, baffle for R but not for T sphere\n");
fprintf(stdout, " # = 2, baffle for T but not for R sphere\n");
fprintf(stdout, " # = 3, baffle for both R and T spheres (default)\n");
fprintf(stdout, " -i # incident angle in degrees\n");
fprintf(stdout, " -j # constrain reduced scattering coefficient \n");
fprintf(stdout, " -J generate grid after inverse calculation\n");
fprintf(stdout, " -l # wavelength limits\n");
fprintf(stdout, " -L # specify the wavelength lambda\n");
fprintf(stdout, " -M # limit number of Monte Carlo iterations\n");
fprintf(stdout, " -n # specify index of refraction of slab\n");
fprintf(stdout, " -N # specify index of refraction of slides\n");
fprintf(stdout, " -o filename explicitly specify filename for output\n");
fprintf(stdout, " -p # # of Monte Carlo photons (default 100000)\n");
fprintf(stdout, " a negative number is max MC time in milliseconds\n");
fprintf(stdout, " -q # number of quadrature points (default=8)\n");
fprintf(stdout, " -r # total reflection measurement\n");
fprintf(stdout, " -R # actual reflectance for 100%% measurement \n");
fprintf(stdout, " -s # specify type of search to do\n");
fprintf(stdout, " -S # number of spheres used\n");
fprintf(stdout, " -t # total transmission measurement\n");
fprintf(stdout, " -T # actual transmission for 100%% measurement \n");
fprintf(stdout, " -u # unscattered transmission measurement\n");
fprintf(stdout, " -v version information\n");
fprintf(stdout, " -V 0 verbosity low --- no output to stdout\n");
fprintf(stdout, " -V 1 verbosity moderate \n");
fprintf(stdout, " -V 2 verbosity high\n");
fprintf(stdout, " -w # wall reflectivity for reflection sphere\n");
fprintf(stdout, " -W # wall reflectivity for transmission sphere\n");
fprintf(stdout, " -x # set debugging level\n");
fprintf(stdout, " -X dual beam configuration\n");
fprintf(stdout, " -Y ignore the diffuse lost-light correction (temporary)\n");
fprintf(stdout, " -z do forward calculation\n");
fprintf(stdout, "Examples:\n");
fprintf(stdout, " iad file.rxt Results will be put in file.txt\n");
fprintf(stdout, " iad file Same as above\n");
fprintf(stdout, " iad -c 0.9 file.rxt \
Assume M_R includes 90%% of unscattered reflectance\n");
fprintf(stdout, " iad -C 0.8 file.rxt \
Assume M_T includes 80%% of unscattered transmittance\n");
fprintf(stdout, " iad -e 0.0001 file.rxt Better convergence to R & T values\n");
fprintf(stdout, " iad -f 1.0 file.rxt All light hits reflectance sphere wall first\n");
fprintf(stdout, " iad -l '500 600' file.rxt Only do wavelengths between 500 and 600\n");
fprintf(stdout, " iad -o out file.rxt Calculated values in out\n");
fprintf(stdout, " iad -r 0.3 R_total=0.3, b=inf, find albedo\n");
fprintf(stdout, " iad -r 0.3 -t 0.4 R_total=0.3, T_total=0.4, find a,b,g\n");
fprintf(stdout, " iad -r 0.3 -t 0.4 -n 1.5 R_total=0.3, T_total=0.4, n=1.5, find a,b\n");
fprintf(stdout, " iad -r 0.3 -t 0.4 R_total=0.3, T_total=0.4, find a,b\n");
fprintf(stdout, " iad -p 1000 file.rxt Only 1000 photons\n");
fprintf(stdout, " iad -p -100 file.rxt Allow only 100ms per iteration\n");
fprintf(stdout, " iad -q 4 file.rxt Four quadrature points\n");
fprintf(stdout, " iad -M 0 file.rxt No MC (iad)\n");
fprintf(stdout, " iad -M 1 file.rxt MC once (iad -> MC -> iad)\n");
fprintf(stdout, " iad -M 2 file.rxt MC twice (iad -> MC -> iad -> MC -> iad)\n");
fprintf(stdout, " iad -M 0 -q 4 file.rxt Fast and crude conversion\n");
fprintf(stdout, " iad -G t file.rxt One top slide with properties from file.rxt\n");
fprintf(stdout, " iad -G b -N 1.5 -D 1 file Use 1 bottom slide with n=1.5 and thickness=1\n");
fprintf(stdout, " iad -x 1 file.rxt Show sphere and MC effects\n");
fprintf(stdout, " iad -x 2 file.rxt Show grid decisions\n");
fprintf(stdout, " iad -x 4 file.rxt Show iterations\n");
fprintf(stdout, " iad -x 8 file.rxt Show lost light effects\n");
fprintf(stdout, " iad -x 16 file.rxt Show best grid points\n");
fprintf(stdout, " iad -x 32 file.rxt Show decisions for type of search\n");
fprintf(stdout, " iad -x 64 file.rxt Show all grid calculations\n");
fprintf(stdout, " iad -x 128 file.rxt Show sphere calculations\n");
fprintf(stdout, " iad -x 256 file.rxt DEBUG_EVERY_CALC\n");
fprintf(stdout, " iad -x 511 file.rxt Show all debugging output\n");
fprintf(stdout, " iad -X -i 8 file.rxt Dual beam spectrometer with 8 degree incidence\n\n");
fprintf(stdout, " iad -z -a 0.9 -b 1 -i 45 Forward calc assuming 45 degree incidence\n\n");
fprintf(stdout, " iad * Process all .rxt files in current directory\n");
fprintf(stdout, " iad x.rxt y.rxt Process multiple files\n\n");
fprintf(stdout, "Report bugs to <[email hidden]>\n\n");
}
static char *strdup_together(char *s, char *t)
{
char *both;
if (s == NULL) {
if (t == NULL)
return NULL;
return strdup(t);
}
if (t == NULL)
return strdup(s);
both = malloc(strlen(s) + strlen(t) + 1);
if (both == NULL)
fprintf(stderr, "Could not allocate memory for both strings.\n");
strcpy(both, s);
strcat(both, t);
return both;
}
static double my_strtod(const char *str)
{
char *endptr;
errno = 0;
double val = strtod(str, &endptr);
if (endptr == str) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " No conversion could be performed for `%s`.\n", str);
exit(EXIT_FAILURE);
}
if (*endptr != '\0') {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " Partial conversion of string = '%s'\n", str);
exit(EXIT_FAILURE);
}
if (errno == ERANGE) {
fprintf(stderr, "Error in the command line\n");
printf(" The value '%s' is out of range of double.\n", str);
exit(EXIT_FAILURE);
}
return val;
}
static double seconds_elapsed(clock_t start_time)
{
clock_t finish_time = clock();
return (double) (finish_time - start_time) / CLOCKS_PER_SEC;
}
static void print_error_legend(void)
{
if (Debug(DEBUG_ANY))
return;
fprintf(stderr, "----------------- Sorry, but ... errors encountered ---------------\n");
fprintf(stderr, " * ==> Success ");
fprintf(stderr, " 0-9 ==> Monte Carlo Iteration\n");
fprintf(stderr, " R ==> M_R is too big ");
fprintf(stderr, " r ==> M_R is too small\n");
fprintf(stderr, " T ==> M_T is too big ");
fprintf(stderr, " t ==> M_T is too small\n");
fprintf(stderr, " U ==> M_U is too big ");
fprintf(stderr, " u ==> M_U is too small\n");
fprintf(stderr, " ! ==> M_R + M_T > 1 ");
fprintf(stderr, " + ==> Hit iteration limit\n");
fprintf(stderr, " x ==> No solution found");
fprintf(stderr, " m ==> Lost-light correction failed\n");
fprintf(stderr, " P ==> Port geometry is wrong");
fprintf(stderr, " L ==> Lost light puts data out of reach\n\n");
}
static char what_char(int err)
{
if (err == IAD_NO_ERROR)
return '*';
if (err == IAD_TOO_MANY_ITERATIONS)
return '+';
if (err == IAD_MR_TOO_BIG)
return 'R';
if (err == IAD_MR_TOO_SMALL)
return 'r';
if (err == IAD_MT_TOO_BIG)
return 'T';
if (err == IAD_MT_TOO_SMALL)
return 't';
if (err == IAD_MU_TOO_BIG)
return 'U';
if (err == IAD_MU_TOO_SMALL)
return 'u';
if (err == IAD_TOO_MUCH_LIGHT)
return '!';
if (err == IAD_SEARCH_STALLED)
return 'x';
if (err == IAD_MC_DID_NOT_CONVERGE)
return 'm';
if (err == IAD_UNREACHABLE_WITH_LOST_LIGHT)
return 'L';
if (err == IAD_BEAM_TOO_BIG_FOR_SAMPLE_PORT)
return 'P';
if (err == IAD_BEAM_TOO_BIG_FOR_ENTRANCE_PORT)
return 'P';
if (err == IAD_SAMPLE_PORTS_DIFFER)
return 'P';
if (err == IAD_BEAM_TOO_BIG_FOR_EXIT_PORT)
return 'P';
if (err == IAD_BEAM_NOT_VALID)
return 'P';
return '?';
}
static void print_long_error(int err)
{
switch (err) {
case IAD_NO_ERROR:
fprintf(stderr, "Successful Search\n");
break;
case IAD_TOO_MANY_ITERATIONS:
fprintf(stderr, "Failed Search, hit the limit of %d iterations\n", IAD_MAX_ITERATIONS);
break;
case IAD_SEARCH_STALLED:
fprintf(stderr, "Failed Search, stopped early without matching the measurements\n");
fprintf(stderr, " no combination of a, b, and g reproduced M_R and M_T;\n");
fprintf(stderr, " compare the measured and fitted values above\n");
break;
case IAD_BEAM_TOO_BIG_FOR_SAMPLE_PORT:
fprintf(stderr, "Failed Search, beam is wider than the sample port\n");
fprintf(stderr, " part of the beam missed the sample entirely;\n");
fprintf(stderr, " check the beam diameter and the sample port size\n");
break;
case IAD_BEAM_TOO_BIG_FOR_ENTRANCE_PORT:
fprintf(stderr, "Failed Search, beam is wider than the entrance port\n");
fprintf(stderr, " the beam was clipped entering the sphere;\n");
fprintf(stderr, " check the beam diameter and the entrance port size\n");
break;
case IAD_BEAM_NOT_VALID:
fprintf(stderr, "Failed Search, the beam has no width\n");
fprintf(stderr, " a beam diameter must be positive; the default is 1 mm\n");
fprintf(stderr, " and -B sets it explicitly\n");
break;
case IAD_BEAM_TOO_BIG_FOR_EXIT_PORT:
fprintf(stderr, "Failed Search, beam is wider than the exit port\n");
fprintf(stderr, " only part of the direct beam would reach the standard;\n");
fprintf(stderr, " use an exit port of zero to close it off entirely\n");
break;
case IAD_SAMPLE_PORTS_DIFFER:
fprintf(stderr, "Failed Search, the two spheres disagree about the sample port\n");
fprintf(stderr, " the sample sits in one hole, so both sphere blocks\n");
fprintf(stderr, " must give it the same diameter\n");
break;
case IAD_MC_DID_NOT_CONVERGE:
fprintf(stderr, "Failed Search, lost-light correction did not settle\n");
fprintf(stderr, " the Monte Carlo re-inversion stopped converging;\n");
fprintf(stderr, " try -M 0 to skip the lost-light correction\n");
break;
case IAD_UNREACHABLE_WITH_LOST_LIGHT:
fprintf(stderr, "Failed Search, lost light puts the data out of reach\n");
fprintf(stderr, " once the estimated lost light is subtracted, no sample\n");
fprintf(stderr, " reproduces these measurements: even one that absorbs\n");
fprintf(stderr, " nothing is too dim. Either the measurements or the\n");
fprintf(stderr, " sphere and port description do not describe the same\n");
fprintf(stderr, " experiment. Use -x 8 to see the shortfall, or -M 0 to\n");
fprintf(stderr, " invert without the lost-light correction\n");
break;
case IAD_MR_TOO_BIG:
fprintf(stderr, "Failed Search, M_R is too big\n");
break;
case IAD_MR_TOO_SMALL:
fprintf(stderr, "Failed Search, M_R is too small\n");
break;
case IAD_MT_TOO_BIG:
fprintf(stderr, "Failed Search, M_T is too big\n");
break;
case IAD_MT_TOO_SMALL:
fprintf(stderr, "Failed Search, M_T is too small\n");
break;
case IAD_MU_TOO_BIG:
fprintf(stderr, "Failed Search, M_U is too big\n");
break;
case IAD_MU_TOO_SMALL:
fprintf(stderr, "Failed Search, M_U is too small\n");
break;
case IAD_TOO_MUCH_LIGHT:
fprintf(stderr, "Failed Search, M_R + M_T exceeds 1\n");
break;
case IAD_RT_LT_MINIMUM:
fprintf(stderr, "Failed Search, M_R + M_T is below the minimum possible\n");
break;
case IAD_EXCESSIVE_LIGHT_LOSS:
fprintf(stderr, "Failed Search, too much light lost out the sides\n");
break;
case IAD_AS_NOT_VALID:
fprintf(stderr, "Failed Search, sample port is too big for the sphere\n");
break;
case IAD_AE_NOT_VALID:
fprintf(stderr, "Failed Search, entrance port is too big for the sphere\n");
break;
case IAD_AD_NOT_VALID:
fprintf(stderr, "Failed Search, detector port is too big for the sphere\n");
break;
case IAD_RW_NOT_VALID:
fprintf(stderr, "Failed Search, sphere wall reflectance is not between 0 and 1\n");
break;
case IAD_RD_NOT_VALID:
fprintf(stderr, "Failed Search, detector reflectance is not between 0 and 1\n");
break;
case IAD_RSTD_NOT_VALID:
fprintf(stderr, "Failed Search, reflectance standard is not between 0 and 1\n");
break;
case IAD_TSTD_NOT_VALID:
fprintf(stderr, "Failed Search, transmittance standard is not between 0 and 1\n");
break;
case IAD_QUAD_PTS_NOT_VALID:
fprintf(stderr, "Failed Search, number of quadrature points is not valid\n");
break;
case IAD_BAD_G_VALUE:
fprintf(stderr, "Failed Search, anisotropy is not between -1 and 1\n");
break;
case IAD_BAD_PHASE_FUNCTION:
fprintf(stderr, "Failed Search, unknown phase function\n");
break;
case IAD_GAMMA_NOT_VALID:
fprintf(stderr, "Failed Search, gamma is not valid\n");
break;
case IAD_F_NOT_VALID:
fprintf(stderr, "Failed Search, sphere wall fraction is not valid\n");
break;
case IAD_TOO_MANY_LAYERS:
fprintf(stderr, "Failed Search, too many layers\n");
break;
case IAD_MEMORY_ERROR:
fprintf(stderr, "Failed Search, out of memory\n");
break;
case IAD_FILE_ERROR:
fprintf(stderr, "Failed Search, error reading the input file\n");
break;
default:
fprintf(stderr, "Failed Search, error %d\n", err);
break;
}
fprintf(stderr, "\n");
}
static void print_dot(clock_t start_time, int err, int points, int final, int verbosity)
{
static int counter = 0;
counter++;
if (verbosity == 0 || Debug(DEBUG_ANY))
return;
if (final)
fprintf(stderr, "%c", what_char(err));
else {
counter--;
fprintf(stderr, "%1d\b", points % 10);
}
if (final) {
if (counter % 50 == 0) {
double rate = (seconds_elapsed(start_time) / counter);
fprintf(stderr, " %3d done (%5.2f s/pt)\n", counter, rate);
}
else if (counter % 10 == 0)
fprintf(stderr, " ");
}
fflush(stderr);
}
static void bright_mr_mt(struct measure_type m, struct invert_type r, double b, double *m_r, double *m_t)
{
r.slab.a = 1.0;
r.slab.b = b;
r.a = 1.0;
r.b = b;
Calculate_MR_MT(m, r, MC_USE_EXISTING, TRUE, m_r, m_t);
}
static double solve_for_b_from_mt(struct measure_type m, struct invert_type r, double target)
{
double lo = 1e-6, hi = 1e4, mid, m_r, m_t;
int i;
for (i = 0; i < 40; i++) {
mid = sqrt(lo * hi);
bright_mr_mt(m, r, mid, &m_r, &m_t);
if (m_t > target)
lo = mid;
else
hi = mid;
}
return sqrt(lo * hi);
}
static void calculate_coefficients(struct measure_type m,
struct invert_type r, double *LR, double *LT, double *musp, double *mua)
{
double mus;
*LR = 0;
*LT = 0;
if (r.error == IAD_NO_ERROR || r.error == IAD_TOO_MANY_ITERATIONS ||
r.error == IAD_SEARCH_STALLED || r.error == IAD_MC_DID_NOT_CONVERGE ||
r.error == IAD_UNREACHABLE_WITH_LOST_LIGHT) {
Calculate_MR_MT(m, r, MC_USE_EXISTING, TRUE, LR, LT);
Calculate_Mua_Musp(m, r, &mus, musp, mua);
}
else {
*musp = 0;
*mua = 0;
}
}
static int parse_string_into_array(char *s, double *a, int n)
{
char *t, *last, *r;
int i = 0;
t = s;
last = s + strlen(s);
while (t < last) {
r = t;
while (*r != ' ' && *r != '\0')
r++;
*r = '\0';
if (sscanf(t, "%lf", &(a[i])) == 0)
return 1;
i++;
if (i == n) {
if (i == 5) {
if (a[i - 1] <= 0 || a[i - 1] > 1) {
fprintf(stderr, "Sphere wall reflectivity (r_w=%g) must be a fraction less than one.\n", a[i - 1]);
exit(EXIT_FAILURE);
}
}
return 0;
}
t = r + 1;
}
return 1;
}
static int has_rxt_extension(char *s)
{
size_t len;
char *extension = ".rxt";
size_t extension_len = strlen(extension);
if (s == NULL)
return 0;
len = strlen(s);
if (len < extension_len)
return 0;
return strcmp(s + len - extension_len, extension) == 0;
}
static void print_results_header(FILE *fp)
{
if (Debug(DEBUG_LOST_LIGHT)) {
fprintf(fp, "# | Meas M_R | Meas M_T | calc calc calc |");
fprintf(fp, " Lost Lost Lost Lost | MC AD Error\n");
fprintf(fp, "# wave | M_R fit | M_T fit | mu_a mu_s' g | ");
fprintf(fp, " UR1 URU UT1 UTU | # # Type\n");
fprintf(fp, "# nm | --- --- | --- --- | 1/mm 1/mm --- |");
fprintf(fp, " --- --- --- --- | --- --- ---\n");
fprintf(fp, "#---------------------------------------------------------");
fprintf(fp, "--------------------------------------------------------\n");
}
else {
fprintf(fp, "# \t Meas \t Fit \t Meas \t Fit \t \t \t \t");
fprintf(fp, "\n");
fprintf(fp, "##wave\t M_R \t M_R \t M_T \t M_T \t mu_a \t mu_s'\t g \t");
fprintf(fp, "\n");
fprintf(fp, "# [nm]\t[---] \t[---] \t[---] \t[---] \t[1/mm]\t[1/mm]\t[---] \t");
fprintf(fp, "\n");
}
}
void print_optical_property_result(FILE *fp,
struct measure_type m, struct invert_type r, double LR, double LT, double mu_a, double mu_sp, int line)
{
int display_error = (!r.found && r.error == IAD_NO_ERROR) ? IAD_SEARCH_STALLED : r.error;
if (Debug(DEBUG_LOST_LIGHT)) {
if (m.lambda != 0)
fprintf(fp, "%6.1f ", m.lambda);
else
fprintf(fp, "%6d ", line);
if (mu_a >= 200)
mu_a = 199.9999;
if (mu_sp >= 1000)
mu_sp = 999.9999;
fprintf(fp, "%6.4f % 6.4f | ", m.m_r, LR);
fprintf(fp, "%6.4f % 6.4f | ", m.m_t, LT);
fprintf(fp, "%6.3f ", mu_a);
fprintf(fp, "%6.3f ", mu_sp);
fprintf(fp, "%6.3f |", r.g);
fprintf(fp, " %6.4f %6.4f ", m.lost_r.direct, m.lost_r.diffuse);
fprintf(fp, "%6.4f %6.4f | ", m.lost_t.direct, m.utu_lost);
fprintf(fp, "%2d ", r.MC_iterations);
fprintf(fp, "%3d", r.AD_iterations);
fprintf(fp, " %c \n", what_char(display_error));
}
else {
if (m.lambda != 0)
fprintf(fp, "%6.1f\t", m.lambda);
else
fprintf(fp, "%6d\t", line);
if (mu_a >= 200)
mu_a = 199.9999;
if (mu_sp >= 1000)
mu_sp = 999.9999;
fprintf(fp, "%6.4f\t%6.4f\t", m.m_r, LR);
fprintf(fp, "%6.4f\t%6.4f\t", m.m_t, LT);
fprintf(fp, "%6.4f\t", mu_a);
fprintf(fp, "%6.4f\t", mu_sp);
fprintf(fp, "%6.4f\t", r.g);
fprintf(fp, " %c \n", what_char(display_error));
}
fflush(fp);
}
int main(int argc, char **argv)
{
struct measure_type m;
struct invert_type r;
char *g_out_name = NULL;
char *g_grid_name = NULL;
int c;
long n_photons = 100000;
int MAX_MC_iterations = 19;
int any_error = 0;
int last_error = IAD_NO_ERROR;
int process_command_line = 0;
int params = 0;
int rt_total = 0;
int mc_total = 0;
int file_index = 0;
int file_count = 1;
int cl_quadrature_points = UNINITIALIZED;
int cl_verbosity = 2;
double cl_forward_calc = UNINITIALIZED;
double cl_grid_calc = UNINITIALIZED;
double cl_default_a = UNINITIALIZED;
double cl_default_g = UNINITIALIZED;
double cl_default_b = UNINITIALIZED;
double cl_default_mua = UNINITIALIZED;
double cl_default_mus = UNINITIALIZED;
double cl_default_musp = UNINITIALIZED;
double cl_tolerance = UNINITIALIZED;
double cl_slide_OD = UNINITIALIZED;
double cl_cos_angle = UNINITIALIZED;
double cl_beam_d = UNINITIALIZED;
double cl_sample_d = UNINITIALIZED;
double cl_sample_n = UNINITIALIZED;
double cl_slide_d = UNINITIALIZED;
double cl_slide_n = UNINITIALIZED;
double cl_slides = UNINITIALIZED;
double cl_default_fr = UNINITIALIZED;
double cl_rstd_t = UNINITIALIZED;
double cl_rstd_r = UNINITIALIZED;
double cl_baffle_r = UNINITIALIZED;
double cl_baffle_t = UNINITIALIZED;
double cl_ru_fraction = UNINITIALIZED;
double cl_tu_fraction = UNINITIALIZED;
double cl_lambda = UNINITIALIZED;
double cl_rwall_r = UNINITIALIZED;
double cl_rwall_t = UNINITIALIZED;
double cl_search = UNINITIALIZED;
double cl_mus0 = UNINITIALIZED;
double cl_mus0_pwr = UNINITIALIZED;
double cl_mus0_lambda = UNINITIALIZED;
double cl_UR1 = UNINITIALIZED;
double cl_UT1 = UNINITIALIZED;
double cl_Tc = UNINITIALIZED;
double cl_method = UNINITIALIZED;
int cl_num_spheres = UNINITIALIZED;
double cl_sphere_one[5] = { UNINITIALIZED, UNINITIALIZED, UNINITIALIZED,
UNINITIALIZED, UNINITIALIZED
};
double cl_sphere_two[5] = { UNINITIALIZED, UNINITIALIZED, UNINITIALIZED,
UNINITIALIZED, UNINITIALIZED
};
double cl_wave_limit[2] = { UNINITIALIZED, UNINITIALIZED };
clock_t start_time = clock();
char command_line_options[] = "1:2:a:A:b:B:c:C:d:D:e:E:f:F:g:G:hH:i:j:Jl:L:M:n:N:o:p:q:r:R:s:S:t:T:u:vV:w:W:x:XYz";
char *command_line = NULL;
{
size_t command_line_length = 0;
for (int i = 0; i < argc; ++i) {
command_line_length += strlen(argv[i]) + 3;
}
command_line = (char *) malloc(command_line_length);
if (command_line == NULL) {
fprintf(stderr, "Memory allocation failed\n");
return 1;
}
strcpy(command_line, "");
for (int i = 0; i < argc; ++i) {
if (strchr(argv[i], ' ') != NULL) {
strcat(command_line, "'");
strcat(command_line, argv[i]);
strcat(command_line, "' ");
}
else {
strcat(command_line, argv[i]);
strcat(command_line, " ");
}
}
optind = 1;
}
while ((c = getopt(argc, argv, command_line_options)) != EOF) {
int n;
char cc;
char *tmp_str = NULL;
switch (c) {
case '1':
tmp_str = strdup(optarg);
parse_string_into_array(optarg, cl_sphere_one, 5);
if (cl_sphere_one[4] == UNINITIALIZED) {
fprintf(stderr, "Error in the command-line argument for -1\n");
fprintf(stderr, " the current argument is '%s' but it must have 5 terms: ", tmp_str);
fprintf(stderr, "'d_sphere d_sample d_entrance d_detector r_wall'\n");
exit(EXIT_FAILURE);
}
break;
case '2':
tmp_str = strdup(optarg);
parse_string_into_array(optarg, cl_sphere_two, 5);
if (cl_sphere_two[4] == UNINITIALIZED) {
fprintf(stderr, "Error in the command-line argument for -2\n");
fprintf(stderr, " the current argument is '%s' but it must have 5 terms: ", tmp_str);
fprintf(stderr, "'d_sphere d_sample d_third d_detector r_wall'\n");
exit(EXIT_FAILURE);
}
break;
case 'a':
cl_default_a = my_strtod(optarg);
if (cl_default_a < 0 || cl_default_a > 1) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " albedo '-a %s'\n", optarg);
exit(EXIT_FAILURE);
}
break;
case 'A':
cl_default_mua = my_strtod(optarg);
if (cl_default_mua < 0) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " absorption '-A %s'\n", optarg);
exit(EXIT_FAILURE);
}
break;
case 'b':
cl_default_b = my_strtod(optarg);
if (cl_default_b < 0) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " optical thickness '-b %s'\n", optarg);
exit(EXIT_FAILURE);
}
break;
case 'B':
cl_beam_d = my_strtod(optarg);
if (cl_beam_d < 0) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " beam diameter '-B %s'\n", optarg);
exit(EXIT_FAILURE);
}
break;
case 'c':
cl_ru_fraction = my_strtod(optarg);
if (cl_ru_fraction < 0.0 || cl_ru_fraction > 1.0) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " unscattered refl fraction '-c %s'\n", optarg);
fprintf(stderr, " must be between 0 and 1\n");
exit(EXIT_SUCCESS);
}
break;
case 'C':
cl_tu_fraction = my_strtod(optarg);
if (cl_tu_fraction < 0.0 || cl_tu_fraction > 1.0) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " unscattered trans fraction '-C %s'\n", optarg);
fprintf(stderr, " must be between 0 and 1\n");
exit(EXIT_SUCCESS);
}
break;
case 'd':
cl_sample_d = my_strtod(optarg);
if (cl_sample_d < 0) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " sample thickness '-d %s'\n", optarg);
exit(EXIT_FAILURE);
}
break;
case 'D':
cl_slide_d = my_strtod(optarg);
if (cl_slide_d < 0) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " slide thickness '-D %s'\n", optarg);
exit(EXIT_FAILURE);
}
break;
case 'e':
cl_tolerance = my_strtod(optarg);
if (cl_tolerance < 0) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " error tolerance '-e %s'\n", optarg);
exit(EXIT_FAILURE);
}
break;
case 'E':
cl_slide_OD = my_strtod(optarg);
if (cl_slide_OD < 0) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " slide optical depth '-E %s'\n", optarg);
exit(EXIT_FAILURE);
}
break;
case 'f':
cl_default_fr = my_strtod(optarg);
if (cl_default_fr < 0.0 || cl_default_fr > 1.0) {
fprintf(stderr, "Error in the command-line argument: ");
fprintf(stderr, "'-f %s' The argument must be between 0 and 1.\n", optarg);
exit(EXIT_SUCCESS);
}
break;
case 'F':
if (isdigit(optarg[0])) {
cl_default_mus = my_strtod(optarg);
if (cl_default_mus < 0) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " mus '-F %s'\n", optarg);
exit(EXIT_FAILURE);
}
break;
}
n = sscanf(optarg, "%c %lf %lf %lf", &cc, &cl_mus0_lambda, &cl_mus0, &cl_mus0_pwr);
if (n != 4 || (cc != 'P' && cc != 'p')) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " bad -F option. '-F %s'\n", optarg);
fprintf(stderr, " -F 1.0 for mus =1.0\n");
fprintf(stderr, " -F 'P 500 1.0 -1.3' for mus =1.0*(lambda/500)^(-1.3)\n");
exit(EXIT_FAILURE);
}
break;
case 'g':
cl_default_g = my_strtod(optarg);
if (cl_default_g < -1 || cl_default_g > 1) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " anisotropy '-g %s'\n", optarg);
exit(EXIT_FAILURE);
}
if (cl_default_g == -1)
cl_default_g = -MAX_ABS_G;
if (cl_default_g == 1)
cl_default_g = MAX_ABS_G;
break;
case 'G':
if (optarg[0] == '0')
cl_slides = NO_SLIDES;
else if (optarg[0] == '2')
cl_slides = TWO_IDENTICAL_SLIDES;
else if (optarg[0] == 't' || optarg[0] == 'T')
cl_slides = ONE_SLIDE_ON_TOP;
else if (optarg[0] == 'b' || optarg[0] == 'B')
cl_slides = ONE_SLIDE_ON_BOTTOM;
else if (optarg[0] == 'n' || optarg[0] == 'N')
cl_slides = ONE_SLIDE_NEAR_SPHERE;
else if (optarg[0] == 'f' || optarg[0] == 'F')
cl_slides = ONE_SLIDE_NOT_NEAR_SPHERE;
else {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " Argument for '-G %s' must be \n", optarg);
fprintf(stderr, " 't' --- light always hits top slide first\n");
fprintf(stderr, " 'b' --- light always hits bottom slide first\n");
fprintf(stderr, " 'n' --- slide always closest to sphere\n");
fprintf(stderr, " 'f' --- slide always farthest from sphere\n");
exit(EXIT_FAILURE);
}
break;
case 'H':
if (optarg[0] == '0') {
cl_baffle_r = 0;
cl_baffle_t = 0;
}
else if (optarg[0] == '1') {
cl_baffle_r = 1;
cl_baffle_t = 0;
}
else if (optarg[0] == '2') {
cl_baffle_r = 0;
cl_baffle_t = 1;
}
else if (optarg[0] == '3') {
cl_baffle_r = 1;
cl_baffle_t = 1;
}
else {
fprintf(stderr, "Error in the command-line -H argument\n");
fprintf(stderr, " argument is '%s', but ", optarg);
fprintf(stderr, "must be 0, 1, 2, or 3\n");
exit(EXIT_FAILURE);
}
break;
case 'i':
cl_cos_angle = my_strtod(optarg);
if (cl_cos_angle < 0 || cl_cos_angle > 90) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " incident angle '-i %s'\n", optarg);
fprintf(stderr, " must be between 0 and 90 degrees\n");
exit(EXIT_FAILURE);
}
cl_cos_angle = cos(cl_cos_angle * M_PI / 180.0);
break;
case 'j':
cl_default_musp = my_strtod(optarg);
if (cl_default_musp < 0) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " musp '-j %s'\n", optarg);
exit(EXIT_FAILURE);
}
break;
case 'J':
cl_grid_calc = 1;
break;
case 'l':
tmp_str = strdup(optarg);
parse_string_into_array(optarg, cl_wave_limit, 2);
break;
case 'L':
cl_lambda = my_strtod(optarg);
break;
case 'M':
MAX_MC_iterations = (int) my_strtod(optarg);
if (MAX_MC_iterations < 0 || MAX_MC_iterations > 50) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " MC iterations '-M %s'\n", optarg);
exit(EXIT_FAILURE);
}
break;
case 'n':
cl_sample_n = my_strtod(optarg);
if (cl_sample_n < 0.1 || cl_sample_n > 10) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " slab index '-n %s'\n", optarg);
exit(EXIT_FAILURE);
}
break;
case 'N':
cl_slide_n = my_strtod(optarg);
if (cl_slide_n < 0.1 || cl_slide_n > 10) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " slide index '-N %s'\n", optarg);
exit(EXIT_FAILURE);
}
break;
case 'o':
g_out_name = strdup(optarg);
break;
case 'p':
n_photons = (long) my_strtod(optarg);
break;
case 'q':
cl_quadrature_points = (int) my_strtod(optarg);
if (cl_quadrature_points % 4 != 0) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " '-q %s'\n", optarg);
fprintf(stderr, " Quadrature points must be a multiple of 4\n");
exit(EXIT_FAILURE);
}
if ((cl_cos_angle != UNINITIALIZED) && (cl_quadrature_points % 12 != 0)) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " '-q %s'\n", optarg);
fprintf(stderr, " Quadrature points must be multiple of 12 for oblique incidence\n");
exit(EXIT_FAILURE);
}
break;
case 'r':
cl_UR1 = my_strtod(optarg);
process_command_line = 1;
if (cl_UR1 < 0 || cl_UR1 > 1) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " UR1 value '-r %s'\n", optarg);
fprintf(stderr, " must be between 0 and 1\n");
exit(EXIT_FAILURE);
}
break;
case 'R':
cl_rstd_r = my_strtod(optarg);
if (cl_rstd_r < 0 || cl_rstd_r > 1) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " Rstd value '-R %s'\n", optarg);
fprintf(stderr, " must be between 0 and 1\n");
exit(EXIT_FAILURE);
}
break;
case 's':
cl_search = (int) my_strtod(optarg);
break;
case 'S':
cl_num_spheres = (int) my_strtod(optarg);
if (cl_num_spheres != 0 && cl_num_spheres != 1 && cl_num_spheres != 2) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " sphere number '-S %s'\n", optarg);
fprintf(stderr, " must be 0, 1, or 2\n");
exit(EXIT_FAILURE);
}
break;
case 't':
cl_UT1 = my_strtod(optarg);
if (cl_UT1 < 0 || cl_UT1 > 1) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " UT1 value '-t %s'\n", optarg);
fprintf(stderr, " must be between 0 and 1\n");
exit(EXIT_FAILURE);
}
process_command_line = 1;
break;
case 'T':
cl_rstd_t = my_strtod(optarg);
if (cl_rstd_t < 0 || cl_rstd_t > 1) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " transmission standard '-T %s'\n", optarg);
fprintf(stderr, " must be between 0 and 1\n");
exit(EXIT_FAILURE);
}
break;
case 'u':
cl_Tc = my_strtod(optarg);
if (cl_Tc < 0 || cl_Tc > 1) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " unscattered transmission '-u %s'\n", optarg);
fprintf(stderr, " must be between 0 and 1\n");
exit(EXIT_FAILURE);
}
process_command_line = 1;
break;
case 'v':
print_version(cl_verbosity);
exit(EXIT_SUCCESS);
break;
case 'V':
cl_verbosity = my_strtod(optarg);
break;
case 'w':
cl_rwall_r = my_strtod(optarg);
if (cl_rwall_r < 0 || cl_rwall_r > 1) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " refl sphere wall '-w %s'\n", optarg);
fprintf(stderr, " must be between 0 and 1\n");
exit(EXIT_FAILURE);
}
break;
case 'W':
cl_rwall_t = my_strtod(optarg);
if (cl_rwall_t < 0 || cl_rwall_r > 1) {
fprintf(stderr, "Error in the command line\n");
fprintf(stderr, " trans sphere wall '-w %s'\n", optarg);
fprintf(stderr, " must be between 0 and 1\n");
exit(EXIT_FAILURE);
}
break;
case 'x':
Set_Debugging((int) my_strtod(optarg));
break;
case 'X':
cl_method = COMPARISON;
break;
case 'Y':
MC_Include_Diffuse_Loss(0);
break;
case 'z':
cl_forward_calc = 1;
process_command_line = 1;
break;
default:
fprintf(stderr, "unknown option '%c'\n", c);
case 'h':
print_usage();
exit(EXIT_SUCCESS);
}
}
argc -= optind;
argv += optind;
Initialize_Measure(&m);
if (cl_cos_angle != UNINITIALIZED) {
m.slab_cos_angle = cl_cos_angle;
if (cl_quadrature_points == UNINITIALIZED)
cl_quadrature_points = 12;
if (cl_quadrature_points != 12 * (cl_quadrature_points / 12)) {
fprintf(stderr, "If you use the -i option to specify an oblique incidence angle, then\n");
fprintf(stderr, "the number of quadrature points must be a multiple of 12\n");
exit(EXIT_SUCCESS);
}
}
if (cl_sample_n != UNINITIALIZED)
m.slab_index = cl_sample_n;
if (cl_slide_n != UNINITIALIZED) {
m.slab_bottom_slide_index = cl_slide_n;
m.slab_top_slide_index = cl_slide_n;
if (cl_slide_d == UNINITIALIZED && !Column_Label_Present('D')) {
if (m.slab_top_slide_thickness == 0)
m.slab_top_slide_thickness = 1.0;
if (m.slab_bottom_slide_thickness == 0)
m.slab_bottom_slide_thickness = 1.0;
}
}
if (cl_slide_OD != UNINITIALIZED) {
m.slab_bottom_slide_b = cl_slide_OD;
m.slab_top_slide_b = cl_slide_OD;
}
if (cl_sample_d != UNINITIALIZED)
m.slab_thickness = cl_sample_d;
if (cl_beam_d != UNINITIALIZED)
m.d_beam = cl_beam_d;
if (cl_slide_d != UNINITIALIZED) {
m.slab_bottom_slide_thickness = cl_slide_d;
m.slab_top_slide_thickness = cl_slide_d;
}
if (cl_slides == NO_SLIDES) {
m.slab_bottom_slide_index = 1.0;
m.slab_bottom_slide_thickness = 0.0;
m.slab_top_slide_index = 1.0;
m.slab_top_slide_thickness = 0.0;
}
if (cl_slides == ONE_SLIDE_ON_TOP || cl_slides == ONE_SLIDE_NEAR_SPHERE) {
m.slab_bottom_slide_index = 1.0;
m.slab_bottom_slide_thickness = 0.0;
}
if (cl_slides == ONE_SLIDE_ON_BOTTOM || cl_slides == ONE_SLIDE_NOT_NEAR_SPHERE) {
m.slab_top_slide_index = 1.0;
m.slab_top_slide_thickness = 0.0;
}
if (cl_slides == ONE_SLIDE_NEAR_SPHERE || cl_slides == ONE_SLIDE_NOT_NEAR_SPHERE)
m.flip_sample = 1;
else
m.flip_sample = 0;
if (cl_slides == NO_SLIDES) {
m.slab_top_slide_b = 0.0;
m.slab_bottom_slide_b = 0.0;
}
if (cl_slides == ONE_SLIDE_ON_TOP || cl_slides == ONE_SLIDE_NEAR_SPHERE)
m.slab_bottom_slide_b = 0.0;
if (cl_slides == ONE_SLIDE_ON_BOTTOM || cl_slides == ONE_SLIDE_NOT_NEAR_SPHERE)
m.slab_top_slide_b = 0.0;
if (cl_method != UNINITIALIZED)
m.method = (int) cl_method;
if (cl_rstd_r != UNINITIALIZED) {
m.rstd_r = cl_rstd_r;
m.rstd_t = cl_rstd_r;
}
if (cl_rstd_t != UNINITIALIZED) {
m.rstd_t = cl_rstd_t;
if (cl_rstd_r == UNINITIALIZED)
m.rstd_r = cl_rstd_t;
}
if (cl_rwall_r != UNINITIALIZED || cl_rwall_t != UNINITIALIZED) {
if (cl_sphere_one[0] != UNINITIALIZED || cl_sphere_two[0] != UNINITIALIZED) {
fprintf(stderr, "A wall reflectance cannot accompany a sphere description.\n");
fprintf(stderr, " -1 and -2 already carry one as their fifth value\n");
fprintf(stderr, " -w and -W are for overriding the wall in an .rxt file\n");
exit(EXIT_FAILURE);
}
}
if (cl_rwall_r != UNINITIALIZED)
m.rw_r = cl_rwall_r;
if (cl_rwall_t != UNINITIALIZED)
m.rw_t = cl_rwall_t;
if (cl_sphere_one[0] != UNINITIALIZED) {
double d_sample_r, d_third_r, d_detector_r;
m.d_sphere_r = cl_sphere_one[0];
d_sample_r = cl_sphere_one[1];
d_third_r = cl_sphere_one[2];
d_detector_r = cl_sphere_one[3];
m.rw_r = cl_sphere_one[4];
m.as_r = sqr(d_sample_r / m.d_sphere_r / 2);
m.at_r = sqr(d_third_r / m.d_sphere_r / 2);
m.ad_r = sqr(d_detector_r / m.d_sphere_r / 2);
m.aw_r = 1.0 - m.as_r - m.at_r - m.ad_r;
m.d_sphere_t = m.d_sphere_r;
m.as_t = m.as_r;
m.at_t = m.at_r;
m.ad_t = m.ad_r;
m.aw_t = m.aw_r;
m.rw_t = m.rw_r;
if (cl_num_spheres == UNINITIALIZED)
m.num_spheres = 1;
}
if (cl_sphere_two[0] != UNINITIALIZED) {
double d_sample_t, d_third_t, d_detector_t;
m.d_sphere_t = cl_sphere_two[0];
d_sample_t = cl_sphere_two[1];
d_third_t = cl_sphere_two[2];
d_detector_t = cl_sphere_two[3];
m.rw_t = cl_sphere_two[4];
m.as_t = sqr(d_sample_t / m.d_sphere_t / 2);
m.at_t = sqr(d_third_t / m.d_sphere_t / 2);
m.ad_t = sqr(d_detector_t / m.d_sphere_t / 2);
m.aw_t = 1.0 - m.as_t - m.at_t - m.ad_t;
if (cl_num_spheres == UNINITIALIZED)
m.num_spheres = 2;
}
if (cl_num_spheres != UNINITIALIZED) {
m.num_spheres = (int) cl_num_spheres;
if (m.num_spheres > 0 && m.method == UNKNOWN)
m.method = SUBSTITUTION;
}
if (cl_ru_fraction != UNINITIALIZED)
m.fraction_of_ru_in_mr = cl_ru_fraction;
if (cl_tu_fraction != UNINITIALIZED)
m.fraction_of_tu_in_mt = cl_tu_fraction;
if (cl_UR1 != UNINITIALIZED)
m.m_r = cl_UR1;
if (cl_UT1 != UNINITIALIZED)
m.m_t = cl_UT1;
if (cl_Tc != UNINITIALIZED)
m.m_u = cl_Tc;
if (cl_default_fr != UNINITIALIZED)
m.f_r = cl_default_fr;
if (cl_baffle_r != UNINITIALIZED)
m.baffle_r = cl_baffle_r;
if (cl_baffle_t != UNINITIALIZED)
m.baffle_t = cl_baffle_t;
if (cl_lambda != UNINITIALIZED)
m.lambda = cl_lambda;
Initialize_Result(m, &r, TRUE);
if (cl_forward_calc != UNINITIALIZED) {
if (cl_quadrature_points != UNINITIALIZED)
r.method.quad_pts = cl_quadrature_points;
else
r.method.quad_pts = 8;
if (cl_default_a != UNINITIALIZED)
r.default_a = cl_default_a;
if (cl_default_mua != UNINITIALIZED) {
r.default_mua = cl_default_mua;
if (cl_sample_d != UNINITIALIZED)
r.default_ba = cl_default_mua * cl_sample_d;
else
r.default_ba = cl_default_mua * m.slab_thickness;
}
if (cl_default_b != UNINITIALIZED)
r.default_b = cl_default_b;
if (cl_default_g != UNINITIALIZED)
r.default_g = cl_default_g;
if (cl_tolerance != UNINITIALIZED) {
r.tolerance = cl_tolerance;
r.MC_tolerance = cl_tolerance;
}
if (cl_mus0 != UNINITIALIZED) {
if (m.lambda != 0) {
cl_default_mus = cl_mus0 * pow(m.lambda / cl_mus0_lambda, cl_mus0_pwr);
}
else {
fprintf(stderr, "Seems like you want to constrain scattering to a power law.\n");
fprintf(stderr, "Unfortunately, there is no wavelength so this cannot be done.\n");
}
}
if (cl_default_mus != UNINITIALIZED) {
r.default_mus = cl_default_mus;
if (cl_sample_d != UNINITIALIZED)
r.default_bs = cl_default_mus * cl_sample_d;
else
r.default_bs = cl_default_mus * m.slab_thickness;
}
if (cl_default_musp != UNINITIALIZED) {
if (cl_default_g != UNINITIALIZED)
r.default_mus = cl_default_musp / (1.0 - cl_default_g);
else
r.default_mus = cl_default_musp;
if (cl_sample_d != UNINITIALIZED)
r.default_bs = r.default_mus * cl_sample_d;
else
r.default_bs = r.default_mus * m.slab_thickness;
}
if (cl_search != UNINITIALIZED)
r.search = cl_search;
double temp_mus = 1;
double temp_mua = 0;
r.g = 0;
if (cl_default_mus != UNINITIALIZED)
temp_mus = cl_default_mus;
if (cl_default_mua != UNINITIALIZED)
temp_mua = cl_default_mua;
if (cl_default_g != UNINITIALIZED)
r.g = cl_default_g;
if (cl_default_musp != UNINITIALIZED)
temp_mus = cl_default_musp / (1 - r.g);
if (cl_default_a != UNINITIALIZED) {
r.a = cl_default_a;
if (cl_default_b != UNINITIALIZED && cl_sample_d != UNINITIALIZED) {
temp_mus = cl_default_a * cl_default_b / cl_sample_d;
temp_mua = cl_default_b / cl_sample_d - temp_mus;
}
else {
if (cl_default_a == 0) {
temp_mus = 0;
temp_mua = 1;
}
else
temp_mua = temp_mus / cl_default_a - temp_mus;
}
}
else
r.a = temp_mus / (temp_mus + temp_mua);
if (cl_default_b != UNINITIALIZED) {
r.b = cl_default_b;
}
else {
if (cl_sample_d == UNINITIALIZED)
r.b = HUGE_VAL;
else
r.b = (temp_mus + temp_mua) * cl_sample_d;
}
r.slab.a = r.a;
r.slab.b = r.b;
r.slab.g = r.g;
{
double mu_s, mu_sp, mu_a, m_r, m_t;
double ur1, ut1, uru, utu, ru, tu;
double cos_critical, theta_inc, theta_crit;
double thickness, denom;
if (MAX_MC_iterations == 0 || m.num_spheres == 0) {
Calculate_MR_MT(m, r, MC_NONE, TRUE, &m_r, &m_t);
}
else {
Calculate_MR_MT(m, r, MC_REDO, TRUE, &m_r, &m_t);
}
Calculate_Mua_Musp(m, r, &mu_s, &mu_sp, &mu_a);
RT(r.method.quad_pts, &r.slab, &ur1, &ut1, &uru, &utu);
Sp_mu_RT_Flip(m.flip_sample,
r.slab.n_top_slide, r.slab.n_slab, r.slab.n_bottom_slide,
r.slab.b_top_slide, r.slab.b, r.slab.b_bottom_slide, r.slab.cos_angle, &ru, &tu);
thickness = (m.slab_thickness > 0) ? m.slab_thickness : 1.0;
cos_critical = Cos_Critical_Angle(r.slab.n_slab, 1.0);
theta_inc = acos(r.slab.cos_angle) * 180.0 / M_PI;
theta_crit = acos(cos_critical) * 180.0 / M_PI;
if (cl_verbosity > 0) {
printf("Intrinsic Properties\n");
printf(" albedo = %.3f\n", r.slab.a);
if (r.slab.b == HUGE_VAL)
printf(" optical thickness = inf\n");
else
printf(" optical thickness = %.3f\n", r.slab.b);
printf(" anisotropy = %.3f\n", r.slab.g);
printf(" thickness = %.3f mm\n", m.slab_thickness);
printf(" sample index = %.3f\n", r.slab.n_slab);
printf(" top slide index = %.3f\n", r.slab.n_top_slide);
printf(" bottom slide index = %.3f\n", r.slab.n_bottom_slide);
printf(" cos(theta incident) = %.3f\n", r.slab.cos_angle);
printf(" quadrature points = %d\n", r.method.quad_pts);
printf("\n");
printf("Derived quantities\n");
denom = (r.slab.b == HUGE_VAL) ? 1.0 : thickness;
printf(" mu_a = %.4f 1/mm\n",
(r.slab.b == HUGE_VAL) ? ((r.slab.a > 0) ? (1.0 - r.slab.a) / r.slab.a : 1.0)
: (1.0 - r.slab.a) * r.slab.b / denom);
printf(" mu_s = %.4f 1/mm\n",
(r.slab.b == HUGE_VAL) ? 1.0 : r.slab.a * r.slab.b / denom);
printf(" mu_s*(1-g) = %.4f 1/mm\n", (r.slab.b == HUGE_VAL) ? (1.0 - r.slab.g)
: (1.0 - r.slab.g) * r.slab.a * r.slab.b / denom);
printf(" theta incident = %.1f\xc2\xb0\n", theta_inc);
printf(" cos(theta critical) = %.4f\n", cos_critical);
printf(" theta critical = %.1f\xc2\xb0\n", theta_crit);
if (m.num_spheres > 0) {
printf("\n");
printf("\nSphere properties (%d sphere%s)\n", m.num_spheres, (m.num_spheres == 1) ? "" : "s");
{
const char *baffle_text = m.baffle_r ? "has a baffle" : "has no baffle";
printf(" Reflection sphere %s between sample and detector\n", baffle_text);
printf(" sphere diameter = %7.1f mm\n", m.d_sphere_r);
printf(" sample port diameter = %7.1f mm\n", 2 * m.d_sphere_r * sqrt(m.as_r));
printf(" entrance port diameter = %7.1f mm\n", 2 * m.d_sphere_r * sqrt(m.at_r));
printf(" detector port diameter = %7.1f mm\n", 2 * m.d_sphere_r * sqrt(m.ad_r));
printf(" detector reflectance = %7.1f %%\n", m.rd_r * 100);
printf(" wall reflectance = %7.1f %%\n", m.rw_r * 100);
printf(" calibration standard = %7.1f %%\n", m.rstd_r * 100);
}
if (m.num_spheres == 2) {
{
const char *baffle_text = m.baffle_t ? "has a baffle" : "has no baffle";
printf(" Transmission sphere %s between sample and detector\n", baffle_text);
printf(" sphere diameter = %7.1f mm\n", m.d_sphere_t);
printf(" sample port diameter = %7.1f mm\n",
2 * m.d_sphere_t * sqrt(m.as_t));
printf(" third port diameter = %7.1f mm\n",
2 * m.d_sphere_t * sqrt(m.at_t));
printf(" detector port diameter = %7.1f mm\n",
2 * m.d_sphere_t * sqrt(m.ad_t));
printf(" detector reflectance = %7.1f %%\n", m.rd_t * 100);
printf(" wall reflectance = %7.1f %%\n", m.rw_t * 100);
printf(" calibration standard = %7.1f %%\n", m.rstd_t * 100);
}
}
}
printf("Calculated quantities\n");
printf(" R total = %.4f\n", ur1);
printf(" R scattered = %.4f\n", ur1 - ru);
printf(" R unscattered = %.4f\n", ru);
printf(" T total = %.4f\n", ut1);
printf(" T scattered = %.4f\n", ut1 - tu);
printf(" T unscattered = %.4f\n", tu);
if (m.num_spheres > 0) {
printf(" M_R (sphere) = %.4f\n", m_r);
printf(" M_T (sphere) = %.4f\n", m_t);
}
}
(void) mu_s;
(void) mu_sp;
(void) mu_a;
}
exit(EXIT_SUCCESS);
}
{
int i;
for (i = 0; i < argc - 1; i++) {
if (strcmp(argv[i], "-o") == 0 && g_out_name == NULL) {
g_out_name = strdup(argv[i + 1]);
memmove(argv + i, argv + i + 2, (argc - i - 2) * sizeof(char *));
argc -= 2;
i--;
}
}
}
if (argc > 1) {
int file_arg_index;
int write_index = 0;
for (file_arg_index = 0; file_arg_index < argc; file_arg_index++) {
if (has_rxt_extension(argv[file_arg_index]))
argv[write_index++] = argv[file_arg_index];
}
argc = write_index;
if (argc == 0)
exit(EXIT_SUCCESS);
if (argc > 1 && g_out_name != NULL) {
fprintf(stderr, "The -o option cannot be used with multiple input files\n");
exit(EXIT_FAILURE);
}
}
if (argc == 1 && strcmp(argv[0], "-") != 0) {
int n;
char *base_name, *rt_name;
base_name = strdup(argv[0]);
n = (int) (strlen(base_name) - strlen(".rxt"));
if (n > 0 && strstr(base_name + n, ".rxt") != NULL)
base_name[n] = '\0';
rt_name = strdup_together(base_name, ".rxt");
if (freopen(argv[0], "r", stdin) == NULL && freopen(rt_name, "r", stdin) == NULL) {
fprintf(stderr, "Could not open either '%s' or '%s'\n", argv[0], rt_name);
exit(EXIT_FAILURE);
}
if (g_out_name == NULL)
g_out_name = strdup_together(base_name, ".txt");
if (g_grid_name == NULL)
g_grid_name = strdup_together(base_name, ".grid");
free(rt_name);
free(base_name);
process_command_line = 0;
}
if (g_out_name != NULL) {
if (freopen(g_out_name, "w", stdout) == NULL) {
fprintf(stderr, "Could not open file '%s' for output\n", g_out_name);
exit(EXIT_FAILURE);
}
}
if (process_command_line) {
if (cl_num_spheres >= 1 && cl_sphere_one[0] == UNINITIALIZED) {
fprintf(stderr, "Sphere measurements need the sphere described.\n");
fprintf(stderr, " -S %d was given without -1\n", cl_num_spheres);
fprintf(stderr, " -1 'd_sphere d_sample d_entrance d_detector r_wall'\n");
exit(EXIT_FAILURE);
}
if (cl_num_spheres == 2 && cl_sphere_two[0] == UNINITIALIZED) {
fprintf(stderr, "Two spheres need both spheres described.\n");
fprintf(stderr, " -S 2 was given without -2\n");
fprintf(stderr, " -2 'd_sphere d_sample d_third d_detector r_wall'\n");
exit(EXIT_FAILURE);
}
m.num_measures = 3;
if (m.m_r == 0)
m.num_measures--;
if (m.m_t == 0)
m.num_measures--;
if (m.m_u == 0)
m.num_measures--;
params = m.num_measures;
{
double ur1 = 0;
double ut1 = 0;
double uru = 0;
double utu = 0;
double mu_a = 0;
double mu_sp = 0;
double LR = 0;
double LT = 0;
int skip = FALSE;
if (cl_wave_limit[0] != UNINITIALIZED) {
if (m.lambda != 0) {
if (m.lambda < cl_wave_limit[0])
skip = TRUE;
if (m.lambda > cl_wave_limit[1])
skip = TRUE;
}
}
if (Debug(DEBUG_ANY) && !skip) {
fprintf(stderr, "\n-------------------NEXT DATA POINT---------------------\n");
if (m.lambda != 0)
fprintf(stderr, "lambda=%6.1f ", m.lambda);
fprintf(stderr, "MR=%8.5f MT=%8.5f\n\n", m.m_r, m.m_t);
if (skip)
fprintf(stderr, "skipping, wavelength out of range.\n");
}
if (!skip) {
rt_total++;
Initialize_Result(m, &r, FALSE);
if (cl_quadrature_points != UNINITIALIZED)
r.method.quad_pts = cl_quadrature_points;
else
r.method.quad_pts = 8;
if (cl_default_a != UNINITIALIZED)
r.default_a = cl_default_a;
if (cl_default_mua != UNINITIALIZED) {
r.default_mua = cl_default_mua;
if (cl_sample_d != UNINITIALIZED)
r.default_ba = cl_default_mua * cl_sample_d;
else
r.default_ba = cl_default_mua * m.slab_thickness;
}
if (cl_default_b != UNINITIALIZED)
r.default_b = cl_default_b;
if (cl_default_g != UNINITIALIZED)
r.default_g = cl_default_g;
if (cl_tolerance != UNINITIALIZED) {
r.tolerance = cl_tolerance;
r.MC_tolerance = cl_tolerance;
}
if (cl_mus0 != UNINITIALIZED) {
if (m.lambda != 0) {
cl_default_mus = cl_mus0 * pow(m.lambda / cl_mus0_lambda, cl_mus0_pwr);
}
else {
fprintf(stderr, "Seems like you want to constrain scattering to a power law.\n");
fprintf(stderr, "Unfortunately, there is no wavelength so this cannot be done.\n");
}
}
if (cl_default_mus != UNINITIALIZED) {
r.default_mus = cl_default_mus;
if (cl_sample_d != UNINITIALIZED)
r.default_bs = cl_default_mus * cl_sample_d;
else
r.default_bs = cl_default_mus * m.slab_thickness;
}
if (cl_default_musp != UNINITIALIZED) {
if (cl_default_g != UNINITIALIZED)
r.default_mus = cl_default_musp / (1.0 - cl_default_g);
else
r.default_mus = cl_default_musp;
if (cl_sample_d != UNINITIALIZED)
r.default_bs = r.default_mus * cl_sample_d;
else
r.default_bs = r.default_mus * m.slab_thickness;
}
if (cl_search != UNINITIALIZED)
r.search = cl_search;
if (cl_method == COMPARISON && m.d_sphere_r != 0 && m.as_r == 0) {
fprintf(stderr, "A dual-beam measurement is specified, but no port sizes.\n");
fprintf(stderr, "You might forsake the -X option and use zero spheres (which gives\n");
fprintf(stderr, "the same result except lost light is not taken into account).\n");
fprintf(stderr, "Alternatively, bite the bullet and enter your sphere parameters,\n");
fprintf(stderr, "with the knowledge that only the beam diameter and sample port\n");
fprintf(stderr, "diameter will be used to estimate lost light from the edges.\n");
exit(EXIT_SUCCESS);
}
if (cl_method == COMPARISON && m.num_spheres == 2) {
fprintf(stderr, "A dual-beam measurement is specified, but a two sphere experiment\n");
fprintf(stderr, "is specified. Since this seems impossible, I will make it\n");
fprintf(stderr, "impossible for you unless you specify 0 or 1 sphere.\n");
exit(EXIT_SUCCESS);
}
if (cl_method == COMPARISON && m.f_r != 0) {
fprintf(stderr, "A dual-beam measurement is specified, but a fraction of light\n");
fprintf(stderr, "is specified to hit the sphere wall first. This situation\n");
fprintf(stderr, "is not supported by iad. Sorry.\n");
exit(EXIT_SUCCESS);
}
if (rt_total == 1 && cl_verbosity > 0) {
Write_Header(m, r, params, command_line);
if (MAX_MC_iterations > 0) {
if (n_photons >= 0)
fprintf(stdout, "# Photons used to estimate lost light = %ld\n", n_photons);
else
fprintf(stdout, "# Time used to estimate lost light = %ld ms\n", -n_photons);
}
else
fprintf(stdout, "# Photons used to estimate lost light = 0\n");
fprintf(stdout, "#\n");
print_results_header(stdout);
}
m.lost_r.direct = 0;
m.lost_r.diffuse = 0;
m.lost_t.diffuse = 0;
m.lost_t.direct = 0;
m.utu_lost = 0;
Inverse_RT(m, &r);
calculate_coefficients(m, r, &LR, &LT, &mu_sp, &mu_a);
if (m.num_spheres > 0 && r.found && r.error == IAD_NO_ERROR) {
if (Debug(DEBUG_LOST_LIGHT)) {
print_results_header(stderr);
print_optical_property_result(stderr, m, r, LR, LT, mu_a, mu_sp, rt_total);
}
{
int mc_failed = 0;
int mc_unreachable = 0;
int pinned = 0;
struct measure_type good_m = m;
struct invert_type good_r = r;
double mc_prev_a = r.slab.a;
double mc_prev_b = r.slab.b;
double mc_prev_g = r.slab.g;
{
double b0 = (r.slab.b < 1e8) ? r.slab.b : 1.0;
r.mc_simplex_a_step = 1e-3;
r.mc_simplex_b_step = 1e-2 * (b0 > 1.0 ? b0 : 1.0);
r.mc_simplex_g_step = 1e-3;
}
int has_prev_diff_ur1_lost = 0;
double prev_abs_diff_ur1_lost = 0.0;
while (r.MC_iterations < MAX_MC_iterations) {
long n_photons_this;
double last_mu_sp, last_mu_a;
struct lost_type current_r, current_t;
double current_utu_lost;
double diff_ur1_lost, diff_ut1_lost, diff_uru_lost, diff_utu_lost;
double diff_uru_lost_t;
double factor = 0.3;
int too_much_lost;
double tol = r.MC_tolerance;
calculate_coefficients(m, r, &LR, &LT, &mu_sp, &mu_a);
last_mu_sp = mu_sp;
last_mu_a = mu_a;
if (Debug(DEBUG_ITERATIONS) || Debug(DEBUG_A_LITTLE)) {
fprintf(stderr, "\n------------- Monte Carlo Iteration %d -----------------\n",
r.MC_iterations + 1);
}
if (n_photons < 0)
n_photons_this = n_photons;
else if (!has_prev_diff_ur1_lost || prev_abs_diff_ur1_lost > 0.01)
n_photons_this = (n_photons / 10 > 10000) ? n_photons / 10 : 10000;
else if (prev_abs_diff_ur1_lost > 0.001)
n_photons_this = n_photons;
else
n_photons_this = (n_photons * 5 < 10000000) ? n_photons * 5 : 10000000;
MC_Lost(m, r, n_photons_this, &ur1, &ut1, &uru, &utu,
&current_r, &current_t, &current_utu_lost);
diff_ur1_lost = current_r.direct - m.lost_r.direct;
diff_uru_lost = current_r.diffuse - m.lost_r.diffuse;
diff_uru_lost_t = current_t.diffuse - m.lost_t.diffuse;
diff_ut1_lost = current_t.direct - m.lost_t.direct;
diff_utu_lost = current_utu_lost - m.utu_lost;
prev_abs_diff_ur1_lost = fabs(diff_ur1_lost);
has_prev_diff_ur1_lost = 1;
if (fabs(diff_ur1_lost) > 0.001 || fabs(diff_ut1_lost) > 0.001)
too_much_lost = 1;
else
too_much_lost = 0;
m.lost_r.direct += factor * diff_ur1_lost;
m.lost_r.diffuse += factor * diff_uru_lost;
m.lost_t.diffuse += factor * diff_uru_lost_t;
m.lost_t.direct += factor * diff_ut1_lost;
m.utu_lost += factor * diff_utu_lost;
mc_total++;
r.MC_iterations++;
Inverse_RT(m, &r);
{
double new_a = r.slab.a, new_b = r.slab.b, new_g = r.slab.g;
double da = fabs(new_a - mc_prev_a);
double db = fabs(new_b - mc_prev_b);
double dg = fabs(new_g - mc_prev_g);
double b_safe = (new_b < 1e8) ? new_b : 1.0;
double a_fixed = 1e-3;
double b_fixed = 1e-2 * (b_safe > 1.0 ? b_safe : 1.0);
double g_fixed = 1e-3;
r.mc_simplex_a_step = (da < 1e-5) ? 1e-5 : (da > a_fixed ? a_fixed : da);
r.mc_simplex_b_step = (db < 1e-4) ? 1e-4 : (db > b_fixed ? b_fixed : db);
r.mc_simplex_g_step = (dg < 1e-5) ? 1e-5 : (dg > g_fixed ? g_fixed : dg);
mc_prev_a = new_a;
mc_prev_b = new_b;
mc_prev_g = new_g;
}
calculate_coefficients(m, r, &LR, &LT, &mu_sp, &mu_a);
if (r.found && r.error == IAD_NO_ERROR) {
good_m = m;
good_r = r;
}
if (0) {
fprintf(stderr, "%2d %2d %2d | %7.4f %7.4f %7.4f | %7.4f %7.4f %7.4f\n",
r.MC_iterations, too_much_lost, r.found,
m.m_r, current_r.direct, m.lost_r.direct, m.m_t, current_t.direct, m.lost_t.direct);
}
if (Debug(DEBUG_LOST_LIGHT))
print_optical_property_result(stderr, m, r, LR, LT, mu_a, mu_sp, rt_total);
else
print_dot(start_time, r.error, mc_total, FALSE, cl_verbosity);
if (mu_a <= 0.0 && too_much_lost)
pinned++;
else
pinned = 0;
if (pinned >= 2) {
{
double b_thin = 1e-6;
double b_t, mr_at, mt_at, mr_thin, mt_thin, mr_dn, mt_dn, mr_up, mt_up;
int unreachable = 0;
bright_mr_mt(m, r, b_thin, &mr_thin, &mt_thin);
if (m.m_t > mt_thin) {
{
unreachable = 1;
b_t = b_thin;
mr_at = mr_thin;
mt_at = mt_thin;
}
}
else {
b_t = solve_for_b_from_mt(m, r, m.m_t);
bright_mr_mt(m, r, b_t, &mr_at, &mt_at);
bright_mr_mt(m, r, b_t * 0.9, &mr_dn, &mt_dn);
bright_mr_mt(m, r, b_t * 1.1, &mr_up, &mt_up);
if (m.m_r > mr_at && mr_up > mr_at && mr_at > mr_dn)
unreachable = 1;
}
if (unreachable) {
r.error = IAD_UNREACHABLE_WITH_LOST_LIGHT;
if (Debug(DEBUG_LOST_LIGHT) || Debug(DEBUG_A_LITTLE) || Debug(DEBUG_ITERATIONS)) {
fprintf(stderr, "The lost light puts the measurements out of reach.\n");
fprintf(stderr, " with no absorption at b = %.4f the model gives\n",
b_t);
fprintf(stderr, " M_R %8.5f measured %8.5f short by %8.5f\n", mr_at,
m.m_r, m.m_r - mr_at);
fprintf(stderr, " M_T %8.5f measured %8.5f\n", mt_at, m.m_t);
}
}
}
if (r.error == IAD_UNREACHABLE_WITH_LOST_LIGHT) {
mc_unreachable = 1;
if (Debug(DEBUG_ITERATIONS))
fprintf(stderr,
"absorption is spent and the data is out of reach — stopping\n");
break;
}
}
if (r.found) {
if (fabs(last_mu_a - mu_a) > tol) {
if (Debug(DEBUG_ITERATIONS))
fprintf(stderr, "Repeat MC because mua is still changing\n");
continue;
}
if (fabs(last_mu_sp - mu_sp) > tol) {
if (Debug(DEBUG_ITERATIONS))
fprintf(stderr, "Repeat MC because musp is still changing\n");
continue;
}
if (too_much_lost) {
if (Debug(DEBUG_ITERATIONS))
fprintf(stderr, "Repeat MC because mua and musp are still changing\n");
continue;
}
if (Debug(DEBUG_ITERATIONS))
fprintf(stderr, "found!\n");
break;
}
else {
if (Debug(DEBUG_ITERATIONS))
fprintf(stderr, "AD did not converge — stopping MC loop\n");
mc_failed = 1;
break;
}
}
if (mc_unreachable) {
r.found = 0;
r.error = IAD_UNREACHABLE_WITH_LOST_LIGHT;
{
int why = r.error;
struct lost_type lost_r = m.lost_r;
struct lost_type lost_t = m.lost_t;
double utu_lost = m.utu_lost;
m = good_m;
m.lost_r = lost_r;
m.lost_t = lost_t;
m.utu_lost = utu_lost;
r = good_r;
r.error = why;
r.found = 0;
}
}
else if (mc_failed) {
r.found = 0;
{
double b_thin = 1e-6;
double b_t, mr_at, mt_at, mr_thin, mt_thin, mr_dn, mt_dn, mr_up, mt_up;
int unreachable = 0;
bright_mr_mt(m, r, b_thin, &mr_thin, &mt_thin);
if (m.m_t > mt_thin) {
{
unreachable = 1;
b_t = b_thin;
mr_at = mr_thin;
mt_at = mt_thin;
}
}
else {
b_t = solve_for_b_from_mt(m, r, m.m_t);
bright_mr_mt(m, r, b_t, &mr_at, &mt_at);
bright_mr_mt(m, r, b_t * 0.9, &mr_dn, &mt_dn);
bright_mr_mt(m, r, b_t * 1.1, &mr_up, &mt_up);
if (m.m_r > mr_at && mr_up > mr_at && mr_at > mr_dn)
unreachable = 1;
}
if (unreachable) {
r.error = IAD_UNREACHABLE_WITH_LOST_LIGHT;
if (Debug(DEBUG_LOST_LIGHT) || Debug(DEBUG_A_LITTLE) || Debug(DEBUG_ITERATIONS)) {
fprintf(stderr, "The lost light puts the measurements out of reach.\n");
fprintf(stderr, " with no absorption at b = %.4f the model gives\n", b_t);
fprintf(stderr, " M_R %8.5f measured %8.5f short by %8.5f\n",
mr_at, m.m_r, m.m_r - mr_at);
fprintf(stderr, " M_T %8.5f measured %8.5f\n", mt_at, m.m_t);
}
}
}
if (r.error == IAD_NO_ERROR)
r.error = IAD_MC_DID_NOT_CONVERGE;
{
int why = r.error;
struct lost_type lost_r = m.lost_r;
struct lost_type lost_t = m.lost_t;
double utu_lost = m.utu_lost;
m = good_m;
m.lost_r = lost_r;
m.lost_t = lost_t;
m.utu_lost = utu_lost;
r = good_r;
r.error = why;
r.found = 0;
}
}
}
}
if (!r.found && r.error == IAD_NO_ERROR)
r.error = IAD_SEARCH_STALLED;
calculate_coefficients(m, r, &LR, &LT, &mu_sp, &mu_a);
print_optical_property_result(stdout, m, r, LR, LT, mu_a, mu_sp, rt_total);
if (r.error != IAD_NO_ERROR) {
any_error = 1;
last_error = r.error;
}
if (Debug(DEBUG_ANY))
print_long_error(r.error);
else
print_dot(start_time, r.error, mc_total, TRUE, cl_verbosity);
}
}
if (cl_verbosity > 0)
fprintf(stderr, "\n\n");
if (any_error && cl_verbosity > 1)
print_long_error(last_error);
exit(EXIT_SUCCESS);
}
file_count = argc > 1 ? argc : 1;
file_index = 0;
read_next_file:
if (argc > 1) {
{
int n;
char *base_name, *rt_name, *out_name;
base_name = strdup(argv[file_index]);
n = (int) (strlen(base_name) - strlen(".rxt"));
base_name[n] = '\0';
rt_name = strdup_together(base_name, ".rxt");
out_name = strdup_together(base_name, ".txt");
if (freopen(rt_name, "r", stdin) == NULL) {
fprintf(stderr, "Could not open file '%s'\n", rt_name);
exit(EXIT_FAILURE);
}
if (freopen(out_name, "w", stdout) == NULL) {
fprintf(stderr, "Could not open file '%s' for output\n", out_name);
exit(EXIT_FAILURE);
}
free(out_name);
free(rt_name);
free(base_name);
}
}
Initialize_Measure(&m);
if (cl_cos_angle != UNINITIALIZED) {
m.slab_cos_angle = cl_cos_angle;
if (cl_quadrature_points == UNINITIALIZED)
cl_quadrature_points = 12;
if (cl_quadrature_points != 12 * (cl_quadrature_points / 12)) {
fprintf(stderr, "If you use the -i option to specify an oblique incidence angle, then\n");
fprintf(stderr, "the number of quadrature points must be a multiple of 12\n");
exit(EXIT_SUCCESS);
}
}
if (cl_sample_n != UNINITIALIZED)
m.slab_index = cl_sample_n;
if (cl_slide_n != UNINITIALIZED) {
m.slab_bottom_slide_index = cl_slide_n;
m.slab_top_slide_index = cl_slide_n;
if (cl_slide_d == UNINITIALIZED && !Column_Label_Present('D')) {
if (m.slab_top_slide_thickness == 0)
m.slab_top_slide_thickness = 1.0;
if (m.slab_bottom_slide_thickness == 0)
m.slab_bottom_slide_thickness = 1.0;
}
}
if (cl_slide_OD != UNINITIALIZED) {
m.slab_bottom_slide_b = cl_slide_OD;
m.slab_top_slide_b = cl_slide_OD;
}
if (cl_sample_d != UNINITIALIZED)
m.slab_thickness = cl_sample_d;
if (cl_beam_d != UNINITIALIZED)
m.d_beam = cl_beam_d;
if (cl_slide_d != UNINITIALIZED) {
m.slab_bottom_slide_thickness = cl_slide_d;
m.slab_top_slide_thickness = cl_slide_d;
}
if (cl_slides == NO_SLIDES) {
m.slab_bottom_slide_index = 1.0;
m.slab_bottom_slide_thickness = 0.0;
m.slab_top_slide_index = 1.0;
m.slab_top_slide_thickness = 0.0;
}
if (cl_slides == ONE_SLIDE_ON_TOP || cl_slides == ONE_SLIDE_NEAR_SPHERE) {
m.slab_bottom_slide_index = 1.0;
m.slab_bottom_slide_thickness = 0.0;
}
if (cl_slides == ONE_SLIDE_ON_BOTTOM || cl_slides == ONE_SLIDE_NOT_NEAR_SPHERE) {
m.slab_top_slide_index = 1.0;
m.slab_top_slide_thickness = 0.0;
}
if (cl_slides == ONE_SLIDE_NEAR_SPHERE || cl_slides == ONE_SLIDE_NOT_NEAR_SPHERE)
m.flip_sample = 1;
else
m.flip_sample = 0;
if (cl_slides == NO_SLIDES) {
m.slab_top_slide_b = 0.0;
m.slab_bottom_slide_b = 0.0;
}
if (cl_slides == ONE_SLIDE_ON_TOP || cl_slides == ONE_SLIDE_NEAR_SPHERE)
m.slab_bottom_slide_b = 0.0;
if (cl_slides == ONE_SLIDE_ON_BOTTOM || cl_slides == ONE_SLIDE_NOT_NEAR_SPHERE)
m.slab_top_slide_b = 0.0;
if (cl_method != UNINITIALIZED)
m.method = (int) cl_method;
if (cl_rstd_r != UNINITIALIZED) {
m.rstd_r = cl_rstd_r;
m.rstd_t = cl_rstd_r;
}
if (cl_rstd_t != UNINITIALIZED) {
m.rstd_t = cl_rstd_t;
if (cl_rstd_r == UNINITIALIZED)
m.rstd_r = cl_rstd_t;
}
if (cl_rwall_r != UNINITIALIZED || cl_rwall_t != UNINITIALIZED) {
if (cl_sphere_one[0] != UNINITIALIZED || cl_sphere_two[0] != UNINITIALIZED) {
fprintf(stderr, "A wall reflectance cannot accompany a sphere description.\n");
fprintf(stderr, " -1 and -2 already carry one as their fifth value\n");
fprintf(stderr, " -w and -W are for overriding the wall in an .rxt file\n");
exit(EXIT_FAILURE);
}
}
if (cl_rwall_r != UNINITIALIZED)
m.rw_r = cl_rwall_r;
if (cl_rwall_t != UNINITIALIZED)
m.rw_t = cl_rwall_t;
if (cl_sphere_one[0] != UNINITIALIZED) {
double d_sample_r, d_third_r, d_detector_r;
m.d_sphere_r = cl_sphere_one[0];
d_sample_r = cl_sphere_one[1];
d_third_r = cl_sphere_one[2];
d_detector_r = cl_sphere_one[3];
m.rw_r = cl_sphere_one[4];
m.as_r = sqr(d_sample_r / m.d_sphere_r / 2);
m.at_r = sqr(d_third_r / m.d_sphere_r / 2);
m.ad_r = sqr(d_detector_r / m.d_sphere_r / 2);
m.aw_r = 1.0 - m.as_r - m.at_r - m.ad_r;
m.d_sphere_t = m.d_sphere_r;
m.as_t = m.as_r;
m.at_t = m.at_r;
m.ad_t = m.ad_r;
m.aw_t = m.aw_r;
m.rw_t = m.rw_r;
if (cl_num_spheres == UNINITIALIZED)
m.num_spheres = 1;
}
if (cl_sphere_two[0] != UNINITIALIZED) {
double d_sample_t, d_third_t, d_detector_t;
m.d_sphere_t = cl_sphere_two[0];
d_sample_t = cl_sphere_two[1];
d_third_t = cl_sphere_two[2];
d_detector_t = cl_sphere_two[3];
m.rw_t = cl_sphere_two[4];
m.as_t = sqr(d_sample_t / m.d_sphere_t / 2);
m.at_t = sqr(d_third_t / m.d_sphere_t / 2);
m.ad_t = sqr(d_detector_t / m.d_sphere_t / 2);
m.aw_t = 1.0 - m.as_t - m.at_t - m.ad_t;
if (cl_num_spheres == UNINITIALIZED)
m.num_spheres = 2;
}
if (cl_num_spheres != UNINITIALIZED) {
m.num_spheres = (int) cl_num_spheres;
if (m.num_spheres > 0 && m.method == UNKNOWN)
m.method = SUBSTITUTION;
}
if (cl_ru_fraction != UNINITIALIZED)
m.fraction_of_ru_in_mr = cl_ru_fraction;
if (cl_tu_fraction != UNINITIALIZED)
m.fraction_of_tu_in_mt = cl_tu_fraction;
if (cl_UR1 != UNINITIALIZED)
m.m_r = cl_UR1;
if (cl_UT1 != UNINITIALIZED)
m.m_t = cl_UT1;
if (cl_Tc != UNINITIALIZED)
m.m_u = cl_Tc;
if (cl_default_fr != UNINITIALIZED)
m.f_r = cl_default_fr;
if (cl_baffle_r != UNINITIALIZED)
m.baffle_r = cl_baffle_r;
if (cl_baffle_t != UNINITIALIZED)
m.baffle_t = cl_baffle_t;
if (cl_lambda != UNINITIALIZED)
m.lambda = cl_lambda;
Initialize_Result(m, &r, TRUE);
params = 0;
rt_total = 0;
mc_total = 0;
if (Read_Header(stdin, &m, &params) != 0)
exit(EXIT_FAILURE);
start_time = clock();
while (Read_Data_Line(stdin, &m, &r, params) == 0) {
if (cl_cos_angle != UNINITIALIZED) {
m.slab_cos_angle = cl_cos_angle;
if (cl_quadrature_points == UNINITIALIZED)
cl_quadrature_points = 12;
if (cl_quadrature_points != 12 * (cl_quadrature_points / 12)) {
fprintf(stderr, "If you use the -i option to specify an oblique incidence angle, then\n");
fprintf(stderr, "the number of quadrature points must be a multiple of 12\n");
exit(EXIT_SUCCESS);
}
}
if (cl_sample_n != UNINITIALIZED)
m.slab_index = cl_sample_n;
if (cl_slide_n != UNINITIALIZED) {
m.slab_bottom_slide_index = cl_slide_n;
m.slab_top_slide_index = cl_slide_n;
if (cl_slide_d == UNINITIALIZED && !Column_Label_Present('D')) {
if (m.slab_top_slide_thickness == 0)
m.slab_top_slide_thickness = 1.0;
if (m.slab_bottom_slide_thickness == 0)
m.slab_bottom_slide_thickness = 1.0;
}
}
if (cl_slide_OD != UNINITIALIZED) {
m.slab_bottom_slide_b = cl_slide_OD;
m.slab_top_slide_b = cl_slide_OD;
}
if (cl_sample_d != UNINITIALIZED)
m.slab_thickness = cl_sample_d;
if (cl_beam_d != UNINITIALIZED)
m.d_beam = cl_beam_d;
if (cl_slide_d != UNINITIALIZED) {
m.slab_bottom_slide_thickness = cl_slide_d;
m.slab_top_slide_thickness = cl_slide_d;
}
if (cl_slides == NO_SLIDES) {
m.slab_bottom_slide_index = 1.0;
m.slab_bottom_slide_thickness = 0.0;
m.slab_top_slide_index = 1.0;
m.slab_top_slide_thickness = 0.0;
}
if (cl_slides == ONE_SLIDE_ON_TOP || cl_slides == ONE_SLIDE_NEAR_SPHERE) {
m.slab_bottom_slide_index = 1.0;
m.slab_bottom_slide_thickness = 0.0;
}
if (cl_slides == ONE_SLIDE_ON_BOTTOM || cl_slides == ONE_SLIDE_NOT_NEAR_SPHERE) {
m.slab_top_slide_index = 1.0;
m.slab_top_slide_thickness = 0.0;
}
if (cl_slides == ONE_SLIDE_NEAR_SPHERE || cl_slides == ONE_SLIDE_NOT_NEAR_SPHERE)
m.flip_sample = 1;
else
m.flip_sample = 0;
if (cl_slides == NO_SLIDES) {
m.slab_top_slide_b = 0.0;
m.slab_bottom_slide_b = 0.0;
}
if (cl_slides == ONE_SLIDE_ON_TOP || cl_slides == ONE_SLIDE_NEAR_SPHERE)
m.slab_bottom_slide_b = 0.0;
if (cl_slides == ONE_SLIDE_ON_BOTTOM || cl_slides == ONE_SLIDE_NOT_NEAR_SPHERE)
m.slab_top_slide_b = 0.0;
if (cl_method != UNINITIALIZED)
m.method = (int) cl_method;
if (cl_rstd_r != UNINITIALIZED) {
m.rstd_r = cl_rstd_r;
m.rstd_t = cl_rstd_r;
}
if (cl_rstd_t != UNINITIALIZED) {
m.rstd_t = cl_rstd_t;
if (cl_rstd_r == UNINITIALIZED)
m.rstd_r = cl_rstd_t;
}
if (cl_rwall_r != UNINITIALIZED || cl_rwall_t != UNINITIALIZED) {
if (cl_sphere_one[0] != UNINITIALIZED || cl_sphere_two[0] != UNINITIALIZED) {
fprintf(stderr, "A wall reflectance cannot accompany a sphere description.\n");
fprintf(stderr, " -1 and -2 already carry one as their fifth value\n");
fprintf(stderr, " -w and -W are for overriding the wall in an .rxt file\n");
exit(EXIT_FAILURE);
}
}
if (cl_rwall_r != UNINITIALIZED)
m.rw_r = cl_rwall_r;
if (cl_rwall_t != UNINITIALIZED)
m.rw_t = cl_rwall_t;
if (cl_sphere_one[0] != UNINITIALIZED) {
double d_sample_r, d_third_r, d_detector_r;
m.d_sphere_r = cl_sphere_one[0];
d_sample_r = cl_sphere_one[1];
d_third_r = cl_sphere_one[2];
d_detector_r = cl_sphere_one[3];
m.rw_r = cl_sphere_one[4];
m.as_r = sqr(d_sample_r / m.d_sphere_r / 2);
m.at_r = sqr(d_third_r / m.d_sphere_r / 2);
m.ad_r = sqr(d_detector_r / m.d_sphere_r / 2);
m.aw_r = 1.0 - m.as_r - m.at_r - m.ad_r;
m.d_sphere_t = m.d_sphere_r;
m.as_t = m.as_r;
m.at_t = m.at_r;
m.ad_t = m.ad_r;
m.aw_t = m.aw_r;
m.rw_t = m.rw_r;
if (cl_num_spheres == UNINITIALIZED)
m.num_spheres = 1;
}
if (cl_sphere_two[0] != UNINITIALIZED) {
double d_sample_t, d_third_t, d_detector_t;
m.d_sphere_t = cl_sphere_two[0];
d_sample_t = cl_sphere_two[1];
d_third_t = cl_sphere_two[2];
d_detector_t = cl_sphere_two[3];
m.rw_t = cl_sphere_two[4];
m.as_t = sqr(d_sample_t / m.d_sphere_t / 2);
m.at_t = sqr(d_third_t / m.d_sphere_t / 2);
m.ad_t = sqr(d_detector_t / m.d_sphere_t / 2);
m.aw_t = 1.0 - m.as_t - m.at_t - m.ad_t;
if (cl_num_spheres == UNINITIALIZED)
m.num_spheres = 2;
}
if (cl_num_spheres != UNINITIALIZED) {
m.num_spheres = (int) cl_num_spheres;
if (m.num_spheres > 0 && m.method == UNKNOWN)
m.method = SUBSTITUTION;
}
if (cl_ru_fraction != UNINITIALIZED)
m.fraction_of_ru_in_mr = cl_ru_fraction;
if (cl_tu_fraction != UNINITIALIZED)
m.fraction_of_tu_in_mt = cl_tu_fraction;
if (cl_UR1 != UNINITIALIZED)
m.m_r = cl_UR1;
if (cl_UT1 != UNINITIALIZED)
m.m_t = cl_UT1;
if (cl_Tc != UNINITIALIZED)
m.m_u = cl_Tc;
if (cl_default_fr != UNINITIALIZED)
m.f_r = cl_default_fr;
if (cl_baffle_r != UNINITIALIZED)
m.baffle_r = cl_baffle_r;
if (cl_baffle_t != UNINITIALIZED)
m.baffle_t = cl_baffle_t;
if (cl_lambda != UNINITIALIZED)
m.lambda = cl_lambda;
{
double ur1 = 0;
double ut1 = 0;
double uru = 0;
double utu = 0;
double mu_a = 0;
double mu_sp = 0;
double LR = 0;
double LT = 0;
int skip = FALSE;
if (cl_wave_limit[0] != UNINITIALIZED) {
if (m.lambda != 0) {
if (m.lambda < cl_wave_limit[0])
skip = TRUE;
if (m.lambda > cl_wave_limit[1])
skip = TRUE;
}
}
if (Debug(DEBUG_ANY) && !skip) {
fprintf(stderr, "\n-------------------NEXT DATA POINT---------------------\n");
if (m.lambda != 0)
fprintf(stderr, "lambda=%6.1f ", m.lambda);
fprintf(stderr, "MR=%8.5f MT=%8.5f\n\n", m.m_r, m.m_t);
if (skip)
fprintf(stderr, "skipping, wavelength out of range.\n");
}
if (!skip) {
rt_total++;
Initialize_Result(m, &r, FALSE);
if (cl_quadrature_points != UNINITIALIZED)
r.method.quad_pts = cl_quadrature_points;
else
r.method.quad_pts = 8;
if (cl_default_a != UNINITIALIZED)
r.default_a = cl_default_a;
if (cl_default_mua != UNINITIALIZED) {
r.default_mua = cl_default_mua;
if (cl_sample_d != UNINITIALIZED)
r.default_ba = cl_default_mua * cl_sample_d;
else
r.default_ba = cl_default_mua * m.slab_thickness;
}
if (cl_default_b != UNINITIALIZED)
r.default_b = cl_default_b;
if (cl_default_g != UNINITIALIZED)
r.default_g = cl_default_g;
if (cl_tolerance != UNINITIALIZED) {
r.tolerance = cl_tolerance;
r.MC_tolerance = cl_tolerance;
}
if (cl_mus0 != UNINITIALIZED) {
if (m.lambda != 0) {
cl_default_mus = cl_mus0 * pow(m.lambda / cl_mus0_lambda, cl_mus0_pwr);
}
else {
fprintf(stderr, "Seems like you want to constrain scattering to a power law.\n");
fprintf(stderr, "Unfortunately, there is no wavelength so this cannot be done.\n");
}
}
if (cl_default_mus != UNINITIALIZED) {
r.default_mus = cl_default_mus;
if (cl_sample_d != UNINITIALIZED)
r.default_bs = cl_default_mus * cl_sample_d;
else
r.default_bs = cl_default_mus * m.slab_thickness;
}
if (cl_default_musp != UNINITIALIZED) {
if (cl_default_g != UNINITIALIZED)
r.default_mus = cl_default_musp / (1.0 - cl_default_g);
else
r.default_mus = cl_default_musp;
if (cl_sample_d != UNINITIALIZED)
r.default_bs = r.default_mus * cl_sample_d;
else
r.default_bs = r.default_mus * m.slab_thickness;
}
if (cl_search != UNINITIALIZED)
r.search = cl_search;
if (cl_method == COMPARISON && m.d_sphere_r != 0 && m.as_r == 0) {
fprintf(stderr, "A dual-beam measurement is specified, but no port sizes.\n");
fprintf(stderr, "You might forsake the -X option and use zero spheres (which gives\n");
fprintf(stderr, "the same result except lost light is not taken into account).\n");
fprintf(stderr, "Alternatively, bite the bullet and enter your sphere parameters,\n");
fprintf(stderr, "with the knowledge that only the beam diameter and sample port\n");
fprintf(stderr, "diameter will be used to estimate lost light from the edges.\n");
exit(EXIT_SUCCESS);
}
if (cl_method == COMPARISON && m.num_spheres == 2) {
fprintf(stderr, "A dual-beam measurement is specified, but a two sphere experiment\n");
fprintf(stderr, "is specified. Since this seems impossible, I will make it\n");
fprintf(stderr, "impossible for you unless you specify 0 or 1 sphere.\n");
exit(EXIT_SUCCESS);
}
if (cl_method == COMPARISON && m.f_r != 0) {
fprintf(stderr, "A dual-beam measurement is specified, but a fraction of light\n");
fprintf(stderr, "is specified to hit the sphere wall first. This situation\n");
fprintf(stderr, "is not supported by iad. Sorry.\n");
exit(EXIT_SUCCESS);
}
if (rt_total == 1 && cl_verbosity > 0) {
Write_Header(m, r, params, command_line);
if (MAX_MC_iterations > 0) {
if (n_photons >= 0)
fprintf(stdout, "# Photons used to estimate lost light = %ld\n", n_photons);
else
fprintf(stdout, "# Time used to estimate lost light = %ld ms\n", -n_photons);
}
else
fprintf(stdout, "# Photons used to estimate lost light = 0\n");
fprintf(stdout, "#\n");
print_results_header(stdout);
}
m.lost_r.direct = 0;
m.lost_r.diffuse = 0;
m.lost_t.diffuse = 0;
m.lost_t.direct = 0;
m.utu_lost = 0;
Inverse_RT(m, &r);
calculate_coefficients(m, r, &LR, &LT, &mu_sp, &mu_a);
if (m.num_spheres > 0 && r.found && r.error == IAD_NO_ERROR) {
if (Debug(DEBUG_LOST_LIGHT)) {
print_results_header(stderr);
print_optical_property_result(stderr, m, r, LR, LT, mu_a, mu_sp, rt_total);
}
{
int mc_failed = 0;
int mc_unreachable = 0;
int pinned = 0;
struct measure_type good_m = m;
struct invert_type good_r = r;
double mc_prev_a = r.slab.a;
double mc_prev_b = r.slab.b;
double mc_prev_g = r.slab.g;
{
double b0 = (r.slab.b < 1e8) ? r.slab.b : 1.0;
r.mc_simplex_a_step = 1e-3;
r.mc_simplex_b_step = 1e-2 * (b0 > 1.0 ? b0 : 1.0);
r.mc_simplex_g_step = 1e-3;
}
int has_prev_diff_ur1_lost = 0;
double prev_abs_diff_ur1_lost = 0.0;
while (r.MC_iterations < MAX_MC_iterations) {
long n_photons_this;
double last_mu_sp, last_mu_a;
struct lost_type current_r, current_t;
double current_utu_lost;
double diff_ur1_lost, diff_ut1_lost, diff_uru_lost, diff_utu_lost;
double diff_uru_lost_t;
double factor = 0.3;
int too_much_lost;
double tol = r.MC_tolerance;
calculate_coefficients(m, r, &LR, &LT, &mu_sp, &mu_a);
last_mu_sp = mu_sp;
last_mu_a = mu_a;
if (Debug(DEBUG_ITERATIONS) || Debug(DEBUG_A_LITTLE)) {
fprintf(stderr, "\n------------- Monte Carlo Iteration %d -----------------\n",
r.MC_iterations + 1);
}
if (n_photons < 0)
n_photons_this = n_photons;
else if (!has_prev_diff_ur1_lost || prev_abs_diff_ur1_lost > 0.01)
n_photons_this = (n_photons / 10 > 10000) ? n_photons / 10 : 10000;
else if (prev_abs_diff_ur1_lost > 0.001)
n_photons_this = n_photons;
else
n_photons_this = (n_photons * 5 < 10000000) ? n_photons * 5 : 10000000;
MC_Lost(m, r, n_photons_this, &ur1, &ut1, &uru, &utu,
&current_r, &current_t, &current_utu_lost);
diff_ur1_lost = current_r.direct - m.lost_r.direct;
diff_uru_lost = current_r.diffuse - m.lost_r.diffuse;
diff_uru_lost_t = current_t.diffuse - m.lost_t.diffuse;
diff_ut1_lost = current_t.direct - m.lost_t.direct;
diff_utu_lost = current_utu_lost - m.utu_lost;
prev_abs_diff_ur1_lost = fabs(diff_ur1_lost);
has_prev_diff_ur1_lost = 1;
if (fabs(diff_ur1_lost) > 0.001 || fabs(diff_ut1_lost) > 0.001)
too_much_lost = 1;
else
too_much_lost = 0;
m.lost_r.direct += factor * diff_ur1_lost;
m.lost_r.diffuse += factor * diff_uru_lost;
m.lost_t.diffuse += factor * diff_uru_lost_t;
m.lost_t.direct += factor * diff_ut1_lost;
m.utu_lost += factor * diff_utu_lost;
mc_total++;
r.MC_iterations++;
Inverse_RT(m, &r);
{
double new_a = r.slab.a, new_b = r.slab.b, new_g = r.slab.g;
double da = fabs(new_a - mc_prev_a);
double db = fabs(new_b - mc_prev_b);
double dg = fabs(new_g - mc_prev_g);
double b_safe = (new_b < 1e8) ? new_b : 1.0;
double a_fixed = 1e-3;
double b_fixed = 1e-2 * (b_safe > 1.0 ? b_safe : 1.0);
double g_fixed = 1e-3;
r.mc_simplex_a_step = (da < 1e-5) ? 1e-5 : (da > a_fixed ? a_fixed : da);
r.mc_simplex_b_step = (db < 1e-4) ? 1e-4 : (db > b_fixed ? b_fixed : db);
r.mc_simplex_g_step = (dg < 1e-5) ? 1e-5 : (dg > g_fixed ? g_fixed : dg);
mc_prev_a = new_a;
mc_prev_b = new_b;
mc_prev_g = new_g;
}
calculate_coefficients(m, r, &LR, &LT, &mu_sp, &mu_a);
if (r.found && r.error == IAD_NO_ERROR) {
good_m = m;
good_r = r;
}
if (0) {
fprintf(stderr, "%2d %2d %2d | %7.4f %7.4f %7.4f | %7.4f %7.4f %7.4f\n",
r.MC_iterations, too_much_lost, r.found,
m.m_r, current_r.direct, m.lost_r.direct, m.m_t, current_t.direct, m.lost_t.direct);
}
if (Debug(DEBUG_LOST_LIGHT))
print_optical_property_result(stderr, m, r, LR, LT, mu_a, mu_sp, rt_total);
else
print_dot(start_time, r.error, mc_total, FALSE, cl_verbosity);
if (mu_a <= 0.0 && too_much_lost)
pinned++;
else
pinned = 0;
if (pinned >= 2) {
{
double b_thin = 1e-6;
double b_t, mr_at, mt_at, mr_thin, mt_thin, mr_dn, mt_dn, mr_up, mt_up;
int unreachable = 0;
bright_mr_mt(m, r, b_thin, &mr_thin, &mt_thin);
if (m.m_t > mt_thin) {
{
unreachable = 1;
b_t = b_thin;
mr_at = mr_thin;
mt_at = mt_thin;
}
}
else {
b_t = solve_for_b_from_mt(m, r, m.m_t);
bright_mr_mt(m, r, b_t, &mr_at, &mt_at);
bright_mr_mt(m, r, b_t * 0.9, &mr_dn, &mt_dn);
bright_mr_mt(m, r, b_t * 1.1, &mr_up, &mt_up);
if (m.m_r > mr_at && mr_up > mr_at && mr_at > mr_dn)
unreachable = 1;
}
if (unreachable) {
r.error = IAD_UNREACHABLE_WITH_LOST_LIGHT;
if (Debug(DEBUG_LOST_LIGHT) || Debug(DEBUG_A_LITTLE) || Debug(DEBUG_ITERATIONS)) {
fprintf(stderr, "The lost light puts the measurements out of reach.\n");
fprintf(stderr, " with no absorption at b = %.4f the model gives\n",
b_t);
fprintf(stderr, " M_R %8.5f measured %8.5f short by %8.5f\n", mr_at,
m.m_r, m.m_r - mr_at);
fprintf(stderr, " M_T %8.5f measured %8.5f\n", mt_at, m.m_t);
}
}
}
if (r.error == IAD_UNREACHABLE_WITH_LOST_LIGHT) {
mc_unreachable = 1;
if (Debug(DEBUG_ITERATIONS))
fprintf(stderr,
"absorption is spent and the data is out of reach — stopping\n");
break;
}
}
if (r.found) {
if (fabs(last_mu_a - mu_a) > tol) {
if (Debug(DEBUG_ITERATIONS))
fprintf(stderr, "Repeat MC because mua is still changing\n");
continue;
}
if (fabs(last_mu_sp - mu_sp) > tol) {
if (Debug(DEBUG_ITERATIONS))
fprintf(stderr, "Repeat MC because musp is still changing\n");
continue;
}
if (too_much_lost) {
if (Debug(DEBUG_ITERATIONS))
fprintf(stderr, "Repeat MC because mua and musp are still changing\n");
continue;
}
if (Debug(DEBUG_ITERATIONS))
fprintf(stderr, "found!\n");
break;
}
else {
if (Debug(DEBUG_ITERATIONS))
fprintf(stderr, "AD did not converge — stopping MC loop\n");
mc_failed = 1;
break;
}
}
if (mc_unreachable) {
r.found = 0;
r.error = IAD_UNREACHABLE_WITH_LOST_LIGHT;
{
int why = r.error;
struct lost_type lost_r = m.lost_r;
struct lost_type lost_t = m.lost_t;
double utu_lost = m.utu_lost;
m = good_m;
m.lost_r = lost_r;
m.lost_t = lost_t;
m.utu_lost = utu_lost;
r = good_r;
r.error = why;
r.found = 0;
}
}
else if (mc_failed) {
r.found = 0;
{
double b_thin = 1e-6;
double b_t, mr_at, mt_at, mr_thin, mt_thin, mr_dn, mt_dn, mr_up, mt_up;
int unreachable = 0;
bright_mr_mt(m, r, b_thin, &mr_thin, &mt_thin);
if (m.m_t > mt_thin) {
{
unreachable = 1;
b_t = b_thin;
mr_at = mr_thin;
mt_at = mt_thin;
}
}
else {
b_t = solve_for_b_from_mt(m, r, m.m_t);
bright_mr_mt(m, r, b_t, &mr_at, &mt_at);
bright_mr_mt(m, r, b_t * 0.9, &mr_dn, &mt_dn);
bright_mr_mt(m, r, b_t * 1.1, &mr_up, &mt_up);
if (m.m_r > mr_at && mr_up > mr_at && mr_at > mr_dn)
unreachable = 1;
}
if (unreachable) {
r.error = IAD_UNREACHABLE_WITH_LOST_LIGHT;
if (Debug(DEBUG_LOST_LIGHT) || Debug(DEBUG_A_LITTLE) || Debug(DEBUG_ITERATIONS)) {
fprintf(stderr, "The lost light puts the measurements out of reach.\n");
fprintf(stderr, " with no absorption at b = %.4f the model gives\n", b_t);
fprintf(stderr, " M_R %8.5f measured %8.5f short by %8.5f\n",
mr_at, m.m_r, m.m_r - mr_at);
fprintf(stderr, " M_T %8.5f measured %8.5f\n", mt_at, m.m_t);
}
}
}
if (r.error == IAD_NO_ERROR)
r.error = IAD_MC_DID_NOT_CONVERGE;
{
int why = r.error;
struct lost_type lost_r = m.lost_r;
struct lost_type lost_t = m.lost_t;
double utu_lost = m.utu_lost;
m = good_m;
m.lost_r = lost_r;
m.lost_t = lost_t;
m.utu_lost = utu_lost;
r = good_r;
r.error = why;
r.found = 0;
}
}
}
}
if (!r.found && r.error == IAD_NO_ERROR)
r.error = IAD_SEARCH_STALLED;
calculate_coefficients(m, r, &LR, &LT, &mu_sp, &mu_a);
print_optical_property_result(stdout, m, r, LR, LT, mu_a, mu_sp, rt_total);
if (r.error != IAD_NO_ERROR) {
any_error = 1;
last_error = r.error;
}
if (Debug(DEBUG_ANY))
print_long_error(r.error);
else
print_dot(start_time, r.error, mc_total, TRUE, cl_verbosity);
}
}
}
file_index++;
if (file_index < file_count)
goto read_next_file;
if (cl_grid_calc != UNINITIALIZED) {
double aa[] = { 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.85,
0.9, 0.93, 0.95, 0.97, 0.98, 0.99, 0.995, 1.0
};
double bb[] = { 0, 0.05, 0.1, 0.2, 0.3, 0.5, 0.7, 1.0, 1.5, 2.0,
3.0, 5.0, 7.0, 10.0, 20.0, 30.0, 50.0, 100.0
};
int na = (int) (sizeof(aa) / sizeof(aa[0]));
int nb = (int) (sizeof(bb) / sizeof(bb[0]));
double grid_step = 1e-3;
int use_mc = (MAX_MC_iterations != 0 && m.num_spheres != 0);
double g;
int i, j, k;
int count = 0;
FILE *grid;
grid = fopen(g_grid_name, "w");
if (grid == NULL) {
fprintf(stderr, "Could not open grid file '%s' for output\n", g_out_name);
exit(EXIT_FAILURE);
}
m.lost_r.direct = 0;
m.lost_r.diffuse = 0;
m.lost_t.diffuse = 0;
m.lost_t.direct = 0;
m.utu_lost = 0;
if (r.default_g != UNINITIALIZED) {
g = r.default_g;
}
else if (r.found) {
g = r.slab.g;
}
else {
g = 0;
}
fprintf(grid, "# %s (g=%6.4f)\n", command_line, g);
fprintf(grid, "# a' b' g M_R M_T rho12\n");
fprintf(stderr, "\ndoing grid calculation\n");
for (i = 0; i < na; i++) {
for (j = 0; j < nb; j++) {
double ap[3], bp[3], m_r[3], m_t[3];
ap[0] = aa[i];
bp[0] = bb[j];
ap[1] = (aa[i] + grid_step <= 1) ? aa[i] + grid_step : aa[i] - grid_step;
bp[1] = bb[j];
ap[2] = aa[i];
bp[2] = bb[j] * (1 + grid_step);
if (use_mc) {
double ur1, ut1, uru, utu;
r.slab.a = ap[0] / (1 - g + ap[0] * g);
r.slab.b = bp[0] / (1 - r.slab.a * g);
r.slab.g = g;
r.a = r.slab.a;
r.b = r.slab.b;
r.g = g;
MC_Lost(m, r, 100000, &ur1, &ut1, &uru, &utu, &m.lost_r, &m.lost_t, &m.utu_lost);
}
for (k = 0; k < 3; k++) {
r.slab.a = ap[k] / (1 - g + ap[k] * g);
r.slab.b = bp[k] / (1 - r.slab.a * g);
r.slab.g = g;
r.a = r.slab.a;
r.b = r.slab.b;
r.g = g;
Calculate_MR_MT(m, r, use_mc ? MC_USE_EXISTING : MC_NONE, TRUE, &m_r[k], &m_t[k]);
}
fprintf(grid, "%10.5f, %10.5f, %10.5f, %10.5f, %10.5f, ", ap[0], bp[0], g, m_r[0], m_t[0]);
{
double dMR_da, dMT_da, dMR_db, dMT_db, u_r, u_t, v_r, v_t, uu, vv;
u_r = u_t = v_r = v_t = 0;
uu = 0;
vv = 0;
if (bp[0] > 0) {
dMR_da = (m_r[1] - m_r[0]) / (ap[1] - ap[0]);
dMT_da = (m_t[1] - m_t[0]) / (ap[1] - ap[0]);
dMR_db = (m_r[2] - m_r[0]) / (bp[2] - bp[0]);
dMT_db = (m_t[2] - m_t[0]) / (bp[2] - bp[0]);
u_r = dMR_db - ap[0] / bp[0] * dMR_da;
u_t = dMT_db - ap[0] / bp[0] * dMT_da;
v_r = dMR_db + (1 - ap[0]) / bp[0] * dMR_da;
v_t = dMT_db + (1 - ap[0]) / bp[0] * dMT_da;
uu = u_r * u_r + u_t * u_t;
vv = v_r * v_r + v_t * v_t;
}
if (uu > 1e-20 && vv > 1e-20)
fprintf(grid, "%10.5f\n", -(u_r * v_r + u_t * v_t) / sqrt(uu * vv));
else
fprintf(grid, "%10s\n", "nan");
}
count++;
fprintf(stderr, "*");
if (count % 10 == 0)
fprintf(stderr, " ");
if (count % 50 == 0)
fprintf(stderr, "\n");
}
}
fclose(grid);
fprintf(stderr, "\n");
}
if (cl_verbosity > 0)
fprintf(stderr, "\n\n");
if (any_error && cl_verbosity > 1)
print_error_legend();
exit(EXIT_SUCCESS);
}

iad_main.c at commit 8a1da20, under MIT · at the source

Overview

Authors: Silvére Ségaud¹, Charlie Budd¹, Matthew Elliot^1,2, Graeme J Stasiuk³, Jonathan Shapey^1,2, Yijing Xie¹, Tom Vercauteren¹

ORCID iDs: Matthew Elliot, Graeme J Stasiuk, Yijing Xie, Tom Vercauteren

Institutions: King's College London (United Kingdom); King's College Hospital NHS Foundation Trust (United Kingdom)

Journal: Physics in medicine and biology, volume 71, issue 5, article 055006

Dates: received 17 October 2025; accepted 13 February 2026; published online 6 March 2026; in print 14 March 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1088/1361-6560/ae45e6 · PMID 41687253 · PMCID PMC12965099 · OpenAlex W4416056020

Open access: hybrid, a free copy (OpenAlex)

Preprint: arxiv:2510.18387

Status: code verified

Categories: human (organism), other condition (population), clinical / translational (subfield)

Methods: Connectivity, Machine learning

Keywords: glioma, fluorescence guided surgery, quantitative fluorescence, tissue-mimicking phantom, 5-ALA-PpIX, PpIX620, PpIX635

MeSH: Aminolevulinic Acid*, Optical Imaging*, Phantoms, Imaging*, Protoporphyrins*, Brain Neoplasms, Glioma, Humans, Spectrometry, Fluorescence (* major topic)

Topic: Photodynamic Therapy Research Studies (Pulmonary and Respiratory Medicine, Medicine), according to OpenAlex

Funding: Wellcome Trust (WT223880/Z/21/Z); National Institute for Health Research (NIHR) (NIHR202114); National Institute for Health and Care Research (NIHR202114); Wellcome/EPSRC (NS/A000049/1)

Citations: not cited yet (Europe PMC); 55 references in the paper

Abstract

Quantification of protoporphyrin IX (PpIX) fluorescence in human brain tumours has the potential to significantly improve patient outcomes in neuro-oncology, but represents a formidable imaging challenge. Protoporphyrin is a biological molecule which interacts with the tissue micro-environment to form two photochemical states in glioma. Each exhibits markedly different quantum efficiencies, with distinct but overlapping emission spectra that also overlap with tissue autofluorescence. Fluorescence emission is known to be distorted by the intrinsic optical properties of tissue, coupled with marked intra-tumoural heterogeneity as a hallmark of glioma tumours. Existing quantitative fluorescence systems are developed and validated using simplified phantoms that do not simultaneously mimic the complex interactions between fluorophores and tissue optical properties or micro-environment. Consequently, existing systems risk introducing systematic errors into PpIX quantification when used in tissue. In this work, we introduce a novel pipeline for quantification of PpIX in glioma, which robustly differentiates both emission states from background autofluorescence without reliance on a priori spectral information, and accounts for variations in their quantum efficiency. Unmixed PpIX emission forms are then corrected for wavelength-dependent optical distortions and weighted for accurate quantification. Significantly, this pipeline is developed and validated using novel tissue-mimicking phantoms replicating the optical properties of glioma tissues and photochemical variability of PpIX fluorescence in glioma. Our workflow achieves strong correlation with ground-truth PpIX concentrations (R2=0.918±0.002), demonstrating its potential for robust, quantitative PpIX fluorescence imaging in clinical settings.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

scottprahl/iad

License: MIT

State: the link answers, verified on 30 September 2026

Evidence: files inventoried

Commit: 8a1da205b9e198011fa6903bd4d896affce228dd, 26 September 2026

Languages: C (37), C/C++ (31), Shell (16), Python (3), Perl (2)

Size: 230 files, 89 scripts

Software Heritage: archived

Found in: the text, “Footnotes”

Holds: README, license file, CITATION.cff, tests, continuous integration, documentation

Not found: environment file

Availability: 1 check, the latest on 30 September 2026: the link answers

30 September 2026: the link answers

91 files

rxtformat.py, Python, 731 lines
src/ad_bound.c, C, 152 lines
src/ad_bound.h, C/C++, 21 lines
src/ad_chapter.c, C, 181 lines
src/ad_cone.c, C, 223 lines
src/ad_cone.h, C/C++, 15 lines
src/ad_cone_ez.h, C/C++, 13 lines
src/ad_cone_test.c, C, 429 lines
src/ad_doubl.c, C, 199 lines
src/ad_doubl.h, C/C++, 22 lines
src/ad_frsnl.c, C, 205 lines
src/ad_frsnl.h, C/C++, 17 lines
src/ad_globl.c, C, 217 lines
src/ad_globl.h, C/C++, 57 lines
src/ad_layers.c, C, 166 lines
src/ad_layers.h, C/C++, 17 lines
src/ad_layers_test.c, C, 295 lines
src/ad_main.c, C, 440 lines
src/ad_matrx.c, C, 314 lines
src/ad_matrx.h, C/C++, 27 lines
src/ad_oblique_test.c, C, 117 lines
src/ad_phase.c, C, 85 lines
src/ad_phase.h, C/C++, 3 lines
src/ad_prime.c, C, 483 lines
src/ad_prime.h, C/C++, 22 lines
src/ad_radau.c, C, 250 lines
src/ad_radau.h, C/C++, 3 lines
src/ad_start.c, C, 372 lines
src/ad_start.h, C/C++, 11 lines
src/cobweb.pl, Perl, 8 lines
src/iad_agrid.c, C, 780 lines
src/iad_agrid.h, C/C++, 7 lines
src/iad_calc.c, C, 664 lines
src/iad_calc.h, C/C++, 55 lines
src/iad_find.c, C, 865 lines
src/iad_find.h, C/C++, 21 lines
src/iad_io.c, C, 724 lines
src/iad_io.h, C/C++, 9 lines
src/iad_main.c, C, 3,065 lines, 2 matches
src/iad_pub.c, C, 829 lines
src/iad_pub.h, C/C++, 21 lines
src/iad_type.h, C/C++, 191 lines
src/iad_util.c, C, 237 lines
src/iad_util.h, C/C++, 30 lines
src/lib_test.c, C, 401 lines
src/mc_lost.c, C, 670 lines
src/mc_lost.h, C/C++, 21 lines
src/mc_lost_main.c, C, 226 lines
src/mc_test.c, C, 393 lines
src/nr_amoeb.c, C, 220 lines
src/nr_amoeb.h, C/C++, 2 lines
src/nr_amotr.c, C, 27 lines
src/nr_amotr.h, C/C++, 1 line
src/nr_brent.c, C, 87 lines
src/nr_brent.h, C/C++, 1 line
src/nr_gaulg.c, C, 33 lines
src/nr_gaulg.h, C/C++, 1 line
src/nr_hj.c, C, 182 lines
src/nr_hj.h, C/C++, 1 line
src/nr_mnbrk.c, C, 74 lines
src/nr_mnbrk.h, C/C++, 2 lines
src/nr_rtsaf.c, C, 67 lines
src/nr_rtsaf.h, C/C++, 1 line
src/nr_util.c, C, 290 lines
src/nr_util.h, C/C++, 57 lines
src/nr_zbrak.c, C, 22 lines
src/nr_zbrak.h, C/C++, 1 line
src/nr_zbrent.c, C, 81 lines
src/nr_zbrent.h, C/C++, 1 line
src/toDOS.pl, Perl, 3 lines
src/version.h, C/C++, 2 lines
tests/cli/compare_numeric.py, Python, 94 lines
tests/cli/lib.sh, Shell, 111 lines
tests/cli/run_cli_tests.sh, Shell, 39 lines
tests/cli/test_ad_basic.sh, Shell, 32 lines
tests/cli/test_ez_rt_callers.sh, Shell, 65 lines
tests/cli/test_iad_basic.sh, Shell, 95 lines
tests/cli/test_iad_files.sh, Shell, 104 lines
tests/cli/test_iad_forward.sh, Shell, 23 lines
tests/cli/test_iad_no_sphere.sh, Shell, 19 lines
tests/cli/test_iad_one_sphere.sh, Shell, 10 lines
tests/cli/test_iad_one_sphere_no_mc.sh, Shell, 14 lines
tests/cli/test_iad_options.sh, Shell, 257 lines
tests/cli/test_iad_order.sh, Shell, 65 lines
tests/cli/test_iad_slides.sh, Shell, 99 lines
tests/cli/test_iad_two_sphere.sh, Shell, 10 lines
tests/cli/test_iad_two_sphere_no_mc.sh, Shell, 14 lines
tests/cli/test_mc_slides.sh, Shell, 120 lines
tools/clean_generated.py, Python, 213 lines
License, License, 7 lines
README.rst, Text, 203 lines

At the source: github.com/scottprahl/iad

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
89 scripts, each with its path and the digest of its content;
2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability statement

All data that support the findings of this study are included within the article (and any supplementary information files). Data will be available from 20 March 2026.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 7 keywords, 8 MeSH terms, 4 funders, 54 references.

Cite

This paper

Ségaud, S., Budd, C., Elliot, M., Stasiuk, G. J., Shapey, J., Xie, Y., & Vercauteren, T. (2026). Quantification of dual-state 5-ALA-induced PpIX fluorescence: methodology and validation in tissue-mimicking phantoms. Physics in medicine and biology, 71(5), 055006. https://doi.org/10.1088/1361-6560/ae45e6

BibTeX

@article{segaud2026quantification, author = {Ségaud, Silvére and Budd, Charlie and Elliot, Matthew and Stasiuk, Graeme J and Shapey, Jonathan and Xie, Yijing and Vercauteren, Tom}, title = {{Quantification of dual-state 5-ALA-induced PpIX fluorescence: methodology and validation in tissue-mimicking phantoms}}, journal = {Physics in medicine and biology}, year = {2026}, month = mar, volume = {71}, number = {5}, pages = {055006}, publisher = {IOP Publishing}, issn = {0031-9155}, doi = {10.1088/1361-6560/ae45e6}, url = {https://doi.org/10.1088/1361-6560/ae45e6}, pmid = {41687253}, pmcid = {PMC12965099} }

RIS

TY - JOUR AU - Ségaud, Silvére AU - Budd, Charlie AU - Elliot, Matthew AU - Stasiuk, Graeme J AU - Shapey, Jonathan AU - Xie, Yijing AU - Vercauteren, Tom TI - Quantification of dual-state 5-ALA-induced PpIX fluorescence: methodology and validation in tissue-mimicking phantoms T2 - Physics in medicine and biology J2 - Phys Med Biol PY - 2026 DA - 2026/03/06 VL - 71 IS - 5 SP - 055006 SN - 0031-9155 PB - IOP Publishing DO - 10.1088/1361-6560/ae45e6 UR - https://doi.org/10.1088/1361-6560/ae45e6 LA - en ER -

CSL-JSON

{ "id": "10.1088/1361-6560/ae45e6", "type": "article-journal", "title": "Quantification of dual-state 5-ALA-induced PpIX fluorescence: methodology and validation in tissue-mimicking phantoms", "container-title": "Physics in medicine and biology", "author": [ { "family": "Ségaud", "given": "Silvére" }, { "family": "Budd", "given": "Charlie" }, { "family": "Elliot", "given": "Matthew" }, { "family": "Stasiuk", "given": "Graeme J" }, { "family": "Shapey", "given": "Jonathan" }, { "family": "Xie", "given": "Yijing" }, { "family": "Vercauteren", "given": "Tom" } ], "container-title-short": "Phys Med Biol", "volume": "71", "issue": "5", "page": "055006", "DOI": "10.1088/1361-6560/ae45e6", "PMID": "41687253", "PMCID": "PMC12965099", "ISSN": "0031-9155", "publisher": "IOP Publishing", "URL": "https://doi.org/10.1088/1361-6560/ae45e6", "language": "en", "issued": { "date-parts": [ [ 2026, 3, 6 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.