Biomedical Image Signal Processing for Reflection-Based Imaging

1. Introduction

Optical coherence tomography (OCT), ultrasound and other reflection based biomedical imaging technologies involve image signal processing that is primarily a filtering, digitizing and summing process so that the tissue cross-section can be visualized. In particular, in OCT, a series of adjacent one dimensional in-vivo axial interferograms (A-scan) are summed to form a two dimensional (B-scan) reflection map or reflectogram. Further graphical combinations can add adjacent B-scan together to form three dimensional C-scans. A physician can make a subjective interpretation and evaluation from the B and/or C-scans that may lead to actions impacting on the patient’s prognosis. More objective information can be obtained by using backwards fitting models (BFM) that fit tissue characteristics, including layer depth and reflectivity, to imaged tissue A-scans, returning values that are not significantly different to the actual values.

BFMs can be either primarily deterministic or stochastic. One example of the latter is a genetic algorithm model (GAM). In this chapter, a GAM is characterised for its degree of precision and accuracy to retrieve the depth and reflectivity profile of a simulated A-scan of a virtual tissue model with defined tissue layer depths and layer refractive indicies. This stochastic model intrinsically evolves successively more precise and accurate generations of solutions, in accordance with certain software defined specific and random selection control parameters. This chapter presents a characterisation of the dependence of the accuracy and precision of the estimated layer depths and reflectivities, returned by this version of the GAM, on the number of generations, and the values of the GAM’s selection-control parameters.

2. Theory

2.1. Optical coherence tomography

Optical coherence tomography (OCT) is a medical imaging technique that is fundamentally an application of low coherence interferometry (Drexler et al 2008), (Fercher 1996), (Izatt et al 1997). In OCT, a low coherent near infrared source is used to generate a reflection intensity map of the tissues cross section. An OCT scan can be 1 to 3 mm deep (Friebel et al 2009), depending on the tissue type and the optical properties of the light source, specifically the wavelength and intensity. Using a conventional solid state broadband light source, the axial resolution is typically 1 to 10 micrometers, depending on the bandwidth, central wavelength (Drexler et al 2008) and the sources Fourier envelope shape (Adie 2007) (Rossetti et al 2005) (Shidlovski 2008). Four key elements impact on the OCT detector signal: the light source, the optical delay line (ODL), type of interferometer used and the sample characteristics.

Figure 1.

Operating principle of a Michelson interferometer type OCT system.

The light source is low-coherent in that it has a broad frequency bandwidth. All the light frequencies interfere with each other resulting in a self modulated light source where the width of the individual peak in the time domain is proportional to the image axial resolution, measured by the coherence length (L_c) of the source (Eq. (1)) divided by the average tissue refractive index. The L_c can be determined from the spectral characteristics of the source using,

$$L_C=\frac{\ln 4.{\lambda }_0^2}{\pi \Delta \lambda },$$E1

where λ₀ is the source central wavelength, and Δλ is the spectral Full Width at Half Maximum (FWHM) of the power spectrum, assuming it has a Gaussian spectral profile. Typical OCT sources have λ₀ at 840 nm and 1320 nm, with Δλ equal to 60 nm. From Equation (1), these values give coherence lengths of 5.2 μm and 12.8 μm, respectively.

Because the inverse Fourier Transform of a perfectly Gaussian spectrum in the frequency domain is itself a Gaussian in the time domain, a perfectly Gaussian spectrum will have one peak and no repeated peaks; ideal for stratified samples. The less Gaussian the source spectrum in the frequency domain, the more frequently that peaks appear in the time domain. A Michelson or Mach Zehnder interferometer can then be used to manipulate the introduced light to acquire an interference pattern of the multi-layered sample’s reflected beam and the reflected reference ODL beam at the detector (Fig. 1).

2.2. The forward model

The use of a Matlab OCT circuit simulation model to produce interferogram A-scans from a defined Gaussian spectral light source for a user defined stratified sample has been previously used to characterise typical OCT light sources (Jansz et al 2012) and optical delay lines (Jansz et al 2011). With this forwards model, backward fitting models (BFM) can now be tested to retro-fit certain sample parameters. These parameters for the present BFM version are sample layer thicknesses and reflectivities.

2.2.1. The source

Suppose that a light source emits a continuous distribution of wavelengths whose amplitude is a function of the wavelength. That is, A = A(λ). For example, if the light source is modelled as a spectrum of continuous wavelengths, with a Gaussian spectral shape defined by the peak amplitude (A₀), the peak wavelength (λ₀), and the spectral bandwidth, FWHM (Δλ), we have,

$$A(\lambda )=A_{0}Exp\left(\frac{-\ln (16){\left(\lambda -{\lambda }_0\right)}^2}{{\left(\Delta \lambda \right)}^2}\right)$$E2

,

with the laser power given by,

$$P=A_0\Delta \lambda \sqrt{\frac{\pi }{\ln (16)}}$$E3

2.2.2. Light amplitude in the interferometer

The distribution of wavelengths is first passed through a 50% mirror, with half of the light reflected (R), and half transmitted (T). The result is a reflected and transmitted wave with amplitude distributions,

$$A_R\left(\lambda \right)=A_T\left(\lambda \right)=\frac{A\left(\lambda \right)}{\sqrt{2}}$$E4

.

Once the original light source is split, it becomes important to keep track of the distance travelled by each wave, as any difference will be associated with a phase shift, and therefore a change in the interference pattern. The transmitted wave, considered to be the sample arm of the interferometer, travels to a multi-layered surface consisting of n partially reflecting interfaces, the phantom structure to be imaged. The transmitted wave goes through transmission and reflection at each interface. Upon returning from the multi-layered interface, the transmitted wave is reflected by the 50% mirror to the detector. The total distance travelled by the transmitted wave which reflects off the interface of the multi-layered structure is denoted d_i for i = 1,2,…,n. The thickness Δd_i of the i^th layer in the sample (i.e. the distance between the i^th and (i+1)^th interface) is then given by

$$\Delta d_i=\frac{d_{i+1}-d_i}{2},i=1,2,\dots ,n-1$$E5

The reflectivity of the interface is denoted by r_i for i = 1,2,…,n. We have assumed that the contribution of waves which are reflected off multiple interfaces within the multi-layered structure is negligible. Therefore it makes sense to decompose the detected component of the transmitted wave into n parts corresponding to each of the reflecting surfaces. Hence, we have,

$$A_T{\left(\lambda \right)}_i=\{\begin{array}{cc}\frac{A\left(\lambda \right)}{2}\sqrt{r_1},& i=1\\ \frac{A\left(\lambda \right)}{2}\sqrt{r_i}{\displaystyle \prod _{j=1}^{i-1}\left(1-r_j\right)},& i=2,3,\dots ,n\end{array}$$E6

The reflected wave, considered to be the reference arm of the interferometer, travels from the 50% mirror to another reflector which is moved incrementally (i.e. no Doppler effect) before being reflected back through the 50% mirror (transmitted) to a detector. The total distance travelled by the reflected wave is denoted d_n+1. Note that a shift of the moving reflector by an amount δ will result in an increase in the total distance travelled by the reflected wave of 2δ. We will denote the reflectivity of the reference arm mirror by r_n+1 (although typically r_n+1 = 1). The final amplitude distribution of the reflected wave is,

$$A_R\left(\lambda \right)=\frac{A\left(\lambda \right)}{2}\sqrt{r_{n+1}}.$$E7

In the interest of notational convenience, let us define a reflectivity factor,

$$R{F}_i=\{\begin{array}{ll}\sqrt{r_1},\hfill & i=1\hfill \\ \sqrt{r_i}{\displaystyle \prod _{j=1}^{i-1}\left(1-r_j\right)},\hfill & i=2,3,\dots ,n\hfill \\ \sqrt{r_{n+1}},\hfill & i=n+1\hfill \end{array}.$$E8

Then the amplitudes of the transmitted and reflected waves can be expressed as,

$$\begin{array}{l}{A}_T{\left(\lambda \right)}_i=\frac{A\left(\lambda \right)}{2}R{F}_i,i=1,2,\dots ,n\\ A_R\left(\lambda \right)=\frac{A\left(\lambda \right)}{2}R{F}_{n+1}.\end{array}$$E9

2.2.3. Light wave interference in the interferometer

Having defined expressions for the amplitudes of the interfering waves, we now consider the interference of these waves. Without loss of generality, we can represent all waves using a sine function with origin at the detector. The expression for the reflected wave is,

$$y_R\left(\lambda ,x,t\right)=A_R\left(\lambda \right)\sin \left(\frac{2\pi }{\lambda }\left(x+d_{n+1}-ct\right)\right),$$E10

while for the (originally) transmitted components we have,

$$y_T{\left(\lambda ,x,t\right)}_i=A_T\left(\lambda \right)\sin \left(\frac{2\pi }{\lambda }\left(x+d_i-ct\right)+{\phi }_i\pi \right),i=1,2,\dots ,n,$$E11

where ϕ_i is a phase shift indicator function given by,

$${\phi }_i=\{\begin{array}{cc}1,& n_{i-1}<n_i\\ 0,& n_{i-1}>n_i\end{array},i=1,2,\dots ,n,$$E12

where n_i denotes the refractive index of the i^th layer of the sample (not to be confused with the number of layers in the sample n). Note that n₀ denotes the refractive index of the medium in front of the sample. Note also that ϕ_n+1 = 1.

For the sake of simplifying notation, let us denote,

$$\omega =\frac{2\pi }{\lambda }\left(x-ct\right),$$E13

$${\theta }_i=\frac{2\pi d_i}{\lambda }+{\phi }_i\pi ,i=1,2,\dots ,n.$$E14

Then we can write,

$$y_R\left(\lambda ,\omega \right)=A_R\left(\lambda \right)\sin \left(\omega +{\theta }_{n+1}\right)$$E15

$$y_T{\left(\lambda ,\omega \right)}_i=A_T\left(\lambda \right)\sin \left(\omega +{\theta }_i\right),i=1,2,\dots ,n$$E16

The resultant wave arriving at the detector is therefore given by,

$$\begin{array}{c}Y\left(\lambda ,\omega \right)=y_R\left(\lambda ,\omega \right)+{\displaystyle \sum _{i=1}^n{y}_T{\left(\lambda ,\omega \right)}_i}\\ ={\displaystyle \sum _{i=1}^{n+1}\frac{A\left(\lambda \right)}{2}}R{F}_i\sin \left(\omega +{\theta }_i\right)\\ ={\displaystyle \sum _{i=1}^{n+1}\frac{A\left(\lambda \right)}{2}}R{F}_i\left[\sin \left(\omega \right)\cos \left({\theta }_i\right)+\cos \left(\omega \right)\sin \left({\theta }_i\right)\right]\\ =C_1\sin \left(\omega \right)+C_2\cos \left(\omega \right)\\ =\sqrt{C_1{}^2+C_2{}^2}\sin \left(\omega +a\tan \left(\frac{C_2}{C_1}\right)\right)\end{array}$$E17

where,

$$\begin{array}{c}{C}_1={\displaystyle \sum _{i=1}^{n+1}\frac{A\left(\lambda \right)}{2}}R{F}_i\cos \left({\theta }_i\right)\\ C_2={\displaystyle \sum _{i=1}^{n+1}\frac{A\left(\lambda \right)}{2}}R{F}_i\sin \left({\theta }_i\right)\end{array}$$E18

It follows that the amplitude of the resultant wave is given by,

$$Amplitude\left(\lambda \right)=\sqrt{C_1{}^2+C_2{}^2},$$E19

and hence that the total intensity of detected light over all wavelengths is given by,

$$I={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\left[C_1{}^2+C_2{}^2\right]}d\lambda .$$E20

Equation (20) can be expressed in a more convenient form by simplifying the integrand as follows,

$$\begin{array}{c}{C}_1{}^2+C_2{}^2={\left({\displaystyle \sum _{i=1}^{n+1}\frac{A\left(\lambda \right)}{2}}R{F}_i\cos \left({\theta }_i\right)\right)}^2+{\left({\displaystyle \sum _{i=1}^{n+1}\frac{A\left(\lambda \right)}{2}}R{F}_i\sin \left({\theta }_i\right)\right)}^2\\ =\frac{A{\left(\lambda \right)}^2}{4}\left[\begin{array}{l}{\displaystyle \sum _{i=1}^{n+1}R{F}_i{}^2{\cos }^2\left({\theta }_i\right)}+2{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=i+1}^{n+1}R{F}_{i}R{F}_j\cos \left({\theta }_i\right)\cos \left({\theta }_j\right)}}\\ +{\displaystyle \sum _{i=1}^{n+1}R{F}_i{}^2{\sin }^2\left({\theta }_i\right)}+2{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=i+1}^{n+1}R{F}_{i}R{F}_j\sin \left({\theta }_i\right)\sin \left({\theta }_j\right)}}\end{array}\right]\\ =\frac{A{\left(\lambda \right)}^2}{4}\left[{\displaystyle \sum _{i=1}^{n+1}R{F}_i{}^2+2{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=i+1}^{n+1}R{F}_{i}R{F}_j\cos \left({\theta }_i-{\theta }_j\right)}}}\right]\end{array}$$E21

$$\begin{array}{c}=\frac{A{\left(\lambda \right)}^2}{4}\left[{\displaystyle \sum _{i=1}^{n+1}R{F}_i{}^2+2{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=i+1}^{n}R{F}_{i}R{F}_j\cos \left({\theta }_i-{\theta }_j\right)+2{\displaystyle \sum _{i=1}^{n}R{F}_{i}R{F}_{n+1}\cos \left({\theta }_i-{\theta }_{n+1}\right)}}}}\right]\\ =\frac{A{\left(\lambda \right)}^2}{4}\left[\begin{array}{l}{\displaystyle \sum _{i=1}^{n+1}R{F}_i{}^2+2{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=i+1}^{n}R{F}_{i}R{F}_j\cos \left(\frac{2\pi }{\lambda }\left(d_i-d_j\right)+\pi \left({\phi }_i-{\phi }_j\right)\right)}}}\\ +2{\displaystyle \sum _{i=1}^{n}R{F}_{i}R{F}_{n+1}\cos \left(\frac{2\pi }{\lambda }\left(d_i-d_{n+1}\right)+\pi \left({\phi }_i-{\phi }_{n+1}\right)\right)}\end{array}\right]\\ =\frac{A{\left(\lambda \right)}^2}{4}\left[\begin{array}{l}{\displaystyle \sum _{i=1}^{n+1}R{F}_i{}^2+2{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=i+1}^{n}R{F}_{i}R{F}_j{\left(-1\right)}^{{\phi }_i+{\phi }_j}\cos \left(\frac{2\pi }{\lambda }\left(d_i-d_j\right)\right)}}}\\ +2{\displaystyle \sum _{i=1}^{n}R{F}_{i}R{F}_{n+1}{\left(-1\right)}^{{\phi }_i+1}\cos \left(\frac{2\pi }{\lambda }\left(d_i-d_{n+1}\right)\right)}\end{array}\right]\end{array}$$E22

Hence, the intensity I at the detector can be expressed as,

$$I=B_0+{\displaystyle \sum _{i=1}^n{B}_{i}F\left(d_i-d_{n+1}\right)},$$E23

where,

$$B_0=\frac{1}{4}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}A{\left(\lambda \right)}^{2}d\lambda }{\displaystyle \sum _{i=1}^{n+1}R{F}_i{}^2+\frac{1}{2}{\displaystyle \sum _{i=1}^{n-1}{\displaystyle \sum _{j=i+1}^{n}R{F}_{i}R{F}_j}}}{\left(-1\right)}^{{\phi }_i+{\phi }_j}F\left(d_i-d_j\right),$$E24

$$B_i=\frac{1}{2}R{F}_{i}R{F}_{n+1}{\left(-1\right)}^{{\phi }_i+1}$$E25

$$F(x)={\displaystyle {\int }_{-\infty }^{\infty }A{\left(\lambda \right)}^2}\cos \left(\frac{2\pi x}{\lambda }\right)d\lambda$$E26

Note that if A(λ) is given by equation (2) then,

$$B_0=\frac{A_0{}^2\Delta \lambda }{4}\sqrt{\frac{\pi }{2\ln (16)}}{\displaystyle \sum _{i=1}^{n+1}R{F}_i{}^2+\frac{1}{2}{\displaystyle \sum _{i=1}^{n-1}{\displaystyle \sum _{j=i+1}^{n}R{F}_{i}R{F}_j}}}{\left(-1\right)}^{{\phi }_i+{\phi }_j}F\left(d_i-d_j\right)$$E27

The expression given in equation (22) separates the intensity into a constant offset component B₀ (i.e. constant for a specific sample structure with fixed distances d₁, d₂,…,d_n and reflectivity’s r₁, r₂,…,r_n), and an interference component for each of the layers given by B_iF(d_i – d_n+1). The coefficient B_i contains only information relating to layer reflectivity’s, while the function F(d_i – d_n+1) contains the information related to the layer distances.

2.2.4. Demonstrating the forwards model functionality

We have demonstrated the use of the forwards model in characterising various OCT optical delay lines (ODLs) as well as typical OCT light sources.

2.2.4.1. Optical delay line simulations

The simulation produced A-scans of a typical moving optical delay line (ODL), as well as a stepped stationary ODL. For the former, the model generated typical A-scans using a Gaussian spectral light source. However, the model showed that there was not enough definition of layer peaks in the A-scan for the stationary stepped ODL. To overcome this, the model showed that the light in the ODL need to be modulated in the light source axis by an amount equal to the source wavelength (Jansz et al 2011).

2.2.4.2. Characterising simulated light sources

The typical OCT light source is a super luminescent light emitting diode (SLD). It is preferred as its spectral shape is Gaussian. The inverse Fourier transform of a Gaussian spectrum in the frequency domain is a single peak Gaussian in the time domain, ideal for multi layer interferometry. Simulated single and multiple SLD light sources were used to characterise A-scans of virtual samples. They demonstrated particular artefacts such as side lobes, whose intensity decreased as the sources central wavelengths moved closer together. Equation (1) was also verified by the simulation demonstrating that broader bandwidth sources generated thinner A-scan peaks (Jansz et al 2012).

2.3. The backward fitting model (BFM) - genetic algorithm approach

We will now propose a method for determining the reflectivity’s and distances $${\left\{\left(r_i,d_i\right)\right\}}_{i=1}^n$$ of each layer in a sample structure, given a discrete set of M observations of intensity at known reference arm distances, $$O={\left\{\left(I^j,d_{n+1}{}^j\right)\right\}}_{j=1}^M$$. It will be assumed that the intensity observations will be made at local maxima of I ; at least one observation is made for d_n+1 such that $$\left|d_i-d_{n+1}\right|>1.5\Delta \lambda$$ for all i = 1,2,…,n; the reference arm distances are indexed such that $$d_{n+1}^j<d_{n+1}^{j+1}$$ for j = 1,2,…,M – 1; and the sample layers are separated by at least Δλ to ensure minimal layer interference effects.

2.3.1. BFM – genetic algorithm method

The problem of reverse fitting an interferogram to a set of M observations $$O={\left\{\left(I^j,d_{n+1}^j\right)\right\}}_{j=1}^M$$ can be formulated as a least squares problem:

$$\underset{s\in S}{\min }LS(s):=\underset{s\in S}{\min }{\displaystyle \sum _{j=1}^M{\left(I^j-B_0-{\displaystyle \sum _{i=1}^n\left|B_i\right|\tilde{F}\left(d_i-d_{n+1}{}^j\right)}\right)}^2}$$E28

where the solution space $$S=\left\{s=\left(B_0,\left|B_1\right|,\dots ,\left|B_n\right|,d_1,d_2,\dots ,d_n\right):s\in \Re \times {(0,0.5)}^n\times {\Re }^n\right\}$$, and the function $$\tilde{F}(x)$$ is the cubic-spline interpolant of the local maxima of F(x). Since this is a non-linear optimization problem there is no obvious systematic/efficient approach to obtaining a solution, so instead we will utilize a meta-heuristic genetic algorithm approach to try to identify a ‘near optimal’ solution.

The basic idea behind a genetic algorithm is as follows (Hillier and Lieberman 2005):

1. Define an initial population $$S^0$$of solutions, and evaluate the objective function (equation 27) for each member of the population.

2. Select out the less optimal solutions from population $$S^0$$ (i.e. those which have a large objective value). Randomly ‘breed’ the remaining solutions to obtain a new solution population $$S^1$$. The optimal solutions from which the new solutions are bred are called ‘parents’, while the new solutions are called ‘children’.

3. Repeat step 2 until a pre-defined termination criterion is satisfied.

By continually selecting out the optimal solutions and ‘breeding’ these to give more optimal solutions, the selected solution set evolves towards a near optimal solution i.e. it evolves towards the specific sample parameter values producing the interferogram, produced by the forwards model. One testable assertion in this mathematically stochastic paradygm, is that, for an optimal solution to be achieved, diversity needs to be maintained in the population at all iterations. This may be achieved in a number of ways. In this approach, two methods will be employed:

1. By introducing a mutation rate in which child solutions randomly inherit a feature not possessed by either parent.

2. The independent addition of new members which were not bred from the initial population.

2.3.2. BFM – The genetic algorithm

Initialise: Specify values for the following parameters:

• The size of the solution population P.

• The non-negative integers P₁, P₂, P₃, which respectively define the number of best solutions retained for the next population, the number of solutions randomly selected from the worst solutions for inclusion in the next population, and the number of new solutions bred from the solutions retained from the current population. Note that P₁, P₂, P₃, must be chosen such that P₁+ P₂+ P₃ ≤ P.

• The maximum number of algorithm cycles or generations, N.

• The termination tolerance $$\epsilon$$.

• The mutation rate MR %

Following are the steps required to perform the algorithm:

Identify all observations $$O^{\prime }$$ which are locally maximum relative to the other observations in $$O$$ (i.e. $$O^{\prime }=\left\{\left(I^j,d_{n+1}{}^j\right):I^j>I^{j-1}\text{ and }{I}^j>I^{j+1}\right\}$$). The size of $$O^{\prime }$$ (denoted $$n^{\prime }$$), is the estimated number of layers in the sample. Denote the index set of $$O^{\prime }$$ as $$J^{\prime }=\left\{i:\left(I^i,d_{n+1}{}^i\right)\in O^{\prime }\right\}$$.

Refine the solution space $$S$$ as follows:

1. Set the offset value $$B_0=\min I^j$$;

2. Note that the location of the layer identified in Step 1 will be such that $$d_{i-1}<d_i<d_{i+1}$$for all $$i\in J^{\prime }$$.

Randomly generate an initial population $$S^0\subset S^{\prime }$$ of size P, and evaluate LS(s) for every $$s\in S^0$$.

Use the existing population $$S^k$$($$k\in \left\{0,1,2,\dots \right\}$$) to generate $$S^{k+1}$$ as follows:

1. $$S_1{}^{k+1}$$$$P_1$$$$s\in S^k$$LS(s)
2. $$S_2{}^{k+1}$$$$P_2$$$$S^k\backslash S_1{}^{k+1}$$

3. Generate /breed a set $$S_3{}^{k+1}$$ of $$P_3$$ new solutions. Each element $$s\in S_3{}^{k+1}$$ is constructed as a convex combination of randomly chosen solutions $$s_1,s_2\in S_1{}^{k+1}\cup S_2{}^{k+1}$$

$$s=\alpha s_1+\left(1-\alpha \right)s_2$$E29

where $$\alpha \in \left(0,1\right)$$ is a uniform random variable.

A mutation is randomly applied to MR % of the elements of the solution vector S.

The mutated element is set to a random value, while ensuring that $$s\in S$$.

1. Randomly generate/introduce a set $$S_4{}^{k+1}$$ of $$P-{\displaystyle \sum _{i=1}^3{P}_i}$$ new solutions.

2. Set

   $$S^{k+1}={\displaystyle \underset{i=1}{\overset{4}{\cup }}{S}_i{}^{k+1}}$$E30
   .

Repeat Step 4 for a fixed number of cycles N, or until there exists a member $$s\in S^k$$ such that $$LS\left(s\right)<\epsilon$$, for a pre-specified tolerance $$\epsilon$$.

The optimal member of the final population $$s=\left(B_0,\left|B_1\right|,\dots ,\left|B_{n^{\prime }}\right|,d_1,d_2,\dots ,d_{n^{\prime }}\right)$$ directly provides the total distance travelled by the laser light off each of the $$n^{\prime }$$ layers identified by the algorithm (i.e. $$d_1,d_2,\dots ,d_{n^{\prime }}$$). The reflectivities of the layers can be obtained recursively as follows:

$$r_1={\left(\frac{2\left|B_1\right|}{R{F}_{n+1}}\right)}^2$$E31

$$r_i={\left(\frac{2\left|B_i\right|}{R{F}_{n+1}{\displaystyle \prod _{j=1}^{i-1}\left(1-r_j\right)}}\right)}^2,i=2,3,\dots ,n^{\prime }$$E32

Note: We will never get negative values of the coefficients because we are using the upper envelope function.

3. Method

3.1. Testing the precision of the genetic algorithm model (GAM)

This characterization of the efficacy of the BFM – Genetic Algorithm compares the spread of the depth and reflectivity of each layer for 2, 3, 4 and 5 layer sample models, for 20, 50, 100, 200, 400 and 800 generations of the GAM. This is in order to gauge the level of precision of the GAM as a function of generation number and number of sample strata. Ideally the precision of the two parameters should remain statistically not significantly different from the actual parameter values for a given layer irrespective of the number of layers in the sample.

The following GAM control variables were used for 2, 3, 4 and 5 strata virtual samples:

• P = 1000 – The total solution population size at the start of each generation.

• P1 = 100 – The number of optimal solutions retained for inclusion in the next generation.

• P2 = 250 – The number of solutions randomly selected from the remaining ranked solutions not in P1, for inclusion in the next generation.

• P3 = 300 – The number of new solutions (children) bred from a convex combination of P1 and P2, for inclusion in the next generation.

• P4 = P – (P1+P2+P3) = 350 – number of extra solutions randomly generated so that the total population quota, P, is achieved for the next generation. This is an implicit variable.

• N= 20, 50, 100, 200, 400, 800.

• Epsilon (ε) = 0.000001 – termination tolerance.

• Mutation rate = 0.01 – for each generation of size P, 1% of P has an introduced mutation.

The thickness was equivalent for each stratum being 100μm. The strata reflectivities used were based on approximate tissue strata refractive indices (n), alternating between 1.45 for the first stratum and 1.49, with n of air being 1, above the first stratum. The expected reflectivities for the i^th interface can then be calculated from,

$$R_i=\frac{{(n_i-n_{i-1})}^2}{{(n_i+n_{i-1})}^2}$$E33

,

given that i = 1 refers to the top most stratum (n₁) interface with the air (n₀). The calculated percentage reflectivities are then 3.37359 for the surface and only 0.0185108 for the lower interfaces. Each generation number is repeated 7 times.

The degree of statistical “equality” of simulated depths to actual strata depths follows from the relative error of the depths introduced by the degree of coherence of the virtual light source. It provides a region of acceptable depths about the actual depths that are experimentally not significantly different from their actual values. The degree of source coherence defines an axial resolution that is measured by the source’s coherence length (L_c) divided by the tissue refractive index. With the average tissue refractive index of 1.47 and the light source’s central wavelength (λ₀) of 1550 nm and the full-width-half-maximum bandwidth (Δλ) of 40 nm, the L_c is 26.5 μm (equation 1). Hence, the in-tissue resolution is 18.03 μm, which is 9.01 μm either side of the actual depth. It follows then, that the relative error (%) for each strata of depth: 100, 200, 300, 400, 500 μm is 9.01, 4.51, 3.00, 2.25, 1.80%, respectively, using the virtual source’s characteristics.

The precision results were compared graphically, with the independent variable as the generation number and the dependent variable a measure of the precision, i.e. the relative 99% confidence interval of the depth and reflectivity, with respect to the average depth and reflectivity, expressed as a percentage. This relative percentage precision was compared as one graph per layer that included each of the two to five layer sample structures. As indicated above, the relative error (%) associated with the source coherence for each layer was included in each of these graphs as an indication of acceptable precision of the simulated depths.

3.2. GAM parameter effect on accuracy and precision

This characterization of the efficacy of the BFM – Genetic Algorithm identifies a trend for which combination of GAM control variable magnitudes achieve optimal convergence speed, i.e. the least number of generations that give a satisfactory precision level for the returned parameters, strata depths and reflectivities. The procedure uses the same five strata sample as in the GAM precision section above (3.1) with the following variation of GAM control parameters in table 1.

Three trials of 7 simulations, each for only 20 generations were undertaken. The pooled mean and pooled standard deviation of the 3 trials for depth and reflectivity were calculated. The relative error (%) between these means and the actual depths and reflectivities were compared graphically for each parameter variation. Also compared graphically, were the relative error (%) between the pooled standard deviations and their pooled means, for each depth and its reflectivity. As indicated above, the relative error (%) associated with the source coherence for each layer was included in these graphs as an indication of acceptable accuracy of the simulated depths.

GAM; Parameters	Varying P1	Varying P2	Varying P3	Varying Mutation Rate (MR)
P	1000	1000	1000	1000
P1	0, 20,... 800	100	100	400
P2	100	0, 20,... 800	100	100
P3	100	100	0, 20,... 800	100
Mutation Rate	0.01	0.01	0.01	0.5, 0.2, 0.1, 0.05, 0.02, 0.01, 0.005, 0.002, 0.001, 0.0005, 0.0002, 0.0001
ε	0.0000001	0.0000001	0.0000001	0.0000001

Table 1.

Variation of GAM parameters

4. Results, Inferences and Implications

4.1. GAM precision

The following section demonstrates the capability of the GAM to provide a precision for layer depth and reflectivity for samples with different number of layers, such that the parameter precisions are not significantly different for any given layer. Results for strata depth precisions is followed by the reflectivity precision results, each demonstrated graphically.

4.1.1. GAM strata depth precision

Fig. 2 shows that, though the GAM precision for each layer depth is dependent on the number of sample layers, all precisions fall within the relative error boundary of the source resolution, i.e. they are all acceptable precision variations. Layer depth precision generally decreases with the increase in the number of layers. The first, more reflective layer, has better precision by more than an order of magnitude, compared to the other layers. This is due to the more prominent peak of the first layer due to its greater reflectivity.

Figure 2.

The precisions (99% confidence intervals) of the GAM calculated strata depths by number of sample layers and number of GAM generations. The precisions for returned (A) layer 1, (B) layer 2, (C) layer 3, and (D) layer 4 depths, are presented.

4.1.1. GAM layer reflectivity precision

The reflectivity precision parallels the depth precision trends for each layer (Fig. 3 A – D), though these precisions are significantly less sparse.

4.2. GAM control parameter effects on performance

Here we investigate the effect of the four GAM parameters on the accuracy and precision of the calculated layer depths and reflectivities, produced from 20 generations, for a five layer virtual sample. Hence the better the accuracy and precision, the more optimal the GAM parameter. The accuracy was measured by the relative error (%) of the resulting pooled mean depth and pooled mean reflectivity compared to the actual, for three trials each of seven 20 generation GAM cycles. The precision was denoted by the relative standard deviation, being the quotient of the pooled standard deviation to the pooled mean for the three trials, expressed as a percentage.

4.2.1. GAM strata depth accuracy and precision

It is important to note, that due to time constraints, the GAM accuracy and precision to predict the layer depth has been characterised using a univariate approach.

Figure 3.

The precisions (99% confidence intervals) of the GAM calculated strata reflectivities by number of sample layers and number of GAM generations. Precision for returned (A) layer 1, (B) layer 2, (C) layer 3, and (D) layer 4 reflectivities, is presented.

4.2.1.1. P1 effect on layer depth

First P1 is the control parameter that indicates the number of best solutions kept for the next generation. It would then seem probable that maximizing this value would arrive at the solution sooner. From the GAM results of three trials, Fig. 4 contradicts this prediction.

Fig. 4 shows that optimal accuracy and precision is not limited to maximizing P1 but rather that in layer 1, optimality is achieved at lower values of P1, while all other layers are optimized randomly across the P1 spectrum. Fig. 4 shows that the percentage of depth relative errors (DREs) that fall inside the percentage error boundaries of the source’s in-tissue resolution for layers 1 to 5 (Fig. 4A to 4E) are 100, 61, 32, 73 and 100% respectively.

Considering layer 1 (Fig. 4A) for the given fixed GAM parameters, for optimal layer 1 depth predictions, P1 needs to be between 100 and 300, with best results below 200. Note also that in this region, the precision is also the greatest, being closest to the null line.

Figure 4.

Effect of P1 varying on the relative errors and standard deviation errors for (A) layer 1, (B) layer 2, (C) layer 3, (D) layer 4, (E) layer 5 depth GAM calculations. These are the Average values of 3 trials. The other GAM control variables, P2, P3, ε and Mutation rate were fixed at 100, 100, 0.0000001, and 0.01, respectively.

For layer 2 (Fig. 4B), P1 delivers best results 61% sporadically across the 0 to 800 spectrum. A constant dip in the precision (red squares) between 120 and 200 may suggest an optimal region of P1 for predicting layer 2. Also P1 from 340 to 580 shows significant accuracy and precision, with the depth relative error within the relative error boundaries of the source’s in-tissue resolution. Layer 2 depth is consistently being underestimated, indicated by the negative relative error values (blue diamonds).

Layer 3 has best depth prediction at P1 = 180, while the accuracy and precision values are not separated. For the accuracy to be not significantly different to zero percent, the precision values (red squares) need to be significantly more positive than the absolute value of the accuracy (relative error) values (blue diamonds).

Layer 4 has best results with P1 between 40 and 200, and 600 to 800. Layer 5 has best results with P1 between 120 and 220.

Layer 2 results show the greatest degree of spread with in a corridor, basin or boundary of attraction. This ‘basin of attraction’ thins as the layer becomes deeper as well as getting progressively more accurate and precise. This difficulty of the GAM to locate the second and third layers is because they falls in the shadow of the much more intense (i.e. more reflective) first layer, while locating the lower layers becomes progressively easier the deeper the GAM searches.

Though the GAM results seem random, their ‘basin of attraction’ localization imparts a degree of confidence in the reproducibity of the GAM solutions which are reliably attracted to some random data corridor or geometry; for layer 1 a ‘wine glass’ geometry, layer 2 a wide corridor, and layers 3 – 5, narrow corridors. These regions of ‘random’ solutions are a product of the interplay between the GAM’s stochasticity and its parametric control.

In summary, for the constant values of P2, P3, Mutation rate and tolerance being 100, 100, 0.01 and 0.0000001 respectively, while P1 is best between 100 and 300 for layer 1, there is no particular value of P1 that is best for the other layers. More refined analysis investigating other values of P2, P3, and MR, for this range of P1, is necessary to establish a broader picture of the effects of P1 and its interaction with the other GAM control parameters.

4.2.1.2. P2 effect on layer depth

P2 is the control parameter that indicates the number of not-best solutions randomly chosen from the sample set remaining after the P1 set has been chosen, which are kept for the next generation. It would then seem probable that minimizing this value would arrive at the solution sooner. From the GAM results of three trials, Fig. 5 suggests this prediction for the layer depth calculations.

Figure 5.

Effect of P2 varying on the relative errors and standard deviation errors for (A) layer 1, (B) layer 2, (C) layer 3, (D) layer 4, (E) layer 5 depth GAM calculations. These are the average results of 3 trials. The other GAM control variables, P1, P3, ε and Mutation rate were fixed at 100, 100, 0.0000001, and 0.01, respectively.

Fig. 5 shows that the percentage of depth relative errors (DREs) that fall inside the percentage error boundaries of the source’s in-tissue resolution for layers 1 to 5 (Fig. 5A to 5E) are 90, 44, 20, 83 and 100% respectively.

In Fig. 5A, the layer 1 accuracy (blue Diamonds) and precision (Red squares) are best for P2 between 0 and 100, being closer to zero. The precision reduces and accuracy is more random beyond these points. Layer 2 results show less spread of the accuracy and precision than for layer 1. However because the blue error boundary of the source resolution is now half compared to layer 1, ~50% less values fall in the +/- 4.51% acceptability zone. As for the previous section, a ‘basin of attraction’ for both the accuracy and precision, thins as the layer becomes deeper as well as getting progressively more accurate and precise. Again the difficulty of locating the second layer is noted due to the magnitude of the first layer’s signal. This also affects the third layer, for which the accuracy values are now larger than the precision, showing that 34% of the depth estimates are significantly different to the actual layer depth value. This shadowing effect is reduced for the lower layers where the majority of the values are not significantly different to the actual depths, with progressively becoming more accurate and precise.

In summary, for the constant values of P1, P3, Mutation rate and tolerance being 100, 100, 0.01 and 0.0000001 respectively, while P2 is best between 0 and 100 for layer 1, there is no particular value of P2 that is best for the other layers. More refined analysis investigating other values of P1, P3, and MR, for this range of P2, is necessary to establish a more accurate picture of the effects of P2 and its interaction with the other GAM control parameters.

4.2.1.3. P3 effect on depth

P3 is the control parameter that indicates the number of children generated by a randomly selected, with replacement, convex combination of P1 and P2. It would then seem probable that neither maximizing or minimizing this value would arrive at a suitable solution sooner. From the GAM results of three trials, Fig. 6 suggests there is no value of P3 that is reproducibly of high accuracy or precision.

Fig. 6 shows that the percentage of depth relative errors (DREs) that fall inside the percentage error boundaries of the source’s in-tissue resolution for layers 1 to 5 (Fig. 5A to 5E) are 100, 63, 39, 88 and 100% respectively.

In Fig. 6A, the layer 1 accuracy values (blue Diamonds) are not significantly different to the actual depth. Layer 2 results (Fig. 6B) show more spread of the accuracy and precision than for layer 1, with the depth being under estimated for all values of P3. As for the previous section, a ‘basin of attraction’ for both the accuracy and precision, thins as the layer becomes deeper as well as getting progressively more accurate and precise. Again due to the magnitude of the first layer’s signal, The GAM has difficulty of locating the second and third layers. For the third layer (Fig. 6C) the accuracy values are now larger than the precision, showing that 22% of the depth estimates are significantly different to the actual value. This shadowing effect is reduced for layers 4 (Fig. 6D) and 5 (Fig. 6E); the majority of the values are not significantly different to the actual depths, though layer 3 and 4 values are over estimated. Layer 5’s accuracy is demonstrated by symmetry about the zero line, with precision comparable to layer 1.

Figure 6.

Effect of P3 varying on the relative errors and standard deviation errors for (A) layer 1, (B) layer 2, (C) layer 3, (D) layer 4, (E) layer 5 depth GAM calculations. These are the average results of 3 trials. The other GAM control variables, P1, P3, ε and Mutation rate were fixed at 100, 100, 0.0000001, and 0.01, respectively.

In summary, for the constant values of P1, P2, Mutation rate and tolerance being 100, 100, 0.01 and 0.0000001 respectively there is no particular value of P3 that will benefit the GAM accuracy or precision. Further, more refined analysis investigating other values of P1, P2, and MR, for this range of P3, is necessary to establish a broader picture of the effects of P3 and its interaction with the other GAM control parameters.

4.2.1.4. Mutation rate effect on layer depth

Mutation rate (MR) is the control parameter that indicates the proportion of each generation of size P, that has (MR x 100) % introduced mutation.

Fig. 7 shows that the percentage of depth relative errors (DREs) that fall inside the percentage error boundaries of the source’s in-tissue resolution for layers 1 to 5 (Fig. 5A to 5E) are 100, 50, 0, 50 and 100% respectively.

In Fig. 7A, the layer 1 accuracy values (blue Diamonds) are not significantly different to the actual depth, with a decreasing general trend towards smaller MR, except for a MR of 0.005, which is as accurate as 0.5. Layer 2 relative error results (Fig. 7B) are all under estimated but not significantly different to zero. Again due to the magnitude of the first layer’s signal, The GAM has difficulty of locating the second and third layers. For the third layer (Fig. 7C) all the accuracy values fall outside the relative error boundary of the source resolution as well as being over estimated. Also 17% of the values are larger than the precision, implying that these depth estimates are significantly different to the actual value. This shadowing effect is reduced for layers 4 (Fig. 7D) and 5 (Fig. 7E); the majority of the values are not significantly different to the actual depths. Layer 5’s accuracy is demonstrated by symmetry about the zero line, with significantly better accuracy and precision than layer 1. As for the previous section, a ‘basin of attraction’ for both the accuracy and precision, widens as the layer becomes deeper as well as getting progressively more accurate and precise.

Figure 7.

Effect of Mutation rate varying on the relative errors and standard deviation errors for (A) layer 1, (B) layer 2, (C) layer 3, (D) layer 4, (E) layer 5 depth GAM calculations. These are the average results of 3 trials. The other GAM control variables, P1, P2, P3 and ε were fixed at 400, 100, 100, and 0.0000001, respectively.

In summary, for the constant values of P1, P2, P3 and tolerance being 400, 100, 100, and 0.0000001 respectively, there is no particular value of MR that will benefit the GAM accuracy or precision for calculating the layer depths. More refined analysis investigating other values of P1, P2, and P3, over this range of MR, is necessary to establish a broader picture of the effects of MR and its interaction effects with the other GAM control parameters on the estimation of layer depth.

4.2.2. GAM layer reflectivity accuracy and precision

It is important to note, that due to time constraints, the GAM accuracy and precision to predict the layer reflectivities has been characterised in only a univariate manner, with each parameter varied singularly, while the other parameters were kept constant, as indicated.

4.2.2.1. P1 effect on layer reflectivity

First P1 is the control parameter that indicates the number of best solutions kept for the next generation. It would then seem probable that by maximizing this parameter, for each generation, the GAM will converge faster to the actual layer reflectivities.

Figure 8.

Effect of P1 varying on the relative errors and standard deviation errors for (A) layer 1, (B) layer 2 reflectivity GAM calculations. These are the Average values of 3 trials. The other GAM control variables, P2, P3, ε and Mutation rate were fixed at 100, 100, 0.0000001, and 0.01, respectively.

However, the two Fig. 8 graphs demonstrate no relationship between P1 and the accuracy and precision of the calculated reflectivity for varying P1 from 0 to 800, given that the other fixed GAM control parameters, P2, P3, MR and ε, are 100, 100, 0.01 and 0.0000001 respectively. Layer 1 (Fig. 8A) reflectivity is the most accurate and precise, with the other layers demonstrating similar accuracy and precision to that of layer 2 (Fig. 8B). However as the pooled standard deviation precision values (red squares) are less than the absolute relative error values (blue diamonds), i.e. accuracy, for all layers, all reflectivity estimates are significantly different from the actual layer reflectivities. Hence, 20 generations is too small to obtain satisfactory accuracy and precision, even after 3 trials, using this GAM version. Also the level of precision is not satisfactory for layers 3 to 5, being each similar to layer 2 (Fig. 8B).

More refined analysis investigating other values of P2, P3, and MR over this range of P1, is necessary to establish a broader picture of the effects of P1 and its interaction effects with the other GAM control parameters on the estimation of layer reflectivity.

4.2.2.2. P2 effect on layer reflectivity

P2 is the control parameter that indicates the number of not-best solutions randomly chosen from the sample set remaining after the P1 set has been chosen, which are kept for the next generation. It would then seem probable that minimizing this value would arrive at the solution sooner. From the GAM results of three trials, Fig. 9 suggests this prediction for the layer reflectivity calculations, but only for P2 = 0, for all layers. Layer 3 to 5 results are similar to layer 2 results (Fig. 9B).

Figure 9.

Effect of P2 varying on the relative errors and standard deviation errors for (A) layer 1, (B) layer 2 reflectivity GAM calculations. These are the Average values of 3 trials. The other GAM control variables, P1, P3, ε and Mutation rate were fixed at 100, 100, 0.0000001, and 0.01, respectively.

Other than P2 = 0, the two Fig. 9 graphs demonstrate no general relationship between P2 and the accuracy and precision of the calculated reflectivity for varying P1 from 0 to 800, given that the other fixed GAM control parameters, P1, P3, MR and ε, are 100, 100, 0.01 and 0.0000001 respectively. Again Layer 1 (Fig. 9A) reflectivity is the most accurate and precise, with the other layers, 3 to 5, demonstrating similar accuracy and precision to that of layer 2 (Fig. 8B). However as the pooled standard deviation precision values (red squares) are less than the absolute relative error values (blue diamonds), i.e. accuracy, for all layers, all reflectivity estimates are significantly different to the actual values. Hence, 20 generations is too low, to obtain satisfactory accuracy and precision, even after 3 trials, with this GAM version.

Figure 10.

Effect of P3 varying on the relative errors and standard deviation errors for (A) layer 1, (B) layer 2 reflectivity GAM calculations. These are the Average values of 3 trials. The other GAM control variables, P1, P2, ε and Mutation rate were fixed at 100, 100, 0.0000001, and 0.01, respectively.

4.2.2.3. P3 effect on layer reflectivity

P3 is the control parameter that indicates the number of children generated by a randomly selected, with replacement, convex combination of P1 and P2. It would then seem probable that neither maximizing or minimizing this value would arrive at a suitable solution sooner. From the GAM results of three trials, Fig. 10 suggests there is no value of P3 that is reproducibly of high accuracy or precision for the reflectivity, and all values are significantly different to the actual reflectivities for layers 1 and 2. This also applies to the lower layers as they are similar to layer 2 results (Fig. 10B).

For the same reason that the depth results were best for layer 1, the reflectivity results are best for layer 1 because its reflectivity is orders of magnitude greater than the lower layer reflectivities. Layer 2 results show no relationship at all to P3, with the accuracy values significantly different to 0 %.

4.2.2.4. Mutation rate effect on layer reflectivity

Mutation rate (MR) is the control parameter that indicates the proportion of each generation of size P, that has (MR x 100) % introduced mutation.

Figure 11.

Effect of Mutation rate varying on the relative errors and standard deviation errors for (A) layer 1, (B) layer 2 reflectivity GAM calculations. These are the average results of 3 trials. The other GAM control variables, P1, P2, P3 and ε were fixed at 400, 100, 100, and 0.0000001, respectively.

In Fig. 11A, the layer 1 reflectivity accuracy relative error values (blue Diamonds) are significantly different to zero percent. Layer 2 relative errors (Fig. 11B) are all over estimated and significantly different to zero percent, based on the large precision values (red squares). This is also true for the lower layers, 3 to 5. Also the same non-relationship trend of accuracy and precision for layers 3 to 5, is similar to layer 2. Again due to the large first layer’s reflectivity and very small reflectivity of the lower layers, this version of the GAM has difficulty of locating the lower layer reflectivities, given only 20 generations of evolution.

In summary, for the constant values of P1, P2, P3 and tolerance being 400, 100, 100, and 0.0000001 respectively, there is no particular value of MR that will benefit the GAM accuracy or precision for calculating the layer reflectivities. More refined analysis investigating other values of P1, P2, and P3, over this range of MR, is necessary to establish a broader picture of the effects of MR and its interaction effects with the other GAM control parameters on the estimation of layer reflectivity.

4.2.3. GAM layer depth and reflectivity accuracy and precision summary

Table 2 summarises the results of the GAM estimates of depth and reflectivity. With the exception of P2, all other parameters show no significant relationship to the relative errors of the GAM estimated depth and reflectivity calculated.

(1) The percentage of results that are within the relative error boundary of the source resolution. NSV = no specific value; AV = all values; SD = significantly different.

Layer #; (P = 1000)	Optimal Depth				Optimal Reflectivity
P1; 0 to 800	P2; 0 to 800	P3; 0 to 800	MR; 0.0001 - 0.5	P1; 0 to 800	P2; 0 to 800	P3; 0 to 800	MR; 0.0001 - 0.5
1	100 - 300	0-100	AV	AV	All SD; to 0%	0	All SD; to 0%	All SD; to 0%
2	NSV- 61%^(1)	NSV -44%	NSV -; 63%	NSV -; 50%	All SD; to 0%	0	All SD; to 0%	All SD; to 0%
3	NSV - 32%	NSV -; 20%	NSV -; 39%	NSV -; 0%	All SD; to 0%	0	All SD; to 0%	All SD; to 0%
4	NSV -; 73%	NSV -; 83%	NSV -; 88%	NSV -; 50%	All SD; to 0%	0	All SD; to 0%	All SD; to 0%
5	AV	AV	AV	AV	All SD; to 0%	0	All SD; to 0%	All SD; to 0%

Table 2.

General and relative error results summary

In the case of P2 = 0, the reflectivity accuracy and precision are significantly improved for all layers, while the layer depth accuracy and precision are significantly improved for the first layer. Except for P2 = 0, all other values of P1, P2, P3, and MR return reflectivities that are significantly different to the actual reflectivities as indicated in the Table 2 as “All SD to 0%”, meaning that all the values of the control parameter give reflectivities that are significantly different to the actual reflectivity values.

5. Discussion

5.1. GAM precision

The results have demonstrated that the precision of GAM layer depth estimation worsens for any given layer as the number of layers in the sample increases. This difference in results between different layered samples reduces as the GAM generation number increases. One can now gauge the degree of result reliability with a given A-scan, which indicates number of layers, so that an appropriately long generation number calculation can be undertaken if a sample with a larger number of layers is predicted.

Though the layer depth and reflectivity precision increases with increasing number of generations, so does the time to arrive at a more optimal solution. This increase in time is linear: if one generation takes 1 min, then 800 generations will consume 800 min. Certainly not practical for returning specific layer depths from a B-scan consisting of 600 A-scans, if each A-scan result was generated from 800 GAM generations. Clearly it is necessary to optimize the speed by minimizing the generation number to the point that the relative error of the depth is just within the light source’s resolution boundary. Furthermore, The back fitting software’s application will be assisting the clinician, by returning depth and reflectivity of points of interest on A-scans in the B-scans, not necessarily all depths and reflectivities of the B-scan.

The reflectivity precision parallels the depth precision trends for each layer (Fig. 3 A – D) because locating the correct depth is enhanced by greater layer reflectivity. Also, the more the layers the more there are peaks in the A-scan for the GAM to locate, requiring more generations to achieve equivalent precision to A-scans with less peaks.

5.2. GAM control parameter effects on depth estimation

For all the parameters, P1, P2, P3 and Mutation rate, layer 3 has least accurate results with most or all values falling outside the relative error boundary of the expected source resolution. layer 2 results show the greatest degree of spread with in a corridor, basin or boundary of attraction. This ‘basin of attraction’ thins as the layer becomes deeper as well as getting progressively more accurate and precise. This difficulty of the GAM to locate the second and third layers is because it falls in the shadow of the much larger (i.e. more reflective) first layer, while locating the lower layers becomes progressively easier the deeper the GAM searches.

The term ‘basin of attraction’ is used because, for the three trials, each trial produced the same region of graph values without the actual values being the same. This is a product of the interplay between the GAM’s stochasticity and its parametric control. This is important to note as it imparts a degree of confidence in the reproducibity of the GAM solutions which are reliably attracted to some random data corridor or geometry. This ‘basin of attraction’ construct is applicable and evident in all of the GAM depth estimates, because the GAM is a parameter controlled stochastic algorithm, and not solely stochastic.

5.3. GAM control parameter effects on reflectivity estimation

For all the parameters, P1, P2, P3 and Mutation rate, layer 1 is the most accurate and precise. However all GAM estimated reflectivities are significantly different to the actual values for all five layers. A larger number of evolution generations are needed, using this version of the GAM. This will require either repeating the GAM evolution to run for a much larger number of generations or extrinsically evolving the GAM by altering the software so that it evolves faster, i.e. alter the existing GAM version so that it arrives at suitable accuracy and precision sooner.

5.4. Future backwards fitting model research

Further trials with realistic non-Gaussian light sources (Rossetti et al 2005), (Shidlovski 2008), (Unterhuber et al 2008) that generated A-scans with side lobes and other artefacts, will require the programming of the GAM to detect and disregard such false information, resulting in the return of accurate parameter values. This will involve further extrinsic evolution of the GAM as well as exploring other types of Backwards Fitting Models that though more deterministic may be much faster due to loss of the random element. However further investigation is needed to establish this.

6. Conclusion

This chapter has demonstrated that a genetic algorithm approach suggests some success at extracting layer depth and reflectivity values from virtual samples with realistic biological layer depths and refractive indices. The GAM’s efficacy has been challenged by these realistically small reflectivities, returning values that are an order of magnitude larger than the actual values. More extensive testing of the GAM using a multivariate approach is needed to fine tune the GAM control parameters to deliver accurate and precise solutions at optimal speed for the least number of generations possible. Further extrinsic evolution of the GAM to extend its application to non Gaussian light sources is envisaged. This will enable the deconvolution of more realistic A-scans. This may provide the opportunity to give a more objective measure of sample layer characteristics to better inform physicians, assisting in their diagnosis and monitoring of particular tissue pathologies.