Method and Apparatus for Multi-spectral Imaging and Analysis of Skin Lesions and Biological Tissues

[0001] This application claims the benefit of U.S. Provisional Patent Application No. 61/088,170, filed Aug. 12, 2008, the entirety of which is incorporated herein by reference.

[0002] This invention relates to methods and apparatus for imaging and analysis of skin lesions and biological tissues.

[0003] Skin cancer is a significant health problem in the United States. It has been reported that one of five Americans will get some form of skin cancer in their lifetime. Currently, nearly half of new cancers are diagnosed as skin cancers. Malignant melanoma, the most fatal skin cancer, first forms at the upper layers of skin. When metastasized, cancerous cells from melanoma enter blood vessels and proliferate throughout the body. Malignant melanoma is highly fatal if not detected in early stages. However, it can be cured with nearly 100% survival rate if removed at an early stage.

[0004] Physicians usually use the “ABCD” rule to determine if a lesion under investigation is malignant melanoma. The acronym “ABCD” refers to asymmetry, border, color and diameter, respectively. Malignant melanoma typically has an asymmetrical shape, an uneven border, varied colors and a large diameter. Once a suspicious lesion is excised a diagnosis can be confirmed by other instruments. However, neither visual inspection using the “ABCD” rule nor examination of the excised lesion can provide depth information of the skin cancer, which is a crucial signature to grade the degree of invasion of a skin lesion. Angiogenesis, or increased blood flow, plays a very important role in detection of melanomas in early curable stage. Specific patterns of distribution of melanin, oxy-hemoglobin and de-oxy-hemoglobin can lead to characterization of dysplastic nevi and their potential for transformation into malignant melanoma in very early phases.

[0005] Various light transportation models have been used by researchers to reconstruct information to characterize skin-lesions. For example, a Kubelka-Munk model was used to simulate the formation of images of melanoma and presented a method to recover blood and melanin distribution in various skin layers. Claridge et al., An inverse method for recovery of tissue parameters from colour images, Information Processing in Medical Imaging. Springer, Berlin, LNCS2732, pp. 306-317. However, the Kubelka-Munk model is theoretically established in a one-dimensional system with point-based measurements. For more complex geometries, Monte Carlo simulation or Diffusion Approximation has been used in optical tomographic modalities for more accurate reconstructions. The commonly adopted strategy for reconstruction involves dividing the field of view into a number of voxels and assuming constant optical properties in each voxel. The optical properties are then estimated voxel-by-voxel by matching model predicted measurements to the actual measurements. This is a typically under-determined and ill-posed inverse problem as the number of measurements is usually much less than the number of voxels to be reconstructed. In general, the forward process is a mapping from high dimensional space (unknown optical properties of voxels) to low dimensional space (limited measurements). Due to the loss of information during the forward process, the solution to the inverse problem is not unique and usually has to be stabilized through various regularization methods. It is therefore difficult to obtain a quantitatively accurate and well-localized solution. In addition, light photons are quickly diffused in a turbid medium such as human skin. As a result, there is a strong dependence or similarity between different measurements such that increasing the number of measurements would not lead to a dramatic change in the characteristic behavior of the inverse problem.

[0006] In recent years, optical medical modalities have drawn significant attention from researchers. Visible and near-infrared light wavelengths have been used in surface reflectance, transillumination and transmission based methods. See, Ganster et al., Computer aided recognition of pigmented skin lesions, Melanoma Research, vol. 7 (1997); Seidenari et al., Digital video-microscopy and image analysis with automatic classification for detection of thin melanomas, Melanoma Research 9(2), 163-171 (1999); Menzies et al., Automated instrumentation and diagnosis of invasive melanoma, Melanoma Research vol. 7, 13 (1997); Claridge et al., From color to tissue histology: Physics-based interpretation of images of pigmented skin lesion, Medical Image Analysis, pp. 489-502 (2003); Tomatis et al., Automated melanoma detection: multi-spectral imaging and neural network approach for classification, Med. Phys. 30(2), pp. 212-221 (2003); Tomatis et al., Spectro-photo-metric imaging of subcutaneous pigmented lesion: Discriminant analysis, optical properties and histological characteristics, J. Photochem. Photobiol., B 42, 32-39 (1998). U.S. Pat. No. 5,146,923 discloses a portable nevoscope which provides a noninvasive means to examine a skin lesion in situ, and provides a means to process and analyze skin lesion data relating to properties such as thickness, color, size, pigmentation, boundary, and texture. Due to the limited view and limited-angle measurements available via the prior art nevoscope, the intrinsic ill-posed and under-determined nature of optical imaging pose problems in reconstructing accurate tomographic information.

[0007] Consequently there is a need for an improved nevoscope device and methods of obtaining improved reconstruction results.

[0008] In accordance with various aspects of the present invention multispectral imaging systems and methods are provided.

[0009] Optical modalities can provide a portable imaging system for routine screening and monitoring of skin-lesions. Multi-spectral optical imaging using visible and infrared light wavelengths as disclosed herein can provide information about physiologically meaningful chromophores such as melanin, oxyhemoglobin and deoxyhemoglobin through utilization of differences in their wavelength dependent absorption and scattering coefficients. The apparatus and methods disclosed herein are generally applicable for optical image reconstruction.

[0010] In accordance with one embodiment an improved multi-spectral nevoscope is disclosed providing transillumination for imaging skin lesions for diagnosing malignant melanoma non-invasively. In one embodiment the device comprises substantially a portable optical imaging device that uses specific wavelengths in the visible and infrared spectrum of electromagnetic radiation to transilluminate a skin-lesion or a biological tissue or specimen for imaging and maps multispectral 2-dimensional images into 3-dimensional virtual space for providing 3-D distributions of pre-defined parameters representing the characteristic properties (such as melanin, hemoglobin and deoxyhemoglobin, etc.) of a skin-lesion. These characteristic distributions of specific parameters can be analyzed and used for detection and management of skin-cancers, or characterization of a biological tissue or specimen. The device allows a background transillumination source for excitation or preparation of background tissue such as the surrounding skin of a skin lesion or the entire tissue itself for fluoroscopy imaging.

[0011] In accordance with an embodiment the device may include multiple transillumination rings for background tissue preparation or excitation for lesion imaging for optimal penetration and subcutaneous illumination of skin lesions. The device may include multiple adaptive combinations of source and receiver channels distributed over the imaging area through fiber-optics cables, optical illuminators and filters, and computer-controlled image sensors such as CCD arrays. Systems in accordance with the present invention are of critical value to characterize skin lesions and biological tissues for optical and/or molecular imaging and analyses of associated distributions of characteristic parameters. A series of images obtained with multiple excitation and source-illumination geometries with multi-spectral filters may be analyzed by visual inspection/diagnosis and/or 3-D mapping of distribution of specific parameters such as oxyhemoglobin, deoxyhemoglobin and melanin for diagnostic evaluation and characterization of skin lesion or tissue.

[0012] The present apparatus may be used for clinical monitoring of skin lesions on patients with high risk of developing malignant melanoma, in addition to monitoring other skin cancers and conditions including those developed from allergic reactions in response to drugs, foods and the like.

[0013] In accordance with a further embodiment a shape-based multi-constrained reconstruction algorithm is disclosed which uses genetic algorithm-based optimization methods to find the best possible reconstruction solution. In one embodiment, a skin lesion such as melanoma is modeled as melanin and blood parts, which are delineated by two cubic tensor-product B-spline surfaces. This reduces the number of unknowns to a few control parameters of the surfaces. The parameters are then coded into a genetic algorithm to find a solution through global optimization.

[0014] In accordance with a further embodiment a multispectral imaging (MSI) method uses plural selected visible wavelengths for transillumination to acquire multiple remittance images. Different wavelengths of light are projected through fiber-optics-directed ring-light sources for transilluminating the skin lesion through the surrounding skin area for background imaging for calibration, tissue excitation, or tissue preparation for imaging. The entire remittance signature spectrum using multiple light wavelengths improves characterization of skin, dysplastic nevi, melanomas and other skin-lesions. Multiple discrete sources used in sequential imaging of the skin lesion provide extended image data for use in characterizing the skin lesion. This characterization may be based on the visual examination of multispectral transillumination images and/or computer based analysis and three-dimensional reconstruction of the skin lesion.

[0015] In accordance with a further embodiment an algorithm is provided in which a skin lesion such as melanoma is modeled as melanin, hemoglobin and deoxyhemoglobin.

[0034] It should be noted that the appended drawings illustrate only typical embodiments of this invention and are therefore not to be construed as limiting of its scope, for the invention may admit to other equally effective embodiments. Where possible, identical reference numerals have been inserted in the figures to denote identical elements.

[0035] In the following description, for purposes of explanation, specific numbers, materials and configurations are set forth in order to provide a thorough understanding of the invention. It will be apparent, however, to one having ordinary skill in the art that the invention may be practiced without these specific details. In some instances, well-known features may be omitted or simplified so as not to obscure the present invention. Furthermore, reference in the specification to phrases such as “one embodiment” or “an embodiment” means that a particular feature, structure or characteristic described in connection with the embodiment is included in at least one embodiment of the invention. The appearances of phrases such as “in one embodiment” in various places in the specification are not necessarily all referring to the same embodiment.

[0036] Now referring to FIG. 1 a schematic of a prior art Nevoscope apparatus 2 employing white light-based transillumination is shown including an optical lens 10, lens mount bracket 20, surface light fibers 30, transillumination fibers 40, a ring light 50 and monocoil 60. As shown the nevoscope is positioned over skin 100.

[0037] Now referring to FIG. 2, an embodiment of an improved Nevoscope apparatus 200 is shown. Apparatus 200 includes a housing 210, face plate 220, transillumination outer ring fiber cable 240, main imaging area fiber cables 250 and 252, focusing lens/polarizer 260, image sensor CCD 270, source mask 280 and multispectral filter 290. Source mask 280 and multispectral filter 290 may be rotatable and/or selectable, such as by microprocessor-based computer control. Apparatus 200 may include a transillumination mask and/or filter which may be microprocessor-based computer controlled. A face plate height adjuster 230 such as for epi-illumination may be included in the apparatus 200.

[0038] Now referring to FIG. 2A, contiguous fiber-optic bundles with receiver channels 306 and illuminating channels 308 with N-to-1 geometry are employed in an alternative nevoscope apparatus. Instead of an incandescent illuminator with a set of filters to produce visible and near-infrared wavelengths, multiple multispectral surface mount light emitting diodes (LED) (such as are commercially available from Lumex Inc.) with an appropriate multiplexed LED driver (such as LTC3219 from Linear Technology) to drive the LEDs can be used. A square pulse signal may be used to turn on the LEDs. A control signal can be used to direct the turn-on time (preferably in the order of a few milliseconds) to a particular LED at a time in a sequential manner. The control and the turn-on signal pulses can be synchronized to the camera frame rate in order to ensure each picture frame corresponds to a single LED illumination.

[0039] Fiber optic bundles (such as those made from high quality silica with better than 99% transmission from 450 nm to 960 nm, available through SCHOTT North America) can be used also in a contiguous manner rather than distributed manner for illumination and receiver channels in a particular geometry (such as in alternate mode or N-to-1 mode in which there are N number of receiver channels for each illumination channel; where N is positive integer preferably N=1, 2, . . . 64; when N=1, it becomes the alternate mode).

[0040] Fiber optic bundles can be directly divided into illumination and receiver multi-core channels where multi-core illumination channels are connected with multispectral LEDs through multiplexed LED drivers (preferably LTC3219, Linear Technology) and receiver channels are connected to a CCD camera (such as Sony ICX415 CCD ) through a focusing optical lens.

[0041] Now referring to FIG. 3, face plate 220 includes outer ring illuminators 222 and 224 and imaging area 226 with distributed sources in alternate positions in a matrix corresponding to the source mask. Outer ring illuminator 222 includes shorter wavelength transillumination fiber channels 223 oriented at 45 degree convergent beam for multispectral or specific wavelength-based background image excitation. Outer ring illuminator 224 includes longer (relative to illuminator 222) wavelength transillumination fiber channels 225 oriented at 45 degree convergent beam for multispectral or specific wavelength-based background image excitation. Imaging area 226 includes plural illumination fiber channels 227 operably connected to an illuminator-filter assembly and plural receiver fiber channels 228 operably connected to the CCD image sensor 270.

[0042] Transillumination imaging is achieved through 45 degree convergent beams through fibers distributed along outer ring(s) 222 and 224 as described above with separate illumination, transilumination masks 280 and multispectral filter 290 selection and control. The described transillumination using any selected optical wavelength can be used for tissue excitation such as for fluoroscopy and/or simple transillumination imaging of background skin/medium.

[0043] Face plate 220 may be removably attached to housing 210 and is operably connected to a removable lens/polarizer 260 that provides an interface to the imaging areas. The fibers 227 and 228 are distributed over the imaging area in a matrix that can be radially symmetric or rectangular (as shown). The fibers 227 and 228 from the imaging area 226 split into fiber cables 250 and 252 with the same access to imaging area 226.

[0044] Now referring to FIG. 2, for imaging a specific imaging mask may be created and used in the illumination path of fiber cable 252. The mask has openings for the desired fibers to be used as source locations to send light into the specimen/skin-tissue 100. As will be apparent to the skilled artisan, many possible schemes for illumination may be employed by fiber placement.

[0045] For multi-spectral imaging, any optical filter of a specific wavelength pass or band filter can be used in the illumination pathway of fiber cable 252 and selected through a computer controlled interface.

[0046] For recording images, an appropriate mask is used to receive the light from the fibers that are not used for illumination. The received light is passed through a focusing lens/polarizer 260 to a CCD sensor 270 to form images and record measurements. The corresponding receiver mask can be appropriately selected through a computer controlled interface.

[0047] Now referring to FIG. 4A, in an epi-illumination mode, light beams L are reflected from above the skin surface 100 and diffused reflected light RL is collected by CCD sensor 270 (not shown) through the optical assembly of cross-polarizer and magnifying lens. This image carries apparent characteristic of a lesion that can be used for automatic diagnostic algorithm in terms of the “ABCD” rule. Now referring to FIG. 4B, in a trans-illumination mode, photons of white-light spectrum or a specific wavelength are directed by a transilluminator ring light such as described in FIG. 3 providing an optical interface for light L to enter into the surrounding area of a skin-lesion with a cone-beam making a forty-five degree angle with respect to the normal of skin 100. The back-scattered diffused light DL re-emerges from the skin and captured by the CCD sensor 270 (not shown) such as a camera through the optical assembly. This image contains the information about absorption and scattering properties of the chromophores of underlying skin layers and lesion. Within the optical tomography framework, it is possible to retrieve the distribution of melanin and blood as two key signature variable for early detection of malignant melanoma. In one embodiment, a shape-based multi-constrained algorithm is applied for reconstruction results, which algorithm overcomes the intrinsic problems associated with reconstructing accurate tomographic information.

[0048] Now referring to FIG. 5 a method of reconstruction is described in terms of Nevoscope transillumination images. The method minimizes the difference between the actual (ΔM^real) and predicted (ΔM^cal); where ΔM^real is the actual measurement vector obtained from multispectral images and ΔM^cal is the computed measurement vector obtained from the reconstructed images measurement by employing genetic algorithms providing global searching. A linearized forward model is adopted and evaluated by Monte Carlo simulation in terms of typical optical properties of normal skin. In one aspect, malignant melanoma is represented by shapes of its melanin part and blood part. These parameters are grouped into genetic algorithms.

[0049] To develop a reconstruction strategy, a forward model is required to relate the measurement to the optical properties of tissue under investigation. Regardless of what kind of imaging geometry is used, an optical system may be described as

M=F(x)   (1)

where M is the measurement and F is a forward model. x is a distribution of unknown optical properties. Given a reasonable initial guess x₀ of the background optical properties, Equation (1) may be expanded into

$\begin{array}{cc}M=F\left(x_0\right)+F^{\prime }\left(x_0\right)\left(x-x_0\right)+\frac{1}{2}F^{″}\left(x_0\right)\left(x-x_0\right)+\dots & \left(2\right)\end{array}$

where F′ and F″ are first order and second order Frechet derivatives respectively.

[0050] Let ΔM=M−F(x₀) and Δx=x−x₀, Equation (2) may be re-arranged as

$\begin{array}{cc}\Delta M=F^{\prime }\Delta x+\frac{1}{2}F^{″}\Delta x+\dots & \left(3\right)\end{array}$

[0051] The discrete form of Equation (3) turns out to be

$\begin{array}{cc}\Delta \stackrel{r}{M}=J\Delta \stackrel{r}{x}+\frac{1}{2}H\Delta \stackrel{r}{x}+\dots & \left(4\right)\end{array}$

Here, J is the Jacobian matrix and H is the Hessian Matrix.

$\Delta \stackrel{1}{M}$

is the measurement vector and

$\Delta \stackrel{1}{x}$

is the vector that gives the variations from the background

${\stackrel{v}{x}}_0.$

Neglecting higher order terms in Equation (4), a simplified linear system can be expressed as

$\begin{array}{cc}\Delta \stackrel{1}{M}=J\Delta \stackrel{r}{x}& \left(5\right)\end{array}$

[0052] The formulation in terms of Equation (5) leads to linear optical tomography which is also known as “difference imaging” with two measurements. One is for background tissue (that is, x₀) and another is for abnormal tissue (that is, unknown x). The difference is then fed to the reconstruction algorithm to obtain the optical properties. In this method, such a linear approach is adopted for nevoscope and the Jacobian matrix is extracted by Monte Carlo simulation.

[0053] Monte Carlo (MC) simulation is viewed as the “gold standard” to modeling light transportation and has been extensively used by researchers in optical communities. One of the major types of MC simulations predicts the trajectories of photons in an absorption-free medium then attenuates the weights of photons in terms of microscopic Beer's law. Using this approach, the trajectories of photons can be stored and used repeatedly which speeds up the calculation and absorption properties in the model and can be easily updated. A Monte Carlo procedure based on this approach and tuned to Nevoscope geometry is described as follows.

[0054] A nevoscope uses a ring light with radius r and incident angle 45 degree with respect to the normal of skin surface. In the coordinate system used in the experiments that follow, skin surface is represented by x-y plane at z=0 and z axis points inward into the medium. The light source is realized by launching a photon at a random position (x₀, y₀, 0) on the ring, where

x₀=r cos(2πξ) and y₀=r sin(2πξ)   (6)

where, ξ is a uniform random number from 0 to 1. Because the discrepancy of refractive indices between air and skin, the photon will change its incident direction which is governed by Snell's law

n_i sin(α_i)=n_t sin(α_t)   (7)

[0055] As the photon loses part of its energy due to specular reflection, considering each photon has normalized energy 1, the remaining energy is given by

$\begin{array}{cc}{w}_0=1-\frac{1}{2}\left[\frac{{\sin }^2\left({\alpha }_i-{\alpha }_t\right)}{{\sin }^2\left({\alpha }_i+{\alpha }_t\right)}+\frac{{\tan }^2\left({\alpha }_i-{\alpha }_t\right)}{{\tan }^2\left({\alpha }_i+{\alpha }_t\right)}\right]& \left(8\right)\end{array}$

[0056] For Equations (7) and (8), α_i=45° is the incident angle and α_t is the transmitted angle. n_i=1 is the refractive index of the ambient and n_t is that of stratum corneum.

[0057] After the photon is launched, it undergoes multiple scattering before termination. The probability of path length that a photon would encounter a scattering event is given by Sobol, The Monte Carlo method (Chicago, Ill.: University of Chicago Press), 1974:

p(l)=μ_s exp(−μ_sl)exp(−μ_αl)+μ_αexp(−μ_αl)exp(−μ_sl)   (9)

[0058] To predict the photon trajectory first, absorption may be neglected temporarily. That is, the simulation is done with an absorption free skin model. Thus Equation (9) reduces to

p(l)=μ_s exp(−μ_sl)   (10)

So the path length between scattering events is derived as

$\begin{array}{cc}l=-\frac{\ln \left(\xi \right)}{{\mu }_s}& \left(11\right)\end{array}$

[0059] When a photon completes its one step, it will scatter into a new direction. The direction is characterized by deflection and azimuthal angles. The azimuthal angle is sampled from a uniform distribution from 0 to 2π but the deflection angle has to be sampled from a phase function which represents the typical forward scattering of biological tissue. A phase function may be expressed as

$\begin{array}{cc}p\left(\cos \theta \right)=\frac{1}{4\pi }\frac{1-g^2}{{\left(1+g^2-2g\cos \theta \right)}^{\frac{3}{2}}}& \left(12\right)\end{array}$

[0060] In addition, it is quite likely that a photon will travel across internal refractive index mismatched interfaces or hit the skin surface during its transportation. If a photon hits the surface, its energy can be separated into two parts: one part escapes the surface and the other part is internally reflected to propagate with reduced energy. The ratio between internal reflected energy and the original energy is given by Fresnel reflection coefficient.

$\begin{array}{cc}R\left({\alpha }_i\right)=\{\begin{array}{cc}\frac{{\left(n_t-n_i\right)}^2}{{\left(n_t+n_i\right)}^2}& \left(\mathrm{if}{\alpha }_i=0\right)\\ \frac{1}{2}\left[\frac{{\sin }^2\left({\alpha }_i-{\alpha }_t\right)}{{\sin }^2\left({\alpha }_i+{\alpha }_t\right)}+\frac{{\tan }^2\left({\alpha }_i-{\alpha }_t\right)}{{\tan }^2\left({\alpha }_i+{\alpha }_t\right)}\right]& \left(\mathrm{if}0<{\alpha }_i<{\sin }^{-1}\left(\frac{n_i}{n_t}\right)\right)\\ 1& \left(\mathrm{if}{\sin }^{-1}\left(\frac{n_i}{n_t}\right)<{\alpha }_i<\frac{\pi }{2}\right)\end{array}& \left(13\right)\end{array}$

[0061] In Equation (13), α_i and α_t are angle of incidence and transmittance related by Snell's law. n_i and n_t are medium of incidence and transmittance respectively.

[0062] In the case of a photon crossing internal interfaces, the “Russian roulette” approach proposed by Wang et al., MCML-Monte Carlo modeling of light transport in multi-layered tissues, Computer Methods and Programs in Biomedicine 47, 131-146 (1995), may be employed. Once the photon hits an interface, the Fresnel reflection coefficient R(α_i) is calculated and compared to a uniform-distributed random number ξ. If R(α_i)>ξ, the photon will cross the boundary. Otherwise, it is totally reflected.

[0063] Following the rules described above, a photon is launched from the ring source and propagates through the optical skin model until terminated or detected. If its total path length exceeds some preset value, the photon gets terminated. If it is captured by a detector, its energy is re-weighted in terms of its trajectory as

$\begin{array}{cc}{w}_d=\left(1-R\left({\alpha }_h\right)\right)w_0\sum _{i=1}^{h-1}R\left({\alpha }_i\right),& \left(14\right)\end{array}$

where it is assumed the photon hits the skin surface h times during its transportation.

[0064] After all trajectories of photons are recorded for a specific detector, the total received intensity is therefore given in terms of microscopic Beer's law

$\begin{array}{cc}M=\sum _{i=1}^pw_d^i\sum _{j=1}^{q_p}{}^{-{\mu }_a\left(j\right)l_j^i}& \left(15\right)\end{array}$

where w_dⁱ of ith photon is given by equation (14), μ_α(j) is the absorption coefficient along the path l_jⁱ of ith photon. It is assumed that the detector received p photons while each of them used q_p steps. The wavelength dependent absorption coefficients are obtained from published data. See, R. Anderson et al., The optics of human skin, Journal of Investigative Dermatology, Vol. 77(1), 13-19 (1981); A. Bashkatov et al., Optical properties of human skin, subcutaneous and mucous tissues in the wavelength range from 400 to 2000 nm, Journal of Physics D: Applied Physics, Vol. 38, 2543-2555 (2005); M. Van Gemert et al., Skin Optics, IEEE Transactions on Biomedical Engineering, Vol. 36(12), 1146-1155 (1989). For purposes of the present experiments a seven-layered skin model is employed which includes stratum corneum (20 μm), epidermis (80 μm), dermis (150 μm), upper blood net dermis (80 μm), reticular dermis (1500 μm), deep blood net dermis (100 μm) and subcutaneous fat (6000 μm). See, I. Meglinski et al., Quantitative assessment of skin layers absorption and skin reflectance spectra simulation in the visible and near-infrared spectral regions, Physiol. Meas, 23, 741-753 (2002).

[0065] The absorption coefficient of each layer is a combined effect contributed by several absorbers such as blood, melanin, water and the baseline skin which is free of these absorbers. The absorption coefficient of baseline skin is given by

μ_α^baseline(λ)=7.84×10⁸×λ^−3.255 (cm⁻¹)   (A.1)

[0066] Here, λ is the wavelength measured in nanometers. The absorption of blood is determined by its oxygenation and fraction of red cells in blood. For example, we assume that haematocrit is Ht % and S % hemoglobin is oxy-hemoglobin. Then the absorption of blood is given as

μ_α^blood(λ)=Ht×(1−S)×μ_α^Hb(λ)+Ht×S×μ_α^HbO2(λ)   (A.2)

[0067] Here μ_α^Hb(λ) and μ_α^HbO2(λ) are absorption coefficients of oxy-hemoglobin and deoxy-hemoglobin respectively. The absorption spectra come from compiled data. See, Prahl, Oregon Medical Laser Center, Optical absorption of hemoglobin, http://omlc.ogi.edu/spectra/hemoglobin/index.html (2006).

[0068] The haematocrit has a typical value 45%. The absorption of melanin is given by Meglinski et al.

μ_α^melanin(λ)=5×10¹⁰×λ^3.33 (cm⁻¹)   (A.3)

[0069] In healthy skin, melanin only exists in epidermis layer. Its fraction ranges from 1% to 3% for light skinned Caucasians, from 11% to 16% for well-tanned Caucasians and Mediterraneans, and from 18% to 43% for darkly pigmented Africans. The presence of melanin in deeper skin layers like dermis usually indicates the malignancy of melanoma. The absorption of water is known. Absorption coefficients of two blood free layers stratum comeum and epidermis are predicted by

μ_α^stratum(λ)=((0.1−0.3×10⁻⁴λ)+0.125μ_α^baseline(λ))(1−C_H2O)+C_H2Oμ_α^H2O(λ)   (A.4)

μ_α^epidermis(λ)=(C_melaninμ_α^melanin(λ)+C_H2Oμ_α^H2O(λ)+(1−C_melanin−C_H2O)μ_α^baseline(λ))   (A.5)

[0070] In the above equations, C_melanin is the fraction of melanin in epidermis and C_H2O is the fraction of water in a specific layer. Absorption coefficients of four dermis layers and subcutaneous fat are described by one formula (A.6) while their values depend on how much blood and water they have.

μ_α^bloodlayer(λ)=C_bloodμ_α^blood(λ)+C_H2Oμ_α^H2O+(1−C_blood−C_H2O)μ_α^baseline   (A.6)

[0071] Additional optical properties of skin include scattering coefficients, anisotropy factor and refractive index. Representative values of these optical properties and fractions of absorbers in each skin layer are reported by Meglinski, supra.

[0072] Now referring to FIG. 6A, in accordance with the experiments, the detector space of the nevoscope is divided into k number of rings 202 (for example, three rings 202 as shown in the FIG. 6A) with equal width. The innermost circle 204 contains four elements or detectors 206 with equal size and other rings 202 are split into a number of elements which maintain the same area as the ones in the innermost circle 204. Using this discretization strategy makes use of the symmetrical geometry of a nevoscope. When evaluating Jacobian matrix, instead of obtaining one MC simulation for each detector, since MC is extremely time-consuming, and noting that all detectors 206 on the same ring actually have similar trajectories, the trajectories for one detector 206 on that ring 202 is recorded and trajectories of other detectors 206 are generated by rotating recorded trajectories. So eventually, only k independent simulations are required to evaluate Jacobian matrix.

[0073] Resolution of the reconstructed image is dependent on the size of the Jacobian matrix. Now referring to FIG. 6B, in the experiments, the interrogated volume is divided into N×N×M voxels 300. Each voxel 400 has the size Δc×←y×Δz and has homogeneous absorption coefficient. Where the total number of detectors is T, the measurements of background normal skin is then given by

$\begin{array}{cc}M\left(s\right)=\sum _{i-1}^{p_s}w_{\mathrm{skin}}^i=\sum _{i=1}^{p_s}\left(w_d^i\sum _{j=1}^{N\times N\times M}{}^{-{\mu }_a\left(j\right)l_j^i}\right)s=1,\dots T& \left(16\right)\end{array}$

where p_s is total number of photons received by the detectors. l_jⁱ is the path length of photon i in the voxel j. μ_α(j) is the absorption coefficient of voxel j and is extracted from optical properties of normal skin.

[0074] For difference imaging, the Jacobian matrix is given by

$\begin{array}{cc}{J=\left[\begin{array}{ccccc}\frac{\partial \Delta M\left(1\right)}{\partial {\mathrm{\Delta \mu }}_a\left(1\right)}& \frac{\partial \Delta M\left(1\right)}{\partial {\mathrm{\Delta \mu }}_a\left(2\right)}& & & \frac{\partial \Delta M\left(1\right)}{\partial {\mathrm{\Delta \mu }}_a\left(N\times N\times M\right)}\\ g& g& & & g\\ g& & g& & g\\ g& & & g& g\\ \frac{\partial \Delta M\left(T\right)}{\partial {\mathrm{\Delta \mu }}_a\left(1\right)}& \frac{\partial \Delta M\left(T\right)}{\partial {\mathrm{\Delta \mu }}_a\left(2\right)}& & & \frac{\partial \Delta M\left(T\right)}{\partial {\mathrm{\Delta \mu }}_a\left(N\times N\times M\right)}\end{array}\right]}_{\Delta {\mu }_a\left(j\right)=0}& \left(17\right)\end{array}$

[0075] Each term of Equation (17) is evaluated by recorded photon history. When a tumor grows on the normal skin, variation in measurements ΔM(s) may be expressed as

$\begin{array}{cc}\Delta M\left(s\right)=\sum _{i=1}^{p_s}\left(w_{\mathrm{skin}}^i\sum _{j=1}^{N\times N\times M}{}^{-\Delta {\mu }_a\left(j\right)l_j^i}\right)& \left(18\right)\end{array}$

[0076] Taking the derivative with respect to Δμ_α(j) and set Δμ_α(j) to zero results in

$\begin{array}{cc}{\frac{\partial M\left(s\right)}{\partial \Delta {\mu }_a\left(j\right)}}_{\Delta {\mu }_a\left(j\right)=0}=\sum _{i=1}^{p_s}\left(w_{\mathrm{skin}}^i·l_j^i\right)& \left(19\right)\end{array}$

[0077] The Jacobian matrix used in the experiments is evaluated on the normal skin model introduced hereinbelow with k=20, N=32 and M=20 which divides the nevoscope detector plane of 1.2 cm diameter into 1588 elements and divides the 1.2 cm×1.2 cm×2000 μm field of view into 32×32×20 voxels.

[0078] Inverse Problem

[0079] The Equation (5) can be solved by Singular Value Decomposition (SVD), Algebraic Reconstruction Technique (ART), Conjugate Gradient method (CG) and their variants. However, due to the ill-posed and under-determined nature of the inverse problem, it is unlikely to reconstruct quantitatively accurate and well-localized absorption coefficients. Therefore, shape-based optical tomography is preferred, wherein the abnormalities are assumed to have piece-wise constant optical properties in the compact supported regions and the boundaries of these regions can be effectively approximated by some shape functions with the limited numbers of control parameters. As to the reconstruction problem, the unknowns just include a few numbers of control parameters and, sometimes, the piece-wise constant optical properties.

[0080] Shape Representation of Malignant Melanoma

[0081] A melanoma is essentially a three-dimensional object with distributed pigmentation and chromophores. Malignant melanoma is a result of uncontrolled replication of melasome cells sitting on the basal layer of epidermis. The shape of the melanoma is hence bounded by the epidermis layer. What is needed to describe the lesion is therefore reduced to 2D surfaces in the 3D domain.

[0082] Breaking the malignant melanoma into parts, such as the melanin and the blood parts, permits representation of the malignant melanoma with 2D surfaces. The melanin part is a 3D region bounded by a single surface and the epidermis layer. Within the region, the optical properties are constant and the only absorber is melanin. A second surface which sits below the first surface is used to represent the blood part. The region bounded by the first surface and the second surface is blood only. This model mimics a skin-lesion and its x-z intersection is shown in FIG. 7.

[0083] The first surface is represented as ƒ₁(x, y) which corresponds to the depth of lesion from the epidermis layer at the position (x,y). The continuous surface is represented with limited parameters. First, a N×N rectangular grid is positioned to lie over the epidermal layer. Second, the function ƒ₁(x, y) is sampled to N×N discrete values ƒ_d1(X, Y). Here, (x, y) is continuous and (X, Y) is N×N numbers of discrete sampling positions. Third, the discrete values are interpolated by the cubic tensor-product B-spline which satisfies the following condition:

$\begin{array}{cc}{f}_{d1}\left(X,Y\right)=\sum _{i=1}^N\sum _{j=1}^Nc_1\left(i,j\right){\beta }^3\left(X-i,Y-j\right)& \left(20\right)\end{array}$

and the original function ƒ₁(x, y) can then be approximated by

$\begin{array}{cc}{f}_{B1}\left(x,y\right)=\sum _{i=1}^N\sum _{j=1}^Nc_1\left(i,j\right){\beta }^3\left(x-i,y-j\right)& \left(21\right)\end{array}$

where

β³(x−i,y−j)=β³(x−i)gβ³(y−j)   (22)

is the tensor product of one-dimensional cubic B-spline basis β³(x−i) and β³(y−j), and c₁(i, j) is B-spline coefficient.

[0084] The basis function of B-spline comes from repeated convolution of the box function. It has compact support so that the computation burden is alleviated. The cubic B-spline is continuous up to second derivative which leads to the smoothness of the interpolated surface and in essence puts regularization effect on the reconstructed surface in our shape-based approach. In the end, for the melanin part, the parameters to determine its shape are reduced to N×N discrete values ƒ_d1(X, Y) or, equivalently, the corresponding B-spline coefficients c₁(i, j).

[0085] Similarly, the second surface can be defined by N×N discrete values ƒ_d2(X, Y) or, equivalently, z_d(X, Y)=ƒ_d2(X, Y)−ƒ_d1(X, Y) which is the thickness of blood region between the first surface and the second surface.

[0086] Shape-Based Reconstruction

[0087] To reconstruct the surfaces and piecewise constant optical properties, the continuous surface representation should be incorporated into the forward photon transportation model. In our study, the linearized forward model is kept intact but the continuous representation is sampled into the discrete vector of unknowns

$\Delta \stackrel{1}{x}.$

[0088] $\begin{array}{cc}\Delta \stackrel{1}{M}=\mathrm{Jg}\Delta \stackrel{r}{x}=\mathrm{JgS}\left(f_{d1}\left(X,Y\right),\mathrm{mp}1,z_d\left(X,Y\right),\mathrm{bp}2\right)& \left(23\right)\end{array}$

[0089] where S(g) is the sampling function which converts the continuous shape representation to the voxel-based optical properties in forward model, mp1 is the fraction of melanin and bp2 is fraction of blood. The inverse problem can therefore be formulated as minimizing the objective function

$\begin{array}{cc}{F}_{\mathrm{obj}}=\frac{1}{2}{\Delta \stackrel{r}{M}-\mathrm{Jg}\Delta \stackrel{r}{x}}^2=\frac{1}{2}\Delta \stackrel{r}{M}-{\mathrm{JgS}(f_{d1}\left(X,Y\right),\mathrm{mp}1,z_d\left(X,Y\right),\mathrm{bp}2}^2& \left(24\right)\end{array}$

[0090] The unknowns of this inverse problem are reduced to ƒ_d1(x, y), z_d(X, Y), mp1 and bp2. The multi-spectral shape reconstruction using N wavelengths can be formulated as a multi-objective optimization problem and its objective function is given as

F_obj=α₁gF_obj^λ1+α₂gF_obj^λ2+ . . . +α_NgF_obj^λN   (25)

[0091] where, {α₁, α₂, . . . α_N} is a set of coefficients to balance the contributions from different single-wavelength objective functions.

[0092] The wavelength dependent objective functions are given as

$\begin{array}{cc}\begin{array}{c}\begin{array}{c}\begin{array}{c}{F}_{\mathrm{obj}}^{{\lambda }_1}=\frac{1}{2}{\Delta \stackrel{r}{M^{{\lambda }_1}}-J^{{\lambda }_1}g\Delta {\stackrel{r}{x^{{\lambda }_1}}}^{}}^2=\frac{1}{2}{\Delta \stackrel{r}{M^{{}_{{\lambda }_1}}}-J^{{\lambda }_1}{\mathrm{gS}}^{{\lambda }_1}\left(f_{d1}\left(X,Y\right),\mathrm{mp}1,z_d\left(X,Y\right),\mathrm{bp}2\right)}^2\\ F_{\mathrm{obj}}^{{\lambda }_2}=\frac{1}{2}{\Delta \stackrel{r}{M^{{\lambda }_2}}-J^{{\lambda }_2}g\Delta \stackrel{r}{x^{{\lambda }_2}}}^2=\frac{1}{2}{\Delta \stackrel{r}{M^{{}_{{\lambda }_2}}}-J^{{\lambda }_2}{\mathrm{gS}}^{{\lambda }_2}\left(f_{d1}\left(X,Y\right),\mathrm{mp}1,z_d\left(X,Y\right),\mathrm{bp}2\right)}^2\end{array}\\ g\end{array}\\ g\\ g\\ F_{\mathrm{obj}}^{{\lambda }_N}=\frac{1}{2}{\Delta \stackrel{r}{M^{{\lambda }_N}}-J^{{\lambda }_N}g\Delta \stackrel{r}{x^{{\lambda }_N}}}^2=\frac{1}{2}{\Delta \stackrel{r}{M^{{\lambda }_N}}-J^{{\lambda }_N}{\mathrm{gS}}^{{\lambda }_N}\left(f_{d1}\left(X,Y\right),\mathrm{mp}1,z_d\left(X,Y\right),\mathrm{bp}2\right)}^2\end{array}& \left(26\right)\end{array}$

[0093] Various optimization techniques have been reported to solve the shape-based reconstruction problems but all have drawbacks. Hence, a genetic algorithm is used because the gradient need not be evaluated which simplifies the computation and provides reliability. Genetic algorithm is one of the most popular methods used in seeking global minimal. Among the global optimization techniques, genetic algorithm provides a reasonable convergence rate due to its implicit parallel computation. See, Whitely, A Genetic Algorithm Tutorial, Technical Report CS-93-103, Colorado State University (1993); Busetti, Genetic algorithms overview, http://www.geocities.com/francorbusetti/gaweb.pdf; Beasley et al., An overview of genetic algorithms: Part 1, Fundamentals, University Computing, 15(2), 58-69 (1993); Beasley et al., An Overview of Genetic Algorithms: Part 2, Research topics, University Computing, University of Cardiff, 15(4), 170-181 (1993); Ingber et al., Genetic Algorithms and Very Fast Simulated Annealing: Comparison, Mathematical and Computer Modelling 16(11), 87-100 (1992); Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learning, Addison-Wesley Pub. Co. (1989); Houck et al., A Genetic Algorithm for Function Optimization: a Matlab Implementation, NCSU-IE TR 95-09 (1995).

[0094] The following is an implementation of the genetic algorithm as it applies to the present invention. As an example, an optimization problem may seek the maximal value of the objective function

Val=ƒ(α₁, α₂, α₃)   (27)

[0095] This objective function itself is an appropriate fitness function because its value measures how good the solution is. The parameters in the objective function are consequently encoded to the chromosome as genes. It is believed a real number representation is more efficient compared to the binary version. The chromosome with a number of real-numbered genes is given as

{α₁, α₂, α₃}  (28)

[0096] The first step in genetic algorithm is to give an initial population with a number of individuals (that is, with different chromosomes) to reflect the variety. This population is usually generated randomly. This step in fact means to randomly sample the parameter space. Different from random searching technique, the genetic algorithm realizes implicit parallel computation so that it has a faster convergence rate. In terms of schema theory, the fitness of an individual not only represents its optimality at that specific position but gives fitness of some underlying pattern known as schema. After the population has been initialized, the fitness score of each individual is evaluated to facilitate the reproduction step.

[0097] Reproduction reflects the natural reality that better fitted individuals have higher probabilities to appear in the next generation. Given the fitness scores of individuals, the next generation is selected in terms of several selection rules. A common method is to select the next generation through the “Roulette Wheel” rule. If the individual i has the fitness score Val_i, then its probability to be selected into next generation is given as

$\begin{array}{cc}P=\frac{{\mathrm{Val}}_i}{\stackrel{\_}{V}\mathrm{al}}& \left(29\right)\end{array}$

where Val is the average fitness of the population with M individuals

$\begin{array}{cc}\stackrel{\_}{V}\mathrm{al}=\frac{\sum _{i=1}^M{\mathrm{Val}}_i}{M}& \left(30\right)\end{array}$

[0098] The selection process can be viewed as putting the individuals on a wheel and each individual occupies a pie representing its probability. Randomly selecting M individuals from the wheel results in the next generation.

[0099] Crossover is the process that one offspring is generated by two parents. The offspring could become a worse-fitted individual. However, there is a good chance that the offspring has better fitness than both of its parents. Crossover is believed to be the major operator that makes the genetic algorithm have the capability to investigate the parameter space efficiently. One special and efficient real crossover operator is the mathematical crossover and is defined as

$\begin{array}{cc}\left\{{\alpha }_1^1,{\alpha }_2^1,{\alpha }_3^1\right\}\Rightarrow \left\{p·{\alpha }_1^1+\left(1-p\right)·{\alpha }_1^2,p·{\alpha }_2^1+\left(1-p\right)·{\alpha }_2^2,p·{\alpha }_3^1+\left(1-p\right)·{\alpha }_3^2\right\}\left\{{\alpha }_1^2,{\alpha }_2^2,{\alpha }_3^2\right\}\Rightarrow \left\{\left(1-p\right)·{\alpha }_1^1+p·{\alpha }_1^2,\left(1-p\right)·{\alpha }_2^1+p·{\alpha }_2^2,\left(1-p\right)·{\alpha }_3^1+p·{\alpha }_3^2\right\}& \left(31\right)\end{array}$

where p is a uniform random number from zero to one.

[0100] Mutation is another operator to change the chromosome of an individual. It gives each gene in the chromosome a small probability to be set to a random number. This operator gives any region in the parameter space some probability to be investigated. It therefore helps to escape the local minima. However, over-emphasis on the mutation operator leads to a stochastic random search.

[0101] Starting from the initial population and going over one generation after another, finally the genetic algorithm converges to the optimal solution. As to the optimization problem occurred in shape-based reconstruction, the objective function is selected as the fitness function. Its parameters is coded into chromosome as

$\begin{array}{cc}\left\{\begin{array}{c}{f}_{d1}\left(X_1,Y_1\right)-f_{d1}\left(X_2,Y_2\right)-\dots -f_{d1}\left(X_{N\times N},Y_{N\times N}\right)-\mathrm{mp}1-\\ z\left(X_1,Y_1\right)-z\left(X_2,Y_2\right)-\dots -z\left(X_{N\times N},Y_{N\times N}\right)-\mathrm{bp}2\end{array}\right\}& \left(32\right)\end{array}$

[0102] Reasonable constraints are added to the parameters in preparation of the optimization algorithm. First, the region of support of the lesion in x-y plane is readily available in terms of the epi-illumination nevoscope image (FIG. 4A). This further reduces the number of parameters to represent a surface from N×N to a smaller set. As a consequence, the optimization algorithm has a faster convergence rate. Second, the blood region typically comprises thin layers within a few hundreds of micrometers which put a constraint on z(X, Y). Third, the fractions of melanin and blood are not free parameters. They can also be bounded according to the appearance of melanoma and the clinical experience. Last, multi-spectral imaging provides implicit constraints. Due to the distinct absorption spectra of blood and melanin, a reasonable solution must also satisfy the measurements of all involved wavelengths.

[0103] In summary, the method of solving the shape-based reconstruction problem using genetic algorithm to involves the following steps:

[0104] Step 1: Initialize the population. The values of genes are randomly generated as well as subject to some specific constraints. The size of the population is sufficient large to reflect randomness.

[0105] Step 2: Evaluate the fitness. Each chromosome in the population is evaluated in terms of the objective function and, as a result, associated with a fitness score. This step implicitly includes sampling the shape-representation to fill into the forward model.

[0106] Step 3: Evolve the population. This step includes reproduction in terms of the fitness score, arithmetic crossover and mutation.

[0107] Step 4: Check the termination condition. The termination condition in the experiment is the number of iterations of the algorithm. Once it exceeds the pre-defined number, the algorithm stops. Otherwise, it returns to step 2 and continues to evolve.

[0108] To validate the shape-based multi-spectral algorithm, a double-surface model was created to represent malignant melanoma. The first and second surfaces are described by a mixed Gaussian function which are given as

$\begin{array}{cc}{f}_1\left(x,y\right)=\mathrm{MAX}\left(\mathrm{peak}1\mathrm{gG}\left(x,y,{\mu }_{1a},{\mu }_{2a},{\sigma }_a\right),\mathrm{peak}2\mathrm{gG}\left(x,y,{\mu }_{1b},{\mu }_{2b},{\sigma }_b\right)\right)f_2\left(x,y\right)=\mathrm{MAX}\left(\mathrm{peak}3\mathrm{gG}\left(x,y,{\mu }_{1a},{\mu }_{2a},{\sigma }_a\right),\mathrm{peak}4\mathrm{gG}\left(x,y,{\mu }_{1b},{\mu }_{2b},{\sigma }_b\right)\right)& \left(33\right)\end{array}$

where the Gaussian function is

$\begin{array}{cc}G\left(x,y,{\mu }_1,{\mu }_2,\sigma \right)=\frac{1}{2{\mathrm{\pi \sigma }}^2}\exp \left(-\frac{{\left(x-{\mu }_1\right)}^2+{\left(y-{\mu }_2\right)}^2}{2{\sigma }^2}\right)& \left(34\right)\end{array}$

The parameters used in Equation (33) are

$\begin{array}{cc}\{\begin{array}{cccc}{\mu }_{1a}=0.0375\times 2& \mathrm{cm}& {\mu }_{2a}=-0.0375\times 3& \mathrm{cm}\\ {\mu }_{1b}=-0.0375\times 2& \mathrm{cm}& {\mu }_{2b}=0.0375\times 3& \mathrm{cm}\\ {\sigma }_a=0.0375\times 2.5& \mathrm{cm}& {\sigma }_b=0.0375\times 2& \mathrm{cm}\\ \mathrm{peak}1=100\times 6& \mathrm{\mu m}& \mathrm{peak}2=100\times 4& \mathrm{\mu m}\\ \mathrm{peak}3=100\times 8& \mathrm{\mu m}& \mathrm{peak}4=100\times 6& \mathrm{\mu m}\end{array}& \left(35\right)\end{array}$

[0109] The fraction of melanin is set to 5% between the epidermal layer and the first surface ƒ₁(x, y). The fraction of blood is set to 20% between the first surface ƒ₁(x, y) and the second surface ƒ₂(x, y). This model has sufficient variation in order to verify the reliability of the reconstruction algorithm. FIG. 8A displays the 3D view of this model.

[0110] The double-surface continuous model is sampled and two measurements ΔM⁵⁸⁰ and ΔM⁸⁰⁰ are calculated by Monte Carlo simulation at 580 nm and 800 nm respectively. Next a 9×9 rectangular grid is overlapped on the epidermal layer. The region of support of the lesion is counted as 10 discrete control points. As a result, the chromosome contains 22 genes and is coded as

$\begin{array}{cc}\left\{\begin{array}{c}{f}_{d1}\left(X_1,Y_1\right)-f_{d1}\left(X_2,Y_2\right)-\dots -f_{d1}\left(X_{10},Y_{10}\right)-\mathrm{mp}1-\\ z\left(X_1,Y_1\right)-z\left(X_2,Y_2\right)-\dots -z\left(X_{10},Y_{10}\right)-\mathrm{bp}2\end{array}\right\}& \left(36\right)\end{array}$

[0111] The fitness function of the genetic algorithm is given as

F_obj=α₁gF_obj⁵⁸⁰+α₂gF_obj⁸⁰⁰   (37)

where α₁ and α₂ can be chosen on the relative weight assigned to melanin contents to the entire volume. About 30% of the entire volume to be pigmented and used α₁=0.3 and α₂=1 was considered in the experiment.

[0112] Four simulations with different constraints were implemented in real number presented genetic algorithm. Results are summarized in Table 1. The recovered surfaces are displayed in FIG. 8B-E. In each case, the left is the first surface and the right is the second surface. There is no constraint for the first surface while the thickness between the first surface and the second surface is set to be 300 μm to represent a thin layer of blood net.

[0113] In an example there are 100 individuals in a population and the genetic algorithm terminates after 150 iterations. Now referring to FIG. 9, the fitness scores are normalized and displayed with respect to the number of iteration. As can be seen the genetic algorithm converges in all cases. Now referring to FIGS. 10A and 10B, the deformation process is shown, suggesting that the genetic algorithm has successfully recovered both surfaces.

[0114] To further evaluate the reconstruction result, volume deviations Volerr1 and Volerr2 were introduced. They are defined as

$\begin{array}{cc}\mathrm{Volerr}1=\frac{\mathrm{vol}1_c-\mathrm{vol}1_m}{\mathrm{vol}1_m}& \left(38\right)\end{array}$

[0115] where vol1_c is the calculated volume bounded by the first surface and vol1_m is the corresponding volume from the model and

$\begin{array}{cc}\mathrm{Volerr}2=\frac{\mathrm{vol}2_c-\mathrm{vol}2_m}{\mathrm{vol}2_m}& \left(39\right)\end{array}$

where, vol2_c is the calculated volume bounded by the first and second surfaces and vol2_m is the corresponding volume from the model.

TABLE 1
Summary of Reconstruction Results
	Melanin	Blood	Recovered	Recovered
	Bounds	Bounds	Melanin	blood	Volerr1	Volerr2
Case	(%)	(%)	(%)	(%)	(%)	(%)
(b)	5-5	20-20	5	20	2.68	16.58
(c)	4.5-5.5	10-30	5.10	18.33	3.81	20.89
(d)	4-6	10-30	4.64	18.97	5.18	24.71
(e)	3-7	10-30	3.37	10.08	44.09	63.50

[0116] As shown in Table 1, the constraints have significant impacts on the reconstructed surfaces. Except the loosest constrained case (e), all cases present reasonable reconstructions which are consistent with the model. Moreover, the reconstructed first surface has a smaller volume error than the second surface. There are several reasons to explain the larger volume error of the second surface. First, the absorption coefficient of blood is smaller than that of melanin. As a result, the change in blood region has less contribution to the fitness function. Second, since a reflectance geometry is adopted in the nevoscope, the sensitivity decreases at deeper layers. This also influences the accurate reconstruction of the second surface. Third, because the two surfaces are attached together, the error resulting from the first surface would inevitably propagate to the second surface.

[0117] The inventive reconstruction techniques are shown to be a good solution to the ill-posed inverse problem. Using additional wavelengths increases the ability to isolate the true solution from a number of possible solutions.

[0118] Monte-Carlo Simulation-Based Mapping of Multispectral Image Information onto a 3-D Volume.

[0119] Now referring to FIG. 11, in another embodiment a Monte-Carlo simulation-based mapping of multispectral image information onto a 3-D volume is provided. The multispectral image information obtained can be projected onto a 3-D discrete reconstruction volume through the knowledge-based model derived from a Monte-Carlo simulation providing probability distribution functions of energy and path length and energy in the 3-D sub-surface space. The Monte-Carlo simulation of the skin-tissue model (without the lesion) provides a probability distribution of each photon path between the source and detector in the simulated medium. An iterative ART algorithm can be used to employ this matrix of initial weights and path-length information to compute the diffused reflectance distribution that is compared with the acquired multispectral images of skin of the human subjects (near the skin-lesion). Thus, the weights and path-length information of the reconstruction space of the skin is updated using an Algebraic Reconstruction Technique (ART) (See, S. Srinath Maganti, “Three-Dimensional Optical Tomographic Image Reconstruction from Nevoscope Images”, Ph.D. Dissertation, Advisor: Atam P. Dhawan, University of Cincinnati (1997); S. Maganti and A. P. Dhawan, “3-D Nevoscope image reconstruction using diverging ray ART”, Proceedings SPIE International Conference on Biomedical Optics (1997)) for minimization of the error between the computed diffused reflectance and acquired multispectral images of skin. This updated weight matrix and path-length information in the 3-D reconstruction space now serves as initial weight matrix for reconstructing the skin-lesion information from multispectral images of the skin-lesion using the iterative ART. Thus, Monte Carlo simulation provides the forward model of imaging for the voxel weights and path-length information which is then used in reconstruction through an iterative Algebraic Reconstruction Technique (ART) (See, S. Srinath Maganti, “Three-Dimensional Optical Tomographic Image Reconstruction from Nevoscope Images”, Ph.D. Dissertation, Advisor: Atam P. Dhawan, University of Cincinnati (1997); S. Maganti and A. P. Dhawan, “3-D Nevoscope image reconstruction using diverging ray ART”, Proceedings SPIE International Conference on Biomedical Optics (1997)) as depicted in FIG. 11.

[0120] It is assumed a skin lesion property (such as melanin, hemoglobin, or de-oxy-hemoglobin concentration) is represented by F_i(x,y,z); i=1, . . . , I where I is the number of skin-lesion properties. It is also assumed that the Monte-Carlo simulation and iterative ART-based reconstruction of skin provides the updated initial weight matrix with 3-D distributions of photon energy density represented by E_j(α,β,γ) with path lengths P_k,i,j((α_l,β_m,γ_n)⁰, E_j;lmn), (α_l,β_m,γ_n)¹, E_j;lmn), . . . , (α_l,β_m,γ_n)^Z, E_j;lmn)); k=1, . . . , K where K is the number of paths between source and detector for a specific imaging wavelength j (j=1, . . . , J) for imaging the skin-lesion property F_i(x,y,z). Each path length is a list of voxel coordinates (α_l,β_m,γ_n)^z; z=1, . . . , Z (where Z is the number of voxels between source and detector) with their respective associated energy level E_j;lmn. The normalized probability distribution of photon energy provides initial weights for reconstructing the object information F_i(x,y,z) as RF_i(x′,y′,z′) where RF_i is the respective reconstructed property in the 3-D reconstruction space (x′,y′,z′). Let NEVO_TRANS function describe the imaging the skin-lesion property, F_i(x,y,z) through Nevoscope Transillumination with the forward model:

[M_j(α,β,γ)]=NEVO_TRANS{F_i(x,y,z)}=[F_i(x,y,z)].[W_j(α,β,γ)]

where [M_j(α,β,γ)] is the measured or acquired image, [W_j(α,β,γ)] is the updated weight matrix, which is approximated and modeled through Monte-Carlo simulation with [E_j(α,β,γ)].

[0121] The objective of the reconstruction process is to obtain an estimate RF_i(x′,y′,z′) of the object property F_i(x,y,z). The ART algorithm starts projecting the acquired image data onto the reconstruction space of voxels using the initial weights and path length information obtained from the Monte-Carlo simulation (forward model). In each iteration, the voxel weight and path length information is used to assign new weight values such that the error between the reconstructed diffused reflectance output from the medium using the new voxel weights and the acquired image is minimized. Thus, the reconstruction process concludes that the incident photon energy radiation on the reconstructed object provides the diffused reflectance output as close as possible to the acquired image as

[RF_j(x′,y′,z′)]=[P_k,i,j((α_l,β_m,γ_n)⁰, E_j;lmn), (α_l,β_m,γ_n)¹, E_j;lmn), . . . , (α_l,β_m,γ_n)^Z, E_j;lmn))]^INV. [M_j(α,β,γ)]=[W_j(α,β,γ)]^INV.[M_j(α,β,γ)]

[0122] Thus, through the ART reconstruction process an inverse of the weight matrix, [W_j(α,β,γ)]^INV is obtained through the iterative process. However, since the matrix is very large to invert, the convergence can be obtained using a hierarchical partitioning. Also, the a priori knowledge of the path length from the Monte-Carlo simulation can be used to reduce the size of the sub-matrices to be inverted. Thus, the assignment of the acquired image data onto a reconstructed 3-D space can be efficiently performed for each property of the skin-lesion through specific wavelength based imaging. Thus, the individual wavelength based reconstruction from each of the acquired multispectral images will be obtained. The specific property—(such as melanin, hemoglobin, etc.) based 3-D reconstructed distributions would be obtained through combining the multispectral reconstructions weighted by their respective volume fraction model described in the equations above. The total melanin fraction would be obtained from multispectral images as

$\left[{\mathrm{RF}}_{\mathrm{melanin}}\left(x^{\prime },y^{\prime },z^{\prime }\right)\right]=\sum _{j=1}^J\left(f_{j,\mathrm{epmel}}+f_{j,\mathrm{pmel}}+f_{j,\mathrm{rmel}}\right)\left[{\mathrm{RF}}_j\left(x^{\prime },y^{\prime },z^{\prime }\right)\right]\left[{\mathrm{RF}}_{\mathrm{oxy}\text{-}\mathrm{hemo}}\left(x^{\prime },y^{\prime },z^{\prime }\right)\right]=\sum _{j=1}^J\left(f_{j,\mathrm{poxy}}+f_{j,\mathrm{roxy}}\right)\left[{\mathrm{RF}}_j\left(x^{\prime },y^{\prime },z^{\prime }\right)\right]\left[{\mathrm{RF}}_{\mathrm{de}\text{-}\mathrm{oxy}\text{-}\mathrm{hemo}}\left(x^{\prime },y^{\prime },z^{\prime }\right)\right]=\sum _{j=1}^J\left(f_{j,\mathrm{pdeoxy}}+f_{j,\mathrm{rdeoxy}}\right)\left[{\mathrm{RF}}_j\left(x^{\prime },y^{\prime },z^{\prime }\right)\right]$

where f_j,epmel, f_j,pmel, and f_j,rmel, are the wavelength dependent volume fractions of melanin in, respectively, epidermis, papillary dermis and reticular dermis; and f_j,poxy, f_j,roxy, f_j,pdeoxy, and f_j,rdeoxy are the wavelength dependent volume fractions of oxy-hemoglobin and de-oxy-hemoglobin in, respectively, papillary dermis and reticular dermis. Validation of the melanin and blood-flow volume information can be compared with the pathology for test cases for validation.

[0123] Experiments

[0124] Different visible light and near-infrared wavelengths in the range of 400 nm 1000 nm (preferably 500 nm to 960 nm) may be used to obtain diffused reflectance-based absorption measurements with different levels of depth of penetration to estimate volumes of selected chromophores such as melanin, oxy-hemoglobin and de-oxy-hemoglobin. A summary of the depth of penetration and characteristic absorption to chromorphores of interest is provided in Table 2.

TABLE 2
Specific light wavelengths and their properties for use in multispectral
skin imaging.
Wavelength	Optical Properties in Skin
405 nm	Superficial penetration (<0.5 mm), used in fluorescence excitation of ALA. High
	melanin absorption.
470 nm	Superficial penetration (<1.0 mm), used in fluorescence excitation of fluorescein,
	highly absorbed in melanin. Excellent for imaging superficial melanin associated
	with sun damage.
500 nm	Superficial penetration (<1.2 mm), used for fluorescence imaging
520 nm	Some fluorescence imaging (<1.5 mm penetration),
560 nm	Superficial penetration (<1.8 mm), highly absorbed by de-oxy-hemoglobin.
	Highlights superficial blood vessels.
580 nm	Deeper penetration in tissue (<2.3 mm), excellent for imaging most pigmented
	lesions without interference from the superficial melanin, peak for oxy-
	hemoglobin absorption.
610 nm	Deeper penetration in tissue (<2.5 mm), excellent for contrast imaging
	deoxygenated hemoglobin and veins due to large difference in absorption for oxy
	and de-oxy-hemoglobin.
660 nm	Deep penetration in tissue (>3.0 mm), shows deep pigmentation; large difference
	in absorption for oxy and de-oxy-hemoglobin.
760 nm	Very deep penetration in tissue, peak for de-oxy-hemoglobin and water
	absorption.
800 nm	Very Deep penetration, almost equal absorption for de-oxy and oxy-hemoglobin.
910 nm	Very Deep penetration in tissue, peak for oxy-hemoglobin absorption.

[0125] Optical wavelength of 560 nm provides superficial penetration (<1.8 mm) with high absorption by deoxy-hemoglobin and is suitable for imaging superficial blood vessels. While 610 nm wavelength provide deeper penetration in tissue (<2.5 mm), it is also excellent for contrast imaging deoxygenated hemoglobin and veins due to large difference in absorption for oxy and de-oxy-hemoglobin. The images obtained from 660 nm and 910 nm provide measurements for the blood oxygen saturation level (SO₂). Thus a combination of selected wavelengths from spectral measurements of specific chromophores can provide useful information about melanin, oxy- and deoxy-hemoglobin.

[0126] For measurements related to melanin, vascular absorption due to de-oxy-hemoglobin, and blood oxygen saturation level, 560 nm, 580 nm, 610 nm, 660 nm, 760 nm, 800 nm and 910 nm wavelengths are selected for imaging through a multispectral nevoscope apparatus 200. The images obtained from 660 nm, 760 nm, 800 and 910 nm can be used to estimate the blood oxygen saturation level (SO₂) of the skin-lesion. These measurements can be used in 3-D reconstruction of Hb and HbO₂ distributions using a Monte-Carlo simulation based reconstruction method as described herein. Other wavelengths (560 nm, 580 nm, and 610 nm) can be used for visual and 2D analysis of melanin and superficial-to-deep deoxy-hemoglobin-based vascular absorption.

[0127] Near-Infrared Imaging for Oxygen Saturation Level Measurement (HbO₂ and Hb Absorption)

[0128] In the absorption spectra of HbO₂ and Hb, as the oxygen saturation goes up, more of the light longer than 800 nm gets absorbed while the shorter wavelength light absorption decreases. To compare the localized absorption difference, we need to obtain the light spectrum with a specific optical path length, and accurately measure how much spectral distortion is caused by the blood absorption, for long and short wavelength range respectfully. In order to analyze the wavelength dependent absorption coefficients, the image obtained from the CCD camera for a fixed wavelength will be analyzed to estimate the absorption coefficients at various wavelength.

[0129] Generally the light signal level, r_s decays with increasing depth. The attenuation coefficient is the summation of scattering coefficient α_s and absorption coefficient α_α(k). To simplify this problem, we assume α_s is constant over the bandwidth of the light used. The cumulative spectral distortion at depth l is determined by the light absorption from the end of probe to l, which can be written as:

${\stackrel{\$}{r}}_s\left(k,l\right)=r_{s0}\exp \left\{-{\int }_0^l\left[{\alpha }_s\left(l^{\prime }\right)+{\alpha }_a\left(k,l^{\prime }\right)\right]l^{\prime }\right\}$

[0130] By averaging {circumflex over (r)}_s(k,l) in long and short wavelength range, the following are obtained:

$\begin{array}{cc}\begin{array}{c}{\stackrel{\$}{r}}_{\mathrm{s\_short}}\left(l\right)=r_{s0}\exp \left(-2{\int }_0^l{\alpha }_s\left(l^{\prime }\right)l^{\prime }\right){\left\{\stackrel{\_}{\exp \left[-{\int }_0^l{\alpha }_a\left(k,l^{\prime }\right)l^{\prime }\right]}\right\}}_{\mathrm{short}}\\ =r_{s0}\exp \left(-{\int }_0^l{\alpha }_s\left(l^{\prime }\right)l^{\prime }\right)A_{\mathrm{short}}\left(l\right)\end{array}& \\ \begin{array}{c}{\stackrel{\$}{r}}_{\mathrm{s\_long}}\left(l\right)=r_{s0}\exp \left(-2{\int }_0^l{\alpha }_s\left(l^{\prime }\right)l^{\prime }\right){\left\{\stackrel{\_}{\exp \left[-{\int }_0^l{\alpha }_a\left(k,l^{\prime }\right)l^{\prime }\right]}\right\}}_{\mathrm{long}}\\ =r_{s0}\exp \left(-{\int }_0^l{\alpha }_s\left(l^{\prime }\right)l^{\prime }\right)A_{\mathrm{long}}\left(l\right)\end{array}& \end{array}$

[0131] { }_short,long indicates that the average value of the quantity in side of the brackets in short (e.g. 660 nm) or long wavelength range (e.g. 910 nm).

[0132] Although the spatial dependence of absorption coefficient is not exactly known, A_short and A_short are related to the oxygen saturation level, because the blood absorption coefficient α_a(k,l) is a function of SO₂

α_α(k,l)=[SO₂(l)]α_HbO2(k)+[1−SO₂(l)]α_Hb(k)

where α_HbO2 is the wavelength dependent absorption coefficient of HbO₂ corresponding to wavevector k; α_Hb is the absorption coefficient of Hb. And SO₂ varies as l. R(l) defined as difference of absorption between long and short wavelength light can be computed as

$R\left(l\right)=\frac{{\stackrel{\$}{r}}_{\mathrm{s\_short}}\left(l\right)-{\stackrel{\$}{r}}_{\mathrm{s\_long}}\left(l\right)}{{\stackrel{\$}{r}}_{\mathrm{s\_short}}\left(l\right)+{\stackrel{\$}{r}}_{\mathrm{s\_long}}\left(l\right)}=\frac{A_{\mathrm{short}}\left(l\right)-A_{\mathrm{long}}\left(l\right)}{A_{\mathrm{long}}\left(l\right)+A_{\mathrm{short}}\left(l\right)}$

[0133] R(l) is not determined by the absolute amplitude of {circumflex over (r)}_s(k,l), but determined by the difference in absorption of the long and the short wavelength range. Therefore, R(l) does not show how much spectral absorption is caused by blood, but rather shows the SO₂ related difference in spectral absorption. By images obtained from the short and long wavelength, processing data with aforementioned algorithm, we can obtain 2D SO₂ image with R(l_x, l_z).

[0134] Computer Classification of Multispectral Nevoscope Images

[0135] The disclosed classification system uses an adaptive fuzzy c-means clustering (AFCM) algorithm to provide an optimal set of classifiable clusters of the correlated features extracted from the white-light and multispectral nevoscope images. Now referring to FIG. 12 a schematic diagram of a computer classification system using adaptive Fuzzy clustering and portioning is detailed.

[0136] Adaptive Fuzzy c-means Clustering Method

[0137] The fuzzy c-means clustering method utilizes an adaptable membership value that can be updated based using the distribution statistics of the data points assigned to the cluster minimizing the following objective function J_m(U, v)

$J_m\left(U,v\right)=\sum _{i=1}^c\sum _{j=1}^nu_{\mathrm{ij}}^md_{\mathrm{ij}}^2=\sum _{i=1}^c\sum _{j=1}^nu_{\mathrm{ij}}^m{x_j-v_i}^2$

[0138] where c is the number of clusters, n the number of data feature vectors, u_ij is the fuzzy membership, m is the fuzziness index, x_j; j=1 . . . , n is the input data feature vectors, v_i; i=1, . . . c are the cluster centroids; and d_ij is the Euclidean distance vector. Based on the constraints defined on the distribution statistics of the data points in the clusters, fuzziness index can be defined between 1 and a very large value (maximum allowable variance within a cluster). The membership values can be defined as:

0≦u_ij≦1 for all i,j

$\sum _{i=1}^cu_{\mathrm{ij}}=1\mathrm{for}\mathrm{all}j\mathrm{and}0<\sum _{j=1}^nu_{\mathrm{ij}}<n\mathrm{for}\mathrm{all}i$

[0139] An adaptive version of the FCM algorithm using a weighted derivative information is described by Pham and Prince, D. L. Pham, J. L. Prince, “Adaptive fuzzy segmentation of magnetic resonance images,” IEEE Trans on Med. Imaging, 18, 9, 737-752, (1999). The adaptive FCM algorithm provides a good regularization in the clustering process through the first- and second-order derivative terms that force the bias field to be smooth. As the objective function is minimized for clustering, the minimization of the derivative information provides smooth homogenous clusters.

[0140] The objective function J_AFCM(U, v) to be minimized for an adaptive fuzzy c-means clustering algorithm can be given as

$J_{\mathrm{ACFM}}\left(U,v\right)=\sum _{i=1}^c\sum _{j=1}^nu_{\mathrm{ij}}^m{x_j-g_jv_i}^2+{\lambda }_1\sum _{j=1}^n\sum _{r=1}^R{\left(D_r\otimes g\right)}_j^2+{\lambda }_2\sum _{j=1}^n\sum _{r=1}^R\sum _{s=1}^R{\left(D_r\otimes D_s\otimes g\right)}_j^2$

where g_j is the unknown gain field to be estimated during the iterative clustering process, m is the fuzzification index (m>1), R is the dimension of the image, D_r is the derivative operator along the r^th dimension of the image, and λ₁ and λ₂ are, respectively, the weights for the first-order and second-order derivative terms. Second and third terms in the objective function provide the smoothness constraints.

[0141] It should be noted that though there may be as many as 14 classes for the lesion classification, the number of clusters c may be much larger than 14 depending on the distribution of data in the feature space. All clusters are required to be homogenous with respect to each category and must satisfy the maximum allowable variance criterion. If a cluster violates the criterion when mapped with the pathology classification labels, it is split and a new value of c is computed. After the mapping and labeling of each cluster is completed, the class models are developed based on the set of clusters and their parameters including mean, variance and fuzzy membership functions. For example, a two-dimensional feature space is shown in FIG. 13A for simplicity to illustrate a hypothetical class distribution. Class A includes two clusters, 1 and 3, while the Class C includes clusters 4 and 6. The model for respective classes must include all clusters and their parameters. These models are stored in the database for each class and updated on a need basis whenever there is disagreement between the pathology and the computed classification from any case in the test set. Thus, the model adapts to the new knowledge and learns new model parameters re-organizing the clusters. The new learning is achieved through the competitive learning approach using a leader-follower clustering algorithm, see, e.g., A. K. Jain and R. C. Dubes, Algorithms for Clustering Data, Prentice Hall, Englewood Cliffs, N.J., (1998). The nearest clusters are polled for the fuzzy membership functions. Based on the “Winner-Take-All” strategy, the cluster with maximum fuzzy membership function is chosen for the inclusion of the new data vector as long as it meets the cluster homogeneity and variance criteria.

[0142] After the models for pathology classification classes are constructed through labeling of corresponding clusters and their parameters, an integrated competitive learning approach for the classification is developed. In one example hyperplanes as shown can be placed by straight lines in FIG. 13A. Hyperplanes can be placed in the feature space through checking the clusters in each model and creating convex sets around them. W. Grohman and A. P. Dhawan, “Fuzzy convex set based pattern classification of mammographic microcalcifications”, Pattern Recognition, vol. 34(7), pp. 119-132 (2001). The creation of convex sets is time consuming but if the cluster parameters are stored for each model class, their respective fuzzy membership functions can be used to partition the feature space efficiently. After the convex sets around the clusters are identified, hyperplanes can be placed for computing the fuzzy classification membership functions to allow final classification using a “Winner-Take-All” strategy as shown in FIG. 13B.

[0143] Classification

[0144] A fuzzy classification membership function M_ƒ would be constructed to reflect the true shape of the convex subset with the greatest precision possible. The function M_ƒ is defined for each subset ƒ (ƒ=1, 2, . . . , K) of the class model as follows:

$M_f\left(x\right)=\sqrt[L_f]{\prod _{i=1}^{L_f}{\theta }_i},{\theta }_i=\frac{1}{\left(1+{}^{{\lambda }_{\mathrm{if}}{\varphi }_ix}\right)}$

where: L_ƒ is the number of separating hyperplanes for the subset ƒ, φ_i is the i-th separating hyperplane function for the subset, in the vector form, x is the input vector in the augmented form, and λ_iƒ is the steepness (scaling) coefficient for the i-th hyperplane in the subset ƒ. The value of λ_iƒ depends on the depth of convex subset ƒ, as projected onto the separating hyperplane H_i (defined by φ_i):

${\lambda }_{\mathrm{if}}=\frac{-\log \left(\frac{1-\chi }{\chi }\right)}{{\mu }_{\mathrm{if}}},{\mu }_{\mathrm{if}}=\frac{1}{n}\sum _{j=1}^n{\varphi }_ix_j$

where: n is the number of training points in the covex subset ƒ, φ_i is the i-th hyperplane equation in the vector form, μ_iƒ is the depth of the convex subset ƒ, as projected onto i-th hyperplane, x_j is an augmented coordinate vector of the j-th point in the subset, and χ is the center value of the membership function.

[0145] The output CO of the classifier is the category C of the convex set fuzzy classification membership function M_i that attains the highest value for the specified input pattern x, i.e.:

$\left(\mathrm{CO}=C1\le \underset{f\ne i}{\overset{\forall }{f}}\le {\mathrm{KM}}_f\left(x\right)<M_i\left(x\right),M_i\in C\right)$

where: CO is the output of the classifier, K is the number of convex sets obtained during training (number of fuzzy function neurons in the fuzzy membership function layer), M_i is the highest fuzzy classification membership function value for the input x, and C is the model class of convex subset used to construct membership function M_i.

[0146] Applicants have attempted to disclose all embodiments and applications of the described subject matter that could be reasonably foreseen. However, there may be unforeseeable, insubstantial modifications that remain as equivalents. While the present invention has been described in conjunction with specific, exemplary embodiments thereof, it is evident that many alterations, modifications, and variations will be apparent to those skilled in the art in light of the foregoing description without departing from the spirit or scope of the present disclosure. Accordingly, the present disclosure is intended to embrace all such alterations, modifications, and variations of the above detailed description.

[0147] All references cited herein are incorporated fully by reference.