Doctorale École

Msii

Enagnon Aguenounon

Thèse Dirigée

M Gioux Sylvain

M Uhring Wilfried

M Gigan Sylvain

M Blondel

M Dadouche Foudil

M Dehghani

Michele Diana

Amir Nahas

Murielle Torregrossa

Manon Schmidt

Swapnesh Panigrahi

Joseph Angelo

Silvère Ségaud

Luca Baratelli

Agathe Marmin

Jason Smith

Julien Lamy

Jesse Schiffler

Amandine Elchinger

Sara Gravelyn

Traitement et visualisation temps réel dans le domaine fréquentiel spatio-temporel pour le guidage du geste chirurgical

Keywords: Diffuse optical imaging, Spatial frequency domain imaging, Multispectral imaging, Image-guided surgery, Real-time system, GPGPU

Optical imaging has a significant potential for image-guided surgery. It is quantitative, fast, lowcost, non-invasive, and non-ionizing. It is in this context that for my thesis I worked on the development of a real-time imaging system based on diffuse optics for surgery. The concept that we have developed revolves around the diffuse optical technique called SFDI for Spatial Frequency Domain Imaging, and more specifically its rapid acquisition implementation SSOP "Single Snapshot imaging of Optical Properties". In practice, the SFDI method uses the analysis of structured light patterns to extract the optical properties of biological tissue. Through multispectral characterization of these properties, it is then possible to quantitatively address physiological and structural information and thus provide information interpretable by the surgeon.

However, the SSOP implementation is subject to artifacts and lower resolution. In addition, the classic multispectral acquisition is time-consuming, oftentimes suffers from low signal-to-noise ratio and is affected by motion. These limitations have long prevented the methodology's ability to successfully translate into the operating room. Therefore, the objective of my thesis was to target the limits of the technique by designing a system capable of real-time acquisition, processing, and visualization (i.e. faster than 25 frames per second). To achieve this goal, two methods were proposed to improve the visual imaging quality and the accuracy of the state of the art SSOP. The first method is based on anisotropic filters design in the Fourier domain and the second on a deep learning approach using convolutional neural networks. Both methods were validated against SFDI in vivo and the deep learning routine was demonstrated to generate a maximum error rate of 9.5% for the extraction of optical properties with consideration of the 3D profile of the sample. Furthermore, we have also developed a real-time multispectral acquisition method based on the temporal modulation of light, which has been successfully validated on phantoms and in vivo.

Finally, to ensure real-time processing and visualization, all algorithms developed were implemented on GPU (Graphics Processing Unit) in low-level programming language C CUDA (Compute Unified Device Architecture). In conclusion, this work contributed to the design and validation of the first imaging system capable of visualizing oxygen saturation levels in real-time (faster than 25 images per second) across a wide field of view (17.5 x 17.5 cm²). Future work following this thesis could consist of increasing the number of wavelengths of the system to six to provide surgeons with an unprecedented accuracy and significantly more physiological

Mean error percentage value for absorption at each spatial frequency using the selected filters combinations. Results are given as % mean ±% standard deviation. Color-coding is used for ease of interpretation (scale on the right) …………………………….…………………………..74 Fig. 4.5. Mean error percentage value for reduced scattering at each spatial frequency using the selected filters combinations. Results are given as % mean ±% standard deviation. Color-coding is used for ease of interpretation (scale on the right) …………………………….………………75 Fig. 4.6. FWHM of a vein on the back of the hand using all filters combinations. Results are given as % mean ±% standard deviation. Color-coding is used for ease of interpretation (scale on the right) ……………………………………………….…………………………………………….76 Filter equations and parameters. For each filter presented, the filter equation is given as well as the definition of the cutoff frequency fc as a function of spatial frequency fx. Note that the filters are also scaled according to the field of view width d ………….………………….……….73 Table 6.1. Mean percentage error in absorption and reduced scattering at 665 nm and 860 nm compared to the reference MATLAB CPU implementation ……………………………………118 A1.3. Mean error percentage over 10 measured hands for absorption at 0.1 mm-1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right) ………………………………………………………….………………………..…..xxix Table A1.4. Mean error percentage over 10 measured hands for absorption at 0.2 mm-1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right) ………………………………………………………………...……………………...xxx Table A1.5. Mean error percentage over 10 measured hands for absorption at 0.3 mm-1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right) ………………………………………………………………………………….....…xxxi | List of tables xviii Table A1.6. Mean error percentage over 10 measured hands for absorption at 0.4 mm -1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right) …………………………………………...………………………………………….xxxii Table A1.7. Mean error percentage over 10 measured hands for reduced scattering at 0.1 mm -1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right) ……………………………………………………………..……………..xxxiii Table A1.8. Mean error percentage over 10 measured hands for reduced scattering at 0.2 mm -1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right) ……………………………………………………………………………xxxiv Table A1.9. Mean error percentage over 10 measured hands for reduced scattering at 0.3 mm-1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right) …………………………………………………………….………………xxxv Table A1.10. Mean error percentage over 10 measured hands for reduced scattering at 0.4 mm-1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right) ……………………………………………………………………………xxxvi Table A1.11. Mean error percentage over 10 measured hands for the measurement of vein full width at half-maximum (FWHM) size. Measurement performed at 0.2 mm-1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right) ……………………………………………………………………………………………….xxxvii Table A1.12. Mean error percentage over 10 measured hands for the measurement of vein full width at half-maximum (FWHM) size. Measurement performed at 0.3 mm-1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right) ……………………………………………………………………………......……………..xxxviii Table A.13. Mean error percentage over 10 measured hands for the measurement of vein full width at half-maximum (FWHM) size. Measurement performed at 0.4 mm-1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right) ..

Introduction

"It is health that is real wealth" said Mahatma Gandhi. Through the history of mankind, progress was truly achieved with fighting diseases and developing healthy living habits resulting in longer life expectancy and higher quality of life. Traditionally, healers from all communities have learned to use plants to treat various diseases and affections. If taking medication can cure many diseases today, it certainly does not solve everything yet. That is why along with development in medical care, the class of surgeons was created; doctors who study the human body and perform invasive operations to treat the disease from the inside. However, despite all technological progresses that are surrounding our lives, most surgical acts are still performed in a subjective way, with the outcome of the surgery depending greatly on the surgeon's skills and experience. In order to improve the outcome of surgeries, technological developments have led the creation of various imaging instruments for the human body, mainly used before the surgery for treatment planning.

Despite this progress, there is still a disparity in the quality of care due to the human nature of surgeons, the incapacity in most cases to "see" the disease and variation in individual anatomy. To overcome such limitations, intraoperative imaging instruments were created to guide surgical procedures. Yet, very few of these instruments allow surgeons to discern unhealthy tissues in situ and in real-time, to resect, or preserve structures of interest.

Recent research has shown that light in the near-infrared domain (from 650 nm to 1000 nm), invisible to the human eye, has the advantage of penetrating deeply into living tissue and providing structural and functional information on the state of these tissues. This type of information has a strong potential for the identification and characterization of tissue types as well as their condition directly during surgery. As such, specific diffuse optical measurement techniques have been developed to extract tissue properties. Fundamentally, when photons propagated in matter, they interact with the tissue constituents, through two major physical interaction phenomena called absorption and diffusion. These two parameters represent the optical characteristics of matter and are conventionally referred to as optical properties. In practice, the measurement of these optical properties using several wavelengths of light allows one to retrieve physiological information needed by surgeons: oxygenation, lipids, melanin, and water levels.

The first diffuse optical imaging technique of interest for our application is Spatial Frequency Domain Imaging (SFDI). This relatively inexpensive imaging method provides quantitative and wide field of view (> 100 cm 2 ) optical properties measuring capabilities. For this purpose, it uses patterns of light that are projected on the tissue. The images of these patterns are collected by a camera and then processed to extract the optical properties of the tissue. However, this method requires the acquisition of several images to trace physiological information and provide data that can be interpreted by the surgeon. This constraint results in a long acquisition time, making this method unsuitable for a real-time application such as surgery. In order to overcome this problem and speed up the process, a real-time acquisition method called Single Snapshot imaging of Optical Properties (SSOP) has been developed. The purpose of this method is to perform the same optical properties measurement than SFDI but by using a single pattern of light projected onto the tissue.

With this technique, the number of scans is reduced to the number of wavelengths necessary for the extraction of physiological information. This is however performed at the cost of a loss of image resolution, unwanted artifacts and a decrease in accuracy in the results.

In this direction, the objective of this doctoral thesis is to develop a new technology for assisting surgical gestures through imaging. Towards this aim, our first contribution presented herein is to combine the SSOP method with the temporal modulation of light in order to compute the optical properties of tissues at several wavelengths simultaneously. The use of this modulation technique should enable fast acquisition across several wavelengths with a slightly longer processing time due to the additional demodulation time step before extracting the optical properties. Hence, a major difficulty of our research is acquiring images using a fast video camera as well as a current lack of real-time SSOP processing techniques necessary for multi-wavelength analysis. In addition to this work, our second contribution is the development of a real-time processing methods for SSOP that ensure better visual quality of images while improving the accuracy of results.

This manuscript covers the scientific work achieved during the three years of the PhD degree and is organized as follows:

Chapter 1 briefly summarizes what are the techniques used currently in image-guided surgery.

We describe how the different imaging systems work and what are their advantages and drawbacks. At the end of this chapter we highlight the contribution we are making in the field of optical imaging for healthcare.

Chapter 2 presents the generalities of optical imaging used in near-infrared spectroscopy. We describe the optical properties phenomena in biological tissue, present the main models used for light propagation in tissue and show different techniques allowing to measure these optical properties.

Chapter 3 summarizes the limitations of the fast Single Snapshot imaging of Optical Properties (SSOP) used in the design of the real-time imaging system. In addition, a brief presentation of our contributions is introduced, leaving the full description to the following chapters of the thesis.

Chapter 4 presents filtering strategies for improving the image quality of optical properties maps obtained with the Single Snapshot imaging of Optical Properties method. To do so, we investigate the effect of anisotropic two-dimensional filtering windows. Both accuracy and image quality of the optical properties are quantified at the end of this chapter in comparison with standard, multiple image acquisitions in the spatial frequency domain.

Chapter 5 introduces real-time, wide-field, and quantitative blood parameters imaging using spatiotemporal modulation of light. The method we propose is presented and is validated by a phantom measurement experimentation and an in vivo experimentation. This chapter ends with the presentation of the limitations of the method.

Chapter 6 presents the General Purpose Graphic Processing Unit (GPGPU) implementation of the approach for real-time multispectral acquisition in the spatial frequency domain using spatiotemporal modulation of light introduced in Chapter 5. We present the comparison of Graphics Processing Units (GPU) and Central Processing Unit (CPU) processing implementations for this method with special emphasis on processing time. This chapter ends with prospective improvement needed for the design of the imaging system.

Chapter 7 presents a direct follow-up to the previous chapter by presenting a novel implementation of SSOP based on a combination of deep learning and GPU processing. The method requirements and implementation are described. A comparison with the multi acquisition SFDI method is also performed to assess its performances. This chapter concludes that simultaneous high-quality image reconstruction, with surface profile correction and accurate optical property (OP) extraction in real-time across large fields of view, is a reality.

The manuscript is concluded with a summary of the achievements obtained during this thesis work, followed by a discussion on the results and the perspectives for future improvements related to the project The design of a medical imaging system needs to be justified and as such, its advantages over existing systems are to be demonstrated. Indeed, the journey leading to the realization of the final product is long and strewn with ambush. In this chapter we provide a quick overview of current surgery image-guided systems in order to clearly emphasize the contribution we are making herein within the field of optical imaging.

Definition and objectives

Visualization and touch sensing are the primary interfaces between a surgeon and his patient.

However, for the human eye and touch, most tumors such as brain, colorectal or breast cancer are not distinguishable from surrounding normal tissues. As such, essential functional information on tissue temperature, tissue perfusion, or metabolism cannot be visualized by the human eye or well sensed by touch and therefore, they cannot be properly used during surgery. This lack of perception of what should be resected and what should be avoided leads to unacceptably high recurrence, morbidity and mortality rates. The purpose of image-guided surgery, also called intraoperative imaging, is to use imaging to enhance the surgeon natural senses. Fundamentally, image-guided surgery aims at improving the localization and targeting of diseased tissue, in order to accurately guide the surgeon to navigate towards it, as well as to monitor and control the patient treatments in real-time during both open and minimally invasive procedures. This need for guidance is even more pronounced during minimally invasive procedures where the target cannot be observed directly, and the touch sense is not possible. In addition, it is critical that image-guidance systems which must necessarily perform in real-time to provide adequate guidance while surgical instruments are manipulated by the operator. Intraoperative imaging must not be confounded with preoperative imaging which enables pre-surgical planning to understand the patient specific anatomy, the location of disease, and develop a specific treatment approach. To develop an appropriate image-guided surgery system, it is important to understand several aspects of the data acquisition, data processing, and visualization methodologies in the context of the interventional procedure. This is a challenging task given the different data acquisition modalities and their specificity. However, with the increasing number of users of image-guided technology and new clinical applications, the integration of these technologies will undoubtedly improve the efficacy of interventions and surgical procedures, expand surgeons' abilities to resect maximal amounts of tumor while preserving essential functions, reduce the risk of human error, and minimize practice disparities between hospitals and surgeons.

Conventional main image-guided system field

Since the discovery of X-rays in 1895 by Röntgen, a wide array of medical imaging modalities is in use currently, and nearly all of them have been integrated into some form of image-guidance system. In this section, they will be briefly described.

X-ray imaging and Computed tomography

The discovery of X-rays by Wilhelm Conrad Röntgen in 1895 represents the birth of medical imaging. For the first time, interior parts of a body could be visualized without cutting into it. The name X-rays was coined in Röntgen's original paper, where it denotes a new unknown ("X") type of radiation. An X-ray image characterizes the whole path pass of radiation through the body and produces a weighted average value of X-ray attenuation. The two dominant mechanisms for the interaction of X-rays with the tissue are photoelectric absorption and Compton scattering.

Photoelectric interactions in the body involve the energy of an incident X-ray being absorbed by an atom in tissue, with a tightly bound electron emitted from the K-or L-shell; the incident X-ray is completely absorbed and does not reach the detector. Compton scattering involves the transfer of a fraction of an incident X-ray's energy to a loosely bound outer shell of an atom in tissue. The Compton effect limits the contrast and quality of X-ray images [1].

Since its introduction, X-ray imaging has played a central role in image-guided surgery. The use of X-ray imaging is traditionally compartmentalized into three modalities:

• Radiography and Fluoroscopy (R&F): radiography refers to simple still 2D X-ray images and fluoroscopy to a special version of radiography where the result can be immediately represented on a screen providing real-time projection of images with high frame rates of up to 60 images per second. Typically, fluoroscopy images are acquired intraoperatively by a C-arm scanner.

• Rotational angiography: mostly used for angiography applications. Angiography is an imaging technique in which an X-ray image is acquired to display the inner open space (lumen) of blood-filled structures and this blood flow can be tracked through contrast agent-enhanced imaging. Since the rotation of the X-ray emitter/detector system is performed on a C-shape mounting, it is also called C-arm angiography and uses a full array of detectors, which allows the measurement of a full cone of rays.

• Computed tomography (CT) is an X-ray-based imaging method used to generate a 3D image volume of patient tissue densities. CT data represents a series of individual X-ray images that compose one volume dataset. These X-ray images are acquired by an emitter/detector system that rotates around the object to be scanned. The recent improvements in detector technology have enabled the development of CT scanners that acquire 256 or even 320 slices simultaneously or in dual source CT, where two detector units are placed in the scanner allowing to considerably increase the speed of CT image acquisition. These new capabilities have allowed to build intraoperative CT such as the Airo (Brainlab, Munich Germany) (Fig. 1.1.).

CT data are frequently acquired and represent a huge advantage over X-ray imaging in its ability to localize pathologies and to represent their shape. A major advantage of CT is that it is able to image bone, soft tissue, and blood vessels all at the same time. However, it is more expensive, and the in-plane resolution is rather low. As such, CT is not seen as a replacement of X-ray. CT data, in general, are rather easy to interpret, despite a few artifacts that need to be considered. Due to its wide availability and fast image acquisition, it is widely used, e.g. in emergency cases, to detect bleedings, to diagnose fractures and diseases related to skeletal structures [2]. A general issue related to X-rays is their potential to impair the genetic material, which may contribute to the development of cancer. Since both fluoroscopy and CT use X-rays for image formation, imageguided procedures using these modalities expose patients, physicians, nurses, technologists, and other staff to potentially hazardous radiation. In general, there is a trade-off between image quality (high SNR) and the dose of radiation. Despite this relation, during the last decades image quality has improved considerably while reducing the necessary radiation.

The most common applications for CT continue to be high-contrast imaging such as for treatment of intracranial aneurysms [3], evaluation of carotid artery stenosis [4], cerebrovascular stenting [5] and treatment of abdominal aortic aneurysms [6].

Magnetic resonance imaging

Magnetic resonance as an imaging modality in medicine was introduced by Lauterbur in 1973.

Lauterbur suggested that the dependency of the resonance frequency of protons on the magnetic field gradients could be used to determine the spatial distribution of protons in water (on a body) by different frequencies. To achieve this measurement, MRI requires the production of a strong, temporally stable, and spatially homogeneous magnetic field within the patient using superconductor technology. The most common field strength for clinical scanning is at least 1.5Tesla. Once protons are placed in the magnetic field, they have a resonant frequency that scales linearly with the magnetic field strength. When these protons are subjected to radiofrequency pulse at the resonant frequency, they absorb this energy and transit to a higher energy state. After the protons are stimulated, they slowly release the energy received. This relaxation is described as free induction decay and is divided into the transverse and longitudinal relaxation. The transversal relaxation is an exponential decay and the time required for a 63% decay compensation for brain shift [12], cardiovascular surgical interventions [13], biopsies guidance [14].

Ultrasound imaging

The first published use of ultrasound (US) as a medical imaging modality was in 1942 when Dr.

Karl Dussik measured transmission attenuation through the head to diagnose brain tumors [15].

Ultrasound as an imaging modality is based on soundwaves that are emitted at very high frequencies (1-15 MHz) from the ultrasound probe ("transducer") in direct contact with the respective body part. The sound waves penetrate human tissue at a speed of 1450-1580 m/s. When the sound waves travel through the medium, they cause periodic changes in density, pressure, and temperature. Indeed, the speed at which ultrasonic waves travel through medium change depending on the specific tissue constituents such as the percentage of protein, collagen, fat, and water. As they travel through different layers of tissue, they encounter different specific acoustic impedances, and therefore a certain fraction of the intensity of the wave is transmitted (viscous losses, heat conduction, and relaxation processes), with the remainder being reflected at the interface between the different tissues. The reflected sound waves are then recorded by sensors located next to the sound sources and used to generate the respective images.

Different ultrasound acquisition modes exist for different applications:

• A-mode: where single lines of pulse and echo ultrasound allow to define the depth of the reflected or backscattered signal.

• B-mode (Brightness mode): where by varying the excitation times of individual elements in a transducer array, an image of the region of interest is taken. This mode typically

shows anatomical details such as kidney, breast, and liver.

• M-mode: specific A-mode to measure motion (e.g. of the heart).

• Doppler-mode: specific A-mode used to measure velocities (e.g. in the carotid artery) using the well-known Doppler effect.

• 3D Ultrasound: allows three-dimensional image acquisition by acquiring and accumulating multiple scans, while tracking the movement of the ultrasound probe.

The first advantage of ultrasound is its real-time imaging capability which is important for intraoperative imaging. Ultrasound is also nonionizing, portable, and inexpensive compared to other clinical imaging modalities. An example of intraoperative US system is shown in Fig.

Ultrasound application areas include:

cardiac and vascular applications for imaging blood flow, ventricular wall motion, valve opening/closing and artery plaque measurements [16] | Chapter 1. Image-guided surgery

14

orthopedic applications for intraoperative imaging of bone or real-time needle guidance [17] neurosurgery applications for intraoperative tumor imaging [18] traditional surgical interventions as a navigation tool and in guided biopsies, such as biopsy of breast tumors [19] abdominal tumor diagnosis and surgeries [20] obstetrics and gynecology [21] Fig. 1.4. BK ultrasound intraoperative system Flex Focus 800 and a few acquired images.

Imaging in nuclear medicine

In contrast to X-ray, ultrasound and MRI, nuclear medicine imaging techniques do not produce an anatomical map of the body but instead, They image the spatial distribution of radioactive materials (radiotracers) introduced into the body in order to detect early biochemical indicators of disease. The radioactive tracer is injected into the patient via inhalation, directly into the bloodstream, or by oral administration. Due to its biochemical properties, the tracer may be distributed in proportion to glucose metabolism or bind to a particular molecule preferentially expressed on cancer cells. Abnormal tissue distribution, an increase or decrease in the rate of the radiotracers are a strong indicator of disease. Upon decay of the radioactive atom attached to the tracer, either a positron or gamma ray is emitted.

In classic planar nuclear medicine (2D), radiation in the form of gamma-rays is detected using an imaging device called a "gamma camera". The gamma camera is based on a large scintillation crystal that transduces the energy of a gamma-ray into light at 415 nm. The intensity of the light is proportional to the energy of the incident gamma-ray.

More advanced nuclear medicine imaging includes two major techniques named Single Photon Emission Computed Tomography (SPECT) and Positron Emission Tomography (PET):

• SPECT is based on gamma cameras, where two or three gamma cameras are rotated around the patient, acquiring multiple images of the 3D distribution of a radiotracer. To compensate for poor anatomical differentiation, SPECT can be combined with CT to form hybrid scanners SPECT/CT that provide truly integrated information of metabolic activity and anatomical properties.

• PET scans generate 3D images of functional processes, such as metabolic activity by using a short-lived radiotracer injected and traced throughout the body to its target area. At the target area, the substance is processed by the metabolism and the radioactive isotope decays and emits positrons. These positrons travel a short distance in the tissue before annihilating with an electron, resulting in the formation of two gamma-rays. The two gamma-rays travel in opposite directions (180°) and are detected by a ring of detectors. The radiotracer spatial distribution can then be used to distinguish diseased from healthy tissues. Like SPECT, PET can be combined with CT or MRI modality to form hybrid scanners such as PET/CT or PET/MRI. An example of such a system is shown in Fig. 1.5..

Both SPECT and PET present the same risks that occur with X-ray imaging because of the radioactive materials they use. Both are designed mostly as diagnostic imaging systems since neither of them is a real-time imaging modality, the number of photons emitted from the patient per second being significantly lower. SPECT can be applied to diagnosis in oncology [22], bone imaging [23], neuroimaging [24], and the investigation of myocardial perfusion [25]. The most common clinical applications of PET are: oncology [26], tumor detection [27], epilepsy treatment | Chapter 1. Image-guided surgery 16 [28], cardiology [29] (to map heart function and contribute to the diagnosis of atherosclerosis) and neurology [30] (to represent brain activity and thus to contribute to the diagnosis of neurodegenerative diseases). Even if they are rarely used for direct image guidance, they can however be used in surgery planning such as a preoperative imaging technique or have their images integrated to a conventional image-guided vision system to enhance its capabilities, such as for hybrid scanners. For example, integrate PET images into image-guided resection PET-/CT-guided cryoablations of soft tissue tumors [31], used intraoperative PET/CT imaging in tumor localization and verification of complete tumor resection in breast cancer [32].

Optical imaging

The human eye can detect information in visible light described through color categories. This color depends on both the physics of the object in its environment and the characteristics of the perceiving eye. Physically, objects can be said to have the color of the light leaving their surfaces, which depends on the spectrum of the incident illumination and the reflectance properties of the surface but also on the light transmitted through or emitted by the object. Medical optical imaging uses the same concept on tissue to enhance limited surgeons' vision. Since many diseases involve changes in vasculature or hemodynamics, light can be used in tissue to infer structural and functional properties. Optical imaging describes various imaging techniques using visible, ultraviolet, and infrared light. In contrast with the conventional imaging modalities, it presents interesting features such as the use of nonionizing radiation, portability, high resolution, high sensitivity, and cost efficiency. However, depth penetration cannot practically reach beyond 1-2 cm due to tissue optical properties. The field of optical imaging includes two main categories of techniques: microscopic and macroscopic. Some of the most widely used techniques are presented in the following section.

White light imaging

White light imaging consists in imaging the reflected white light properties of the tissue. A simple image of the tissue is then provided to the surgeon to help him see usually invisible tissue information. The most common interventional use of white light is for diagnostic and therapeutic endoscopy and laparoscopy, where a camera enables less invasive procedures by accessing the anatomy through natural orifices or through small ports. Another common usage is for procedures that require significant zoom into the surgical area and in which a microscope is used to perform microsurgical tasks.

Multispectral imaging

Multispectral imaging is the extension of white light imaging. It incorporates better spectral resolution resulting in improved ability to distinguish different constituents in the tissue as shown in Fig. 1.6.. For example, it is used in flexible endoscopy to improve the detection of adenomas or hyperplasia through better imaging of the pit pattern or microvascular tortuosity [33] by a technique called narrow band imaging (NBI). There are many other applications of multispectral imaging, such as monitoring of diabetic feet, imaging of tissue ischaemia during transplant surgeries or during tissue monitoring prior to anastomosis, characterization of skin bruises, detection of cancerous changes in the prostate, stomach, breast, etc [34].

Microscopy techniques

The techniques already presented focused more on large area surveillance with wide-field imaging, since they are well-matched with the scale of many surgical interventions. However, surgical use of microscopy dates to the beginning of the twentieth century, originally employed to aid the otolaryngologists, ophthalmologists, as well as vascular, neurological, and plastic surgeons to perform precise small-scale tissue manipulation by simple lens magnification [35]. Nowadays, microscopy is also used at early stages to detect changes in the biochemistry and tissue architecture, thus allowing to reduce healthcare cost (e.g. Fig.

Augmented reality

In basic image-guided surgery, it is common to display images of the surgical or therapeutic environment based on preoperative or intraoperative images on a screen. A link between the patient and the medical images may sometimes be created through image registration and tracking tools.

However, such solutions force the surgeons to look away from the surgical site, thus resulting in a psychophysical decoupling of the actual and virtual therapeutic working spaces. To improve this cumbersome situation, virtual reality techniques are used to integrate reality and virtuality into an Augmented Reality (AR) system [36]. The goal of this combination, is to improve the surgeon's senses and abilities, therefore enhancing his hand-eye coordination. Using augmented reality, preoperative or intraoperative patient image (e.g., PET, CT, MRI, US) appears at the location of the actual patient anatomy, giving the surgeon the perception that he is "seeing through" the patient or organ. Moreover, the used of AR improves work efficiency, safety, cost, and patient outcome.

To achieve these goals, the virtual and the real patient image need to be registered with extreme precision and the virtual representation needs to be interactive and should be updated in real-time while performing a surgical or interventional procedure. If these requirements are not satisfied, the AR system could become dangerous, as the surgeons would be guided to place instruments in incorrect locations. Thus, medical AR requires a complete alignment pipeline, where instruments are tracked, cameras need to be calibrated, patient data need to be registered, display and visualization systems must be adapted for the surgical task. Applications targeted by AR systems include neurosurgery and otolaryngology, cranio-and maxillofacial surgery, breast and abdominal needle biopsies and tumor ablations, orthopedics, and cardiovascular and thoracic surgery [37].

An example of an AR surgery product which received a U.S.

Impact for surgeons

Generally, image-guided surgery improves patient outcome as expected, compared to simple surgery. Image-guided surgery also has many benefits for surgeons. The use of image-guided surgery systems reduces surgery duration, improves the surgeon's senses, abilities, and performances. On the other hand however, surgeons need a training period, preparation time before each surgery is increased, and, for systems that directly intervene in the workflow of surgeons, a physiological effort may be needed at the beginning [39]. The use of image-guided systems may also affect surgeons' skills, thus making them dependent on imaging technology.

Conclusion

This chapter summarizes the currently used techniques in image-guided surgery, starting from their fundamental operating characteristics and concluding with their advantages and drawbacks.

Our point of view is that optical imaging has the best tradeoff compare to other modalities and further research in this field will bring at low cost a considerable difference in healthcare quality.

In accordance with the observations, our scientific contribution is related to macroscopic optical imaging. The objective of our project is to develop a wide-field (> 100 cm2), real-time (25 images per second), quantitative intraoperative image-guided surgery device based on the diffuse reflectance imaging. We propose to use SFDI method, which will allow the extraction of optical properties on a wide field of view using a structured pattern of light. The method is quantitative unlike standard white light imaging. Then perform a multispectral acquisition with the system, in an innovative way using temporal modulation of light to overcome classical limits of multispectral systems. The measurement of the optical properties at several wavelengths of light will allow to retrieve the physiological and structural parameters interpretable by surgeons such as oxygenation, lipids, melanin, and water levels, as well as microstructure density. Real-time being mandatory for image-guided surgery system, special care will be made on the implementation, by choosing low computation algorithm combined with high computation power. More details on the SFDI method, material and its applications are given in Chapter 2. Furthermore, in Chapter 3, the scientific issues related to our project and the specific objectives of our work are presented.

The next chapter will cover near-infrared optical imaging. We will describe the optical properties of the tissues, the different ways to model the propagation of light in a medium and the techniques of measuring the optical properties in the near-infrared spectral range.

Optical imaging brings at low cost numerous advantages in image-guided surgery. This chapter presents the main concept of optical imaging used in near-infrared spectroscopy. Starting by the description of light tissue interaction phenomena and how they are modeled, we present the three main near-infrared techniques used to measure biological tissue optical properties, before describing two widely used near-infrared imaging systems.

2.

Optical properties of biological tissues

Reflection, refraction, and refractive index

The first light-matter interactions to be considered to analyze optical properties are reflection and refraction, which occur where an incident light beam enters the biological tissue, and more generally at the interface between two media with distinct refractive indices (Fig. 2.1.). The refractive index n of a medium is defined as the ratio between the speed of light in vacuum and its speed in the medium. In water, n=1.33, whereas it can reach 1.5 in pure fat. In biological tissues, the refractive index can range from 1.35 to 1.45 [1]. Reflection of photons occurs due to a change of refractive index n between two media which vary according to temperature and wavelength. This part of the incident beam energy which is not transmitted through the interface but is reflected, is referred also as specular reflection. Refraction is the phenomena by which photons are transmitted in the medium with a change in forward speed and direction. The fractions of reflected and refracted photons are determined by the medium reflectivity and transmittance through the Fresnel equations and the direction they take after the

n i n r θ r θ i Specular reflection Refraction Surface θ i

Incident beam

collision is described by the Snell-Descartes Law [2] according to their incident angle θi and the refractive index ratio of the incident medium ni and refractive medium nr.

Snell-Descartes Law:

( ) ()

sin sin i ri r n n qq =× (2.1)

Absorption

Light interaction with biological tissues can result of photon absorption. Absorption of the light photons is a consequence of the partial transition of the electromagnetic light energy into another form of energy, such as thermal motion or vibration of the molecules [3]. In fact, when light travels in biological tissue, since the intensity of the incident electromagnetic wave is attenuated while passing through the medium, it allows to characterize the absorption capacity of the medium.

a dl dl I ml =- (2.2)
where is the light intensity and is the path length, leading to the well-known Beer-Lambert law:

() 0 e a l II ml - = (2.3)
where ( ) a ml is the absorption coefficient which depends on the wavelength λ and is expressed in

cm -1 .

Scattering

The scattering of light in biological tissues is a phenomena of reflection that occurs when photons encounter a particle (scatterer) with a refractive index n different from that of the biological entities (cells, nuclei, membranes), thus producing a deviation of the photon incident direction [3].

Scattered photons that travel through the medium to the opposite side are quasi-ballistic. Photons that experience multiple scattering event through the tissue and come back to the surface contribute to the Diffuse Reflectance (Fig. 2

.3.).

Fig. 2.3. Scattering in a medium.

There are different types of scattering regimes, depending on the relative size of the scattering object with respect to the wavelength of the incident wave (Fig. 2

.4.):

• The Rayleigh scattering describes the isotropic scattering of light by particles significantly smaller than the wavelength.

• The Mie scattering describes the regime in which the sizes of particles are comparable to the wavelength and corresponds to mainly anisotropic forward scattering.

• Finally, for objects larger than the wavelength, the laws of diffraction can be used.

I 0 I dl

The scattering coefficient is determined similarly to the absorption coefficient. If we consider a homogenous non-absorbing medium with a thickness l, the scattering coefficient links the measured output intensity I for a given incident collimated beam of wavelength λ and intensity I0:

() 0 e s l II ml - = (2.4) 1 s
m is then the mean path of the photon before undergoing a scattering event.

The scattering coefficient s m perfectly describes isotropic scattering. However, most biological tissues have a preferential scattering direction and tend to scatter light forward. To describe this, the phase function µ ( )

,' f ss $ is introduced: it gives the probability for a given photon of an incident direction s $ to be scattered in the direction µ ' s . This problem is generally simplified by assuming that the phase function only depends on the angle between s $ and µ ' s (Fig. 2.5.). The anisotropy of scattering g is then defined as the mean value of () cos q weighted by the phase function f:

() ( ) () () ( ) 1 1 cos cos cos g fd q qq - = ò (2.5)
In an isotropic medium, g=0, whereas g=1 for pure forward scattering. Typically, the value of g ranges from 0.7 to 0.99 in biological tissues [5].

The directional effects of scattering can be included in the reduced scattering coefficient which can be interpreted as an equivalent isotropic scattering coefficient:

( ) ' 1 ss g mm =- (2.6)

Fluorescence

Fluorescence is a phenomenon in which a molecule emits light when it returns to the ground state after being excited by an external light source [6]. The light produced by the molecule during relaxing to the ground state is of a longer wavelength relative to the wavelength of the excitation light (Fig. 2.6.). The fluorescence process can be divided into three steps: (1) excitation of molecule to new state by an external light source of specific wavelength within certain femtoseconds, (2) vibrational relaxation from the excited state for certain picoseconds and (3) molecule emitting light photons at a longer wavelength while returning to the ground state, which requires a few nanoseconds [6]. Fluorescence can be endogenous or exogenous. Endogenous fluorescence also called autofluorescence is the emission of light from naturally occurring fluorescent biomolecules due to particular molecules that have absorbed light. In contrast, exogenous fluorescence involves artificially added fluorescent markers called fluorescent dyes (also referred to as fluorescent probes). These contrast agents are molecules that can be targeted or non-targeted, and are conventionally used in the visible and/or the NIR spectral range [7]. A fluorophore can be characterized by three parameters:

• Fluorescence cross section corresponds to the absorption capacity of the fluorophore at a given wavelength.

• Quantum yield is the ratio of the number of photons emitted Ne to the number of photons absorbed Na. Quantum efficiency is always less than 1 and considerably variable depending on the used molecules.

• Lifetime is the duration of time for which the molecule is in the excited state. Let us consider a tub, containing fluorophores which are illuminated by a light pulse. The decrease in the number of excited fluorophores is given by:

() 0 e t N tN t aeö - ç÷ èø = (2.8)
where N0 is the number of photons excited at the initial instant and t the time constant.

The entire molecular fluorescence lifetime from excitation to relaxing to ground state is measured in nanoseconds, mostly between 0.5 and 20 ns. Such phenomenon is used in a number of ways, from fluorescent microscopic imaging devices to optical imaging of biological tissues.

Functional imaging

At the molecular level, the absorption phenomenon is caused by the medium composition in endogenous chromophores (specific light absorbing molecules) (Fig.

2.7.).

Each chromophore is characterized by its molar extinction coefficient in mol -1 .L.cm -1 which links its absorption coefficient to its concentration C in mol.L -1 in the medium. Therefore, according to the law of mixtures, the absorption coefficient at a single wavelength i l is given by:

() () 1 M a i m im m c ml el = = å (2.9)
where ( ) ai ml is the measured optical absorption coefficient, ( ) mi el is the chromophore extinction coefficient at i l , m c is the chromophore concentration, and M is the number of chromophores [9]. In order to determine the functional parameters of the sample, it is necessary to separate the contributions of the intrinsic biological chromophores.

These chromophores can be categorized into those that exhibit oxygen-dependent absorption and those whose absorption is essentially fixed during the course of a clinical measurement. Oxygendependent chromophores are the compounds of primary interest in tissue NIR spectroscopy, as their absorption spectra are directly influenced by the presence of oxygen. Oxygen-independent chromophores are compounds that display significant absorption in the NIR although their absorption spectra do not vary significantly as the blood/tissue oxygen level varies over a physiologically reasonable range. Their presence can be ignored when using a trend or quantified trend monitor. However, their presence must be accounted for when extracting measurements of tissue oxygenation from absolute measurements of a m at various wavelengths. The main chromophores of human tissues are:

• Hemoglobin. Hemoglobin is a complex metalloprotein (molecular weight = 64,500) that consists of four heme groups. Each of these heme groups can bind one molecule of O2 by undergoing a conformal change that causes the O2 molecule to become physically bound.

This conformal change shifts the absorption band exhibited by the Fe atom within each heme group, producing the observed spectral difference between deoxyhemoglobin (HHb)

and oxyhemoglobin (HbO2). Considering that these two forms have different absorption spectra in the NIR range, it is possible to deduce their relative concentrations from a measurement at multiple wavelengths. This is the working principle of commercially available pulse oximeters. A measure of oxygenation of blood is given by: ()

2 2 2
% 100 HbO StO HbO HHb =´+ (2.10)

• Cytochrome. Cytochrome aa3 is an enzyme located within the cell mitochondria that catalyzes the final step of the respiratory chain in which oxygen is reduced to water. The enzyme contains two copper centers, denoted as CuA and CuB. The CuA center is responsible for the distinctive broad absorption spectrum centered at approximately 830 nm when the enzyme is oxidized. This band disappears when the enzyme is reduced, therefore measuring the amplitude of this band can, theoretically allow one to determine the proportion of the enzyme that is oxidized relative to the proportion that is reduced (the "redox" state). This potential application has received enormous attention because it would allow for tracking oxygen utilization directly to the site of adenosine triphosphate ATP synthesis [7].

• Water. Water accounts for more than 80% of the volume of the human body and is a uniquely versatile solvent in which a host of biochemical reactions takes place. For this reason, water is considered to be one of the most important chromophores in tissue spectroscopy measurements. Nevertheless, for the purposes of most clinical measurements, the water concentration in tissue can be thought of as constant, and as such water acts as a fixed constant absorber. This absorption remains low for the wavelengths in the band ranging from Ultra-Violet (UV) to Near-Infrared (NIR).

• Lipids. Lipids are found throughout the body, particularly in adipose tissue. The absorption spectrum of lipids is similar to that of water. However, its overall contribution to absorption is relatively small. Lipids are considered a fixed constant absorber with normally unchanging concentration during clinical measurements and therefore, the measurements of changes in attenuation are not affected. These molecules strongly absorb light above 900 nm.

• Melanin. Melanin is a pigment that is responsible for the coloring of several biological tissues. Its role is to protect the epidermis and hair. It is strongly absorbing in the UV and in the visible. In addition, it can be broken down into two forms, namely eumelanin and pheomelanin [9].

Structural imaging

Optical scattering originates from light interaction with biological structures, which range from cell membranes to whole cells (Fig. 2.8.). Optical scattering can be described either as scattering by particles that have a refractive index different from the surrounding medium, or as scattering by a medium with a continuous but fluctuating refractive index. Photons are thus scattered differently according to the biological structure size. The scattering of a medium can be described as a contribution of Rayleigh scattering and Mie scattering. The review of tissue scattering properties as a function of wavelength which is not exhaustive but sufficient to characterize the behavior of tissues shows that scattering follows a power-law in biological tissue in the visible to near-infrared region [9]. Hence, data on a minimum of three wavelengths is sufficient to allow prediction of scattering at all wavelengths in the UV, visible or near-IR range. From this study two power-law were derived: where the wavelength dependence of scattering is described in terms of the separate contributions of Rayleigh and Mie scattering.

In these equations, the parameters are a, defined as the "scattering amplitude" and is representative of the density number of scatterers, and b(or bMie), the "scattering power" which represents the average size of the scattering particles. One interest of measuring scattering amplitude and scattering power is to use changes in the structure of cells and tissues as a contrast parameter while imaging cells or tissues macroscopically. Such analyses may prove useful in imaging the margins of cancers. Fig. 2.8. Biological structures of various sizes for photon scattering [11].

Physical models for light propagation in diffusive media

Radiative transfer equation

The radiative transfer equation (RTE), also called the Boltzmann equation, is an integraldifferential equation that models the transport of energy in an absorbing and scattering medium.

Following the principle of conservation of energy, it considers the energy transfer in a medium by superposition of energy flux (Fig. 2.9.). Before describing the equation, the radiance ( )

,, L rst r $
quantity must first be explained. It is defined as the average power that, at position r r and time t, flows through the unit area oriented in the direction of the unit vector s $ , due to photons within a unit frequency band centered at ν, that are moving within the unit solid angle around s $ , and is

measured in SI unit in (W•sr -1 •m -2 ).
The RTE with time-dependent sources is written as:

( ) ( ) () ()( ) () µ ( ) µ ( )( ) 4 ,, 1 ,, ,, ,' ,' , ' ,, a ss

L rst s Lrst r r Lrst r f ss Lrs t d S rst ct

p m mm ¶ éù =-×Ñ - + + W+ ëû ¶ ò r $ r r r r rrr $ $ $ $$ (2.13)
where c is the speed of light inside the medium, ' dW is the elementary solid angle in the direction µ ' s . The definition of each term involved in this equation are: : radiance gains from source.

o ( ) ,, 1 L rst ct ¶ ¶ r $ : temporal variations of radiance, o ( ) ,, s L rst -×Ñ r $$ : radiance lost from divergence, o () () ( )
As shown in the equation, solving analytically RTE turns out to be extremely complex in most cases. Therefore, a simplification can be performed or a stochastic approach can be used.

Diffusion approximation equation

The diffusion approximation is a popular approach of the transport model. This method allows the derivation of a simplified version of the RTE after performing successive approximations, leading to the diffusion equation (DE), which gives rise to simple analytical solutions. In the following, we will not enter the details of the derivation of the DE from e RTE, but we will only give the formulation of the DE and specify the various approximations considered to deduct the DE from the RTE. The RTE is an integro-differential equation for the radiance while the DE is a partialdifferential equation for the fluence rate f (w/m -2 ) defined as the integral of the radiance over the solid angle 4 :

( ) ( ) 4 , ,, rt Lrstd p f =W ò rr $ (2.14)
The crucial approximations which are made to deduct the DE from the RTE are the following:

-highly scattering medium:

( ) a r m r << ( ) s r m r
-the light source is mainly isotropic which means that the non-isotropic terms are not modeled accurately -the phase function

( ) ( ) ( ) ( )( ) ( )( ) , 1 , ,, a rt D r rt r rt S rt ct f f mf ¶ -Ñ × Ñ += ¶ r ur r ur r r r r (2.15) with:
-( ) (2.17)

,
In a semi-infinite homogenous medium, the solutions can be expressed as: where, z is the distance normal to the boundary and

( ) () ( ) () ( ) ( ) () ( ) (
( ) ( ) ( ) 1 ' 0 as z rr mm - =+ rr
For a more complex shape and heterogeneous medium, with layers and inclusions, the diffusion approximation equation can be solved with numerical methods such as the finite-elements methods, finite difference or finite volumes. In general, one drawback to the use of the standard diffusion approximation is that the assumption of ( ) a r m r << ( ) s r m r does not always hold for all biological tissues [12].

Monte Carlo Model

Monte Carlo refers to a technique first proposed by Metropolis and Ulam to simulate physical processes using a statistical sampling of the phenomena being modeled [13]. Introduced in the field of biomedical optics and photon migration by Wilson and Adam [14], the Monte Carlo (MC) simulation has become the gold standard for modeling light propagation in tissue optics. It is used to calculate the theoretical diffuse reflectance measured by the system, considering the system geometry and skin optical properties. In practice, his algorithm (Fig. 2.10.) works as follows. Once launched, a photon is moved by a distance Δs, a step size determined probabilistically through the mean free path of the photon depending on absorption and scattering properties. During this propagation, a portion of the weight is lost through absorption (Eq. 2.3.), after which the photon will be scattered with the new direction calculated probabilistically. If the photon encounters an interface it will be partially internally reflected, and partially transmitted out of the tissue (Fresnel reflection). The photon is repeatedly moved until its weight reaches a value that is low enough.

Note that to have no loss of energy, once reaching a low weight, 1 in n photons will continue propagating with a weight multiplied by n (roulette). This process is repeated for several million photons to ensure statistical accuracy of the simulated light transport in skin [15]. based simulations that utilize parallel computing to speed up the simulation time [18]. Due to its flexibility and recent advances in speed, the MC method has been explored to solve both the forward and inverse problems in tissue optics.

Monte Carlo simulations have also been used to generate lookup tables [19]. The lookup table method of determining optical properties involves the generation of a large table of predetermined diffuse reflectance values mapped to known optical properties. The diffuse reflectance measured on a sample can then be searched for or interpolated from the lookup table, along with the corresponding sample optical properties. The primary advantage of using a lookup table is that optical property determination is quite fast, given that the interpolation used is not computationally intensive. Since our goal is to provide real-time processing, we choose the Monte Carlo-based lookup table method for our work.

NIR techniques for biological tissues

Three different methods of using NIR light are currently used:

• continuous wave method

• time-of-flight method (also known as time-domain or time-resolved)

• frequency domain method

The continuous wave method has a rapid response but registers relative change only; it is therefore not possible to make absolute measurements using this technique. The time-of-flight method needs extensive data processing but provides more accurate measurements. It enables one to explore different information provided by the measurement. The third approach, which uses frequency domain or phase modulation technology has a lower resolution than that of the time-of-flight method.

Continuous Wave Near-Infrared Spectroscopy

The term "continuous wave" (CW) refers to the property of the instrumentation to be solely based on a light intensity measurement, i.e. near-infrared light is sent into the tissue and the intensity of the re-emerging (i.e. diffusely reflected) light is measured (Fig.

Frequency Domain Near-Infrared Spectroscopy

Frequency domain spectrometers, also called frequency-resolved or intensity-modulated systems, are based on the modulation of a laser light at given frequencies typically from 100 to 500 MHz.

They generally employ a laser diode or modulated white light sources. They measure the attenuation, phase shift, and modulation depth of the exiting light by either a photon-counting detector or gain-modulated area detector (Fig. 2.13.) [20,21]. Frequency domain instruments determine the absorption coefficient and the reduce scattering coefficient of the tissue by measuring alternating current (AC), direct current (DC), and phase change as functions of the distance through the tissue [22]. Moreover, the phase and amplitude shifts can be used for localization of the signal. Frequency domain spectrometers measure changes in intensity, phase, and modulation by using either (1) single wavelength and a fixed interoptode distance, (2) multiple wavelengths and a fixed interoptode distance, (3) or single wavelength and multiple interoptode distances [20]. Problems intrinsic to frequency domain spectrometers include noise and leakage associated with the high-frequency signal. In order to scan the entire wavelength spectrum of interest, multispectral acquisition is performed using individual light sources at particular wavelengths, such as LEDs or laser diodes, using a filter wheel to acquire each wavelength, using optical tunable filters, or using multispectral cameras which are however expensive and reduce image resolution [24]. Based on the previous description, few commercial devices are available on the market (Fig.

Diffuse Reflectance Imaging: SFDI

2.16.

). An example is the Reflect RS™ system; a research platform system with flexible 51 acquisition, processing, and visualization software commercialized by Modulim. The system architecture guarantees rapid sequential acquisition of SFDI data at multiple wavelengths with coregistered color images and is used by researchers in academia and industry to measure turbid systems. More recently, Modulim also released the Clarifi Imaging system for the measurement of tissue oxygenation in patients to detect possible compromised circulation before visible complications occur. Other systems that have been developed include a custom surgical guidance system that has been translated to first-in-human pilot studies [25].

Fig. 2.16. Illustration of SFDI commercial devices [25].

Over the last decades, the capabilities of SFDI have been demonstrated across a large number of biomedical applications including small animal imaging, burns wound assessment, skin flap monitoring, early stage diabetic foot ulcers diagnosis, among others [26][27][28][29][30]. For example, the use of SFDI in assessing burn wounds has shown that the change in water concentration, deoxyhemoglobin concentration, and optical scattering immediately postburn differed depending on the type of burn (superficial vs. partial thickness). In Port Wine Stain treatment with pulsed lasers, measurement using SFDI, allow the practitioner to quantify changes in PWS after treatment by determining physiological parameter changes [24].

Conclusion

In this chapter, we summarize important background knowledge that will support the following chapters. Optical properties in biological tissue were described and their use in order to obtain functional and structural parameters was presented. We also detailed the main models used for the description of light propagation in biological tissues, together with their application for optical properties estimation. The different approaches for measuring these optical properties of tissue were presented, as well as fluorescence and diffuse reflectance imaging capabilities. In the next chapter, the main challenges limiting diffuse reflectance imaging for image-guided surgery will be discussed, and our plan to overcome those limitations will be outlined. In Chapter 1, a quick overview of imaging methods for surgery was given. We highlighted the advantages of optical imaging over other methods. Then in Chapter 2, optical imaging in the nearinfrared range was presented in detail. Emphasis was placed on the theory of photon propagation and on imaging methods such as fluorescence and diffuse reflectance imaging. In this chapter, we focus on diffuse reflectance imaging, by presenting the limitations of the method called SFDI as part of an image-guided surgery system development. We also review the fast acquisition SFDI method's limitations and present the objectives of our work.

3.

General purpose of the project

The aim of our project is to develop a new intraoperative image-guided surgery device based on the diffuse reflectance imaging method SFDI. To achieve our goal, we plan to combine the fast acquisition method named SSOP for 'Single Snapshot imaging of Optical Properties' with temporal modulation of light, in order to discriminate optical properties of tissues acquired at several wavelengths simultaneously. The measurement of these optical properties at several wavelengths of light will permit to retrieve the physiological and structural information needed by surgeons: oxygenation, lipids, melanin, water levels, scattering power and scattering amplitude.

However, the use of this modulation technique will add a longer processing time, due to the demodulation time step, in addition to the optical properties estimation process implementation, which is not performed in real time at present. It is therefore mandatory to develop innovative processing techniques that guarantee real time acquisition, processing, and visualization.

SFDI/SSOP method limitations

Presented in the previous chapter, SFDI provides a quantitative measurement of the optical properties. The measurement at several wavelengths of these properties enables the retrieval of functional and structural parameters of the tissue. However, this requires the measurement at several spatial frequencies for each optical property. Additionally, many acquisitions are needed for each spatial frequency (at least 3). As such, a significant acquisition and processing times is needed to obtain a rapidly evolving information. As a result, this method is not suitable for a surgical application. In order to approach real-time capabilities, acquisitions are generally reduced 58 to two spatial frequencies, but this is still too slow compared to the surgeons' need. Finally, in order to overcome this problem and speed up the process, a real-time acquisition method has been developed, namely the SSOP method which aims to measure optical properties using a single image projected onto the surgical area [1]. With this technique, the number of acquisitions is reduced, as shown in By measuring these optical properties at several wavelengths, it is possible to extract the functional and structural parameters of a tissue, for this a multispectral acquisition is required [3]. All conventional multispectral acquisition solutions always involve sequentially taking images at each wavelength or requiring a complex optical system. Depending on the approach used, the scanning time can be long, the signal to noise ratio can be low, the acquisition may be affected by motion [4], and the cost of the system can be high.

Finally, if a shorter acquisition time is essential, processing time to extract optical properties, functional and structural parameters is equally important. Indeed, a shorter processing time than the desired visualization rate is of critical importance. In the SFDI/SSOP literature existing at the beginning of our thesis, no publication refers to a processing time that would allow for at least 10 images per second in order to guarantee the real-time visualization of the results.

The main objectives of our thesis following the above-mentioned limitations are:

-SSOP image artifact reduction and accuracy improvement addressed in Chapter 4 and 7.

-Design of a real time multispectral imaging acquisition system presented in Chapter 5.

-Implementation of a real-time processing and visualization imaging system shown in Chapter 6 and 7.

Artifact reduction and accuracy improvement

The artifacts on the SSOP images are due to the filtering strategy used in the Fourier domain and because it is impossible to perfectly recover an image from a single image affected by the sinusoidal pattern. Indeed, the size of the filters influences the resolution and the information that can be extracted while the shape of the filters has a strong impact on the reduction of the Gibbs effects. To improve the visual quality of the images and their accuracy, we conducted a first study presented in Chapter 4 which consisted in designing more efficient filters in Fourier domain. We investigated the effect of a set of anisotropic two-dimensional filtering windows. Both accuracy and image quality of the obtained optical properties were quantified and compared to SFDI results.

To overcome the inherent limits to this filtering technique, the use of convolutional neural network-based algorithm was used to bypass Fourier filtering set, as described in Chapter 7.

Additionally, results are compared to the multi acquisition SFDI method for validation. In Chapter 7, we also implemented a new single pattern real-time 3D profile correction to increase the accuracy of the results. In the same chapter we will demonstrate that simultaneous high-quality image reconstruction, with surface profile correction and accurate optical property (OP) extraction in real-time across large fields of view, is possible.

Real-time multispectral imaging acquisition system design

To overcome the limitations of conventional multispectral measurements, two innovative approaches have been studied: single pixel camera technology and spatiotemporal modulation.

The study on the single pixel camera, which is not presented in this manuscript, has concluded that the technique is interesting because of the high sensitivity of the detectors. However, due to the limitations of the technology in terms of acquisition time, we did not pursue with the investigations [5]. The study on spatiotemporal modulation enabled us to obtain better results. It is presented in Chapter 5 as proof of concept of the method and validated by a phantom measurement experimentation and by in vivo oxygenation measurement. More details about the real-time multispectral imaging acquisition system are also provided in Chapter 6 which deals with the realtime processing implementation.

Real-time processing and visualization imaging system

As presented above, the ability to provide real-time acquisition, processing and visualization is one of the primary goals of image-guided surgery systems. In order to add these functionalities to SSOP-based systems, we studied the massively parallel computing architectures that are Field Programmable Gate Arrays (FPGA) and Graphics Processing Unit (GPU). Indeed, they enable to

Conclusion

This chapter presented the aim of our project and summarized the limitations of the SFDI/SSOP methods we plan to use in the design of the real-time image-guided system. In addition, we provided a brief presentation of our contributions (study or method). A brief reference to the corresponding chapters was presented as well. In the next chapter, our first contribution related to improving visual image quality and accuracy of SSOP are discussed.

In the previous chapter, we have presented the major limitations to the success of our thesis project which can be summarize as the low accuracy of SSOP measurements, the poor visual image quality, the need of a real-time multispectral acquisition method, the need of a real-time processing implementation and the need of a real-time visualization. In fact, to develop an intraoperative image-guided surgery device accuracy and visual quality are mandatory. The lower accuracy and visual image quality of SSOP compared to those of diffuse reflectance imaging method SFDI is the first major limitation of our project and its use in operating rooms. In this chapter, which presents our first contribution, we will describe the first approach we studied to improve accuracy and visual quality of SSOP results, which consisted in designing a new filter mask in Fourier domain. The material for this chapter is our publication "Single snapshot of optical properties image quality improvement using anisotropic two-dimensional windows filtering" in Journal of Biomedical Optics, Volume 24, Issue 7.

4.

Context and motivations

Optical imaging methods capable of providing real-time information related to functional and structural conditions of living tissues are becoming increasingly popular [1][2][3][4][5][6]. Some medical fields in particular, such as surgery, are driving this need for real-time interpretable information content to aid decision-making in a time-constrained environment [7][8][9][10][11][12]. Several solutions have been investigated for fulfill this need, each with their pros and cons, ranging from raster scanning of microscopic information to wide-field acquisitions of diffused information [13,14]. For instance, fluorescence imaging to help surgeon visualizing structures of interest, such as lymph nodes, ureters, or micrometastasis, stands as an example where providing tissue-related information in real-time was key in clinical adoption [15][16][17][18].

Among all existing solutions, spatial frequency domain imaging (SFDI) has been recently developed with the unique capacity to provide rapidly quantitative diffuse information (absorption and reduced scattering) over large fields of view (typically greater than 10 x 10 cm²) [14,19]. Even more recently, the SFDI method has been improved to image optical properties more rapidly using 2 images [20] or even a single image [21,22]. The latter method, termed single snapshot of optical properties (SSOP), has the unique capacity to image fully in real-time with the only limitation being the reflectance signal-to-noise ratio at high camera frame rates. However, improving speed of acquisition comes at the cost of image quality, a factor that cannot be neglected when aiming at translating an imaging method to the clinic. It is therefore necessary to provide at the same time speed and good quality image so they can be fully utilized by clinicians.

In this chapter, we present means to improve image quality when using SSOPs. We investigate in particular two parameters: the filtration strategy in the spatial frequency domain, and the choice of spatial frequency. We present the effect of these parameters onto both recovered optical properties and image quality in comparison with standard, non-real-time SFDI. Finally, we propose an optimal set of parameters allowing to image optical properties accurately while maintaining good image quality.

Materials and Methods

Spatial Frequency Domain Imaging

The theory of SFDI is described in the literature [14,19]. Briefly, it consists of analyzing the spatial modulation transfer function (s-MTF) of a turbid medium for every pixel of an image at once.

SFDI makes measurements relative to a tissue-mimicking calibration phantom of known optical

properties and uses a model-based approach, typically diffusion theory or Monte Carlo simulations, to extract the optical properties. In the case of sub-surface imaging, where optical properties are considered to be independently measured at every location in the image, a fast, precomputed lookup table (LUT) can be used to recover the optical properties from only two spatial frequency images. Typically, one DC image (e.g., 0 mm -1 ), sensitive to changes in both reduced scattering and absorption, and the other AC image (e.g., 0.2 mm -1 ), mainly sensitive to changes in reduced scattering, are used [14,23,24].

The most common approach is to use a projector and a camera to acquire for each spatial frequency (e.g., f0 = 0 and fx = 0.2 mm -1 ), three images phase shifted by 120 degrees (φ1, φ2 and φ3), hence a total of six images (Fig. 4.1). When using the three-phases demodulation technique, an analytical expression can then be used to extract the AC and DC components for each pixel xi of the sample (MAC and MDC).

= ( , ) = 2 ⁄ 3 {[ ( ) -( )] + [ ( ) -( )] + [ ( ) -( )] } ⁄ (4.1)

Single Snapshot of Optical Properties

The principle of SSOP [21,22] consists of acquiring a single image at a high spatial frequency (e.g., fx = 0. or a two-dimensional (2D) Fourier transform (FT) is performed on the image, the modulation frequency is detected in Fourier domain in the same manner as in FT profilometry [25,26], and an ideal rectangular filter is used to separate low spatial frequency components from the high spatial frequency components. The MDC image is obtained using simple inverse FT and the MAC image using Hilbert transform. Following this demodulation step, data are processed in an identical fashion as with SFDI: a calibration reference is used, the medium diffuse reflectance is obtained, and finally, reduced scattering and absorption maps are extracted. There are, however, major issues of this method that impact accuracy and image quality: first, the Gibbs effect due to the ideal rectangular filtering and second, the difficulty to properly extract the AC component with just one phase.

In the following sections, we study different filtering windows to obtain the MAC and MDC images and we investigate means to reduce Gibbs effect, reduce fault value due to edge discontinuity, keep optical properties close to the reference multiphase SFDI values, and improve in the same time the visual aspect.

Design of the SSOP demodulation filters

From the literature, the classic rectangular filtering is not causal, and, even if it gives a most perfect response, is subject to ripple [27]. To reduce the ripple effect, and improve image quality, other well-known windows were investigated, such as Sine, Hann, Blackman, and Gaussian [27]. These windows have known properties in 1D (e.g., maximum side lobe level, side lobe roll-off rate, -3 dB main lobe width, and scalloping loss), but because images are two-dimensional, 2D versions of these windows have to be generated, inherently altering their characteristics. There exists different way to create 2D windows from 1D windows, such as simple product or circular rotation [28], and for simplicity of implementation in Matlab, a simple product of the x and y directions of the filter characteristic was employed.

The design rationale for 2D anisotropic demodulation filters is as follows. The filter should capture as much information as possible in the frequency domain, information that should be related to the modulation frequency of interest (here, DC or AC). The cutoff frequencies should be chosen so information content in the spatial frequency domain is well separated between modulation frequencies (DC and AC). Because there is no known method to theoretically determine the best filter shape or cutoff frequencies, we developed our own method.

In the x direction, defined as the direction that contains spatially modulated information, filters were designed around the zero frequency for the DC filter and around the projected spatial frequency for the AC filter. Filters sizes in this direction were incremented by steps of 0.0067 mm -1 to study the effect of the cutoff frequency onto demodulated images. In the y direction, because there is no spatially modulated information, the filter shape was maximized to contain as much information as possible. Filters sizes in this direction were then kept constant and equal to the maximum spatial frequency (here, 6.83 mm -1 ). Two-dimensional filters were then obtained by multiplying the filters in the x and y directions, resulting in elliptical-shape anisotropic filters.

Finally, for each 2D filter shape (with varying cutoff frequency in the x-direction), a single-phase SSOP image of a complex in vivo sample (a human hand) was demodulated. The demodulated images were compared with those obtained using standard SFDI demodulation as a reference. The best filters parameters were chosen when demodulated images with SSOP best matched the ones obtained with SFDI. These operations were repeated for different projected spatial frequencies.

Using this method, standard filters parameters were defined: center and cutoff frequencies for both DC and AC.

Effect of Spatial Frequency

As our interest is to separate DC and AC content in the frequency domain, the value of the projected spatial frequency is of paramount importance. The higher the spatial frequency, the lower the cross talk between the DC and AC components and the easier it is to separate contributions from the MDC and MAC. Within the limits of our projection system, a range of 4 spatial frequencies were tested: fx = 0.1 mm -1 , 0.2 mm -1 , 0.3 mm -1 and 0.4 mm -1 . An additional acquisition at fx = 0 mm -1 was performed for processing with SFDI. For each spatial frequency, SSOP processing was performed using the AC and DC filters and the resulting images compared both in optical properties value and in image quality.

Instrumental Setup, Phantoms and Sample

The instrumental setup was custom-built using a DMD (Vialux, Germany) for the projection of custom patterns, fiber-coupled to a 665nm laser diode (LDX Optronics, Maryville, Tennessee).

The projection system projects a sine wave pattern over a 150 × 150 mm 2 field of view at 47 cm working distance. About 1024 x1024 images are acquired using a scientific monochrome 16 bits CMOS camera (PCO AG, pco.edge 5.5, Kelheim, Germany). Polarizers (PPL05C; Moxtek, Orem, Utah), arranged in a crossed configuration, are used to minimize the contribution from specular reflections at the surface of the sample.

Silicone-based optical phantoms were built using titanium dioxide (TiO2) as a scattering agent and India ink as an absorbing agent [29,30]. One large calibration phantom was made (210 mm × 210 mm × 20 mm in size) with reduced scattering µs ' = 1.08 mm -1 and absorption µa = 0.012 mm -1 at 665 nm.

Anatomically relevant in vivo samples were used to best evaluate the potential of the proposed method for improved accuracy and image quality in clinically relevant conditions.

Data processing & results analysis

Data was processed using both SSOP and SFDI methods as described previously [21,22]. Custom processing code in Matlab (Math Works, Natick, Massachusetts) using a lookup table (LUT) extraction method was generated from White Monte Carlo and, therefore valid for high spatial frequencies beyond the diffusion limit [14].

For the evaluation of the results accuracy, optical properties maps of N = 10 hands, complex heterogeneous sample, obtained using the SSOP method for each filter set and each spatial frequency were compared with optical properties maps obtained from the SFDI method. Variations of optical properties between the two methods were calculated in terms of quality and fidelity with a custom percentage error formula over the entire image:

%= 1 × 100 × , - , , (4.6)
For the evaluation of the visual quality of the images, we measured the full width at half-maximum (FWHM) of a vein of N = 10 hands onto the optical properties maps obtained from the SSOP method using the AC and DC filters in comparison with the SFDI method. Images were also analyzed visually to take into account degradations, such as ripples and edge artefacts. Mean and standard deviations were calculated for the N = 10 samples.

Results

Filters design

We evaluated a total of 39 common filters (13 DC filters and 26 AC filters) for a total number of filtering combinations of 338. All filters equations and parameters are detailed in Appendix 1. For ease of interpretation, we present only the results obtained with the filters having the best • fx = 0.1 mm -1 : DC3 and AC2 with 17.00 % ± 1.49 % error

• fx = 0.2 mm -1 : DC3 and AC2 with 9.72 % ± 0.87 % error

• fx = 0.3 mm -1 : DC3 and AC2 with 8.84 % ± 1.07 % error

• fx = 0.4 mm -1 : DC3 and AC2 with 9.43 % ± 1.09 % error Globally, these results indicate that DC3 and DC2 filters allow to extract absorption properties more accurately than the rectangular filter (DC1). The DC3 filter consistently shows best results

Filter Name Filter Equation

Cutoff frequency in mm -1

DC1

Rectangular bandpass

( , ) = (1) × (1) for -( × ) ≤ ≤ ( × ) -( × ) ≤ ≤ ( × ) = 1.6 DC2 Blackman bandpass ( , ) = 0 -1 cos 2 n -1 + 2 cos 4 n -1 × 0 -1 cos 2 m -1 + 2 cos 4 m -1 for -( × × 3) ≤ ≤ ( × × 3) -( × ) ≤ ≤ ( × ) = 0.16; 0 = 1- 2 ; 1 = 1 2 ; 2 = 2 = 3 DC3 Sine bandpass ( , ) = N -1 × m -1 for -( × × 1.88) ≤ ≤ ( × × 1.88) -( × ) ≤ ≤ ( × ) = 2.25

AC1

Rectangular highpass ease of interpretation (scale on the right).

( , ) = (1) × (1) for (( × ) -( × × 0.32)) ≤ ≤ ( × ) -( × ) ≤ ≤ ( × ) = 1.3 AC2 Blackman bandpass ( , ) = 0 -1cos 2 n -1 + 2cos 4 n -1 × 0 -1 cos 2 m -1 + 2cos 4 m -1 for (( × ) -( × × 1)) ≤ ≤ (( × ) + ( × × 1)) -( × ) ≤ ≤ ( × ) = 0.16; 0 = 1- 2 ; 1 = 1 2 ; 2 = 2 = 1.3 AC3 Blackman highpass ( , ) = 0 -1cos 2 n -1 + 2cos 4 n -1 × 0 -1cos 2 m -1 + 2 cos 4 m -1 for (( × ) -( × × 1)) ≤ < ( × ) -( × ) ≤ ≤ ( × ) (1) × (1) for ( × )≤ ≤( × ) -( × ) ≤ ≤ ( × ) = 0.16; 0 = 1- 2 ; 1 = 1

Reduced scattering images

Results quantifying accuracy for the selected filters on N = 10 hands are shown in Fig. 4.5., colorcoded for ease of interpretation. Overall, the mean error expressed in percentage value for reduced scattering varies from 14.47 % in the worst case to 7.46 % in the best case. Results for all 338 filters combinations are shown in Appendix 1. More precisely, the best combinations to extract reduced scattering, as a function of spatial frequency, are as follows:

• fx = 0.1 mm -1 : DC2 and AC3 with 12.33 % ± 1.37 % error

• fx = 0.2 mm -1 : DC2 and AC3 with 8.10 % ± 1.12 % error

• fx = 0.3 mm -1 : DC2 and AC2 with 7.46 % ± 1.02 % error

• fx = 0.4 mm -1 : DC2 and AC2 with 8.41 % ± 1.08 % error spatial frequency is increased and ripples around the edges as well as on the background are strongly diminished as filtration quality is improved.

Finally, a movie of a hand was acquired making a quarter turn with tight fingers at the beginning and spread fingers at the end to demonstrate the filtering improvement onto data acquired in realtime (Fig. 4.9.). This movie is acquired at a spatial frequency of 0.4 mm -1 , a true frame rate of 50 images per second and processed with the standard rectangular DC1 and AC1 filter combination

Discussion

In this study, we have investigated different 2D filtering strategies allowing improvement of the SSOP method's visual aspect while preserving good optical properties measurement accuracy.

Additionally, this work provided an opportunity to precisely design equations for 2D anisotropic filters in the spatial frequency domain. First, we were able to define DC and AC filters shape and cutoff frequencies by comparing demodulation accuracy with SFDI as a reference. Then, the filters that give best results in terms of optical properties and image quality were combined, applied onto complex in vivo objects (N = 10 hands), and results presented in this chapter. The extracted optical properties at different spatial frequencies from these filters were compared to SFDI, and their visual aspects were assessed with measurement at mid-height of the veins.

In summary, using an acquisition spatial frequency of 0.4 mm -1 , combinations using the DC1 and good quantitative value while providing an improved visual aspect at the same time. In addition, not quantified in the visual quality, ripples are significantly reduced from both the background phantom and the in vivo sample.

However, this work is not without limitations, in particular regarding the quantitative comparison metrics, the precision of the comparison reference using SFDI, and the image quality metrics.

Regarding the quantitative value metric, the formula used (Eq. 4.6) gives a unique value for all the pixels of the image. This implies that the pixels at the background and at the edges of the hand have the same importance than the pixels at the center of the hand. For example, in Fig. 4.7.

(absorption images), at 0.4 mm -1 , vein structures of the hand using SSOP std are better visualized than using SSOP opt2; however, the mean error percentage of SSOP opt2 is lower since there are more errors on SSOP std edges and background. Therefore, when choosing a filtering combination, a compromise is necessary between the defaults at the edges of the image, the defaults in the background and the accuracy of the values in the region of interest, depending of the imaging application.

Regarding the use of SFDI as a reference for the quantitative comparison, noise related to multiplephases acquisition can result in stripe artefacts, as seen in reduced scattering in It should also be noted that the percent error was computed over the entire hand giving a result that appears to be high. We purposely chose this approach to evaluate errors on both flat zones (back of the hand, for instance) and zones with height changes (edge of the fingers, for instance). The errors resulting from the edge artifacts that we attempted to reduce in this work contribute largely to a high error value for the entire hand. Clinically, the accuracy necessary for proper use of this technology remains unknown and dependent on the clinical scenario.

We chose to perform the image quality assessment by measuring the width of a vein on the back of N = 10 hands. This choice was motivated by the fact that it reflects the image quality metric expected by the end-user of the method. We did investigate other known image quality metrics, such as BRISQUE, NIQE, or SSIM. However, none of these metrics provided satisfactory results

for quantifying image quality. This is partially due to the fact that these algorithms have been developed for a specific context of image quality that does contain features in the entire image. In our case, in vivo images contain both a flat, featureless zone (the reference phantom) and a zone containing features (the hand). This results in an inaccurate assessment of the impact of image artifacts, such as ripples, on flat zones, leading to a disagreement with simple visual evaluation.

We, therefore, preferred to use a simpler metric reflecting the resolution of absorption within the image, by measuring the width a vein.

Overall, this quantitative assessment proves to be effective for absorption images, showing an image quality improvement as spatial frequency increases, as one would expect and subjectively confirms. However, this measurement does not take into account the background of the image or the edge ripples. Finding a way to better quantify image quality would certainly improve such study. In addition, because veins at the back of hands are common absorption features, we were able to perform this measurement on several hands. Unfortunately, scattering features that would be similar between individuals are not common. We, therefore, did not perform the image quality assessment onto scattering images. However, this is not a strong limitation since it is well known that scattering images do not suffer quality degradation as much as absorption images. This is due to the fact that higher spatial frequencies are preserved in scattering images and, as a consequence, image quality is higher.

Of interest as well, this study offers a better understanding of both quantification and image quality as a function of spatial frequency. From a pure image analysis standpoint, it may appear obvious that increasing spatial frequency would improve demodulation in the case of SSOP by providing better separation between DC and AC content. However, image quality suffers when increasing spatial frequency with an increase in noise due to lower signal-to-noise ratio that, depending on the instrumentation used, can lead to inaccuracies in extracting optical properties [14]. This topic is currently being investigated in particular regarding instrumental contributions to both bias and image quality. In addition, it is well known that increasing the spatial frequency results in more superficial measurement in the spatial frequency domain [14]. This is well demonstrated in the reduced scattering maps (Fig. 4.8.) with a noticeable change in properties with spatial frequency.

On the contrary, absorption maps remain nearly identical for all spatial frequencies (Fig.

4.7.).

This result is interesting in the context of surgical guidance where functional parameters, derived from absorption, are particularly useful. This result justifies the use of a high spatial frequency acquisition in the interest of image quality while preserving accurate functional measurements.

More studies are necessary to understand the implication of this result onto measuring reduced scattering in the context in vivo applications.

Finally, it is interesting to note that some filter combinations perform better for absorption accuracy while others perform better for reduced scattering accuracy. One may, therefore, consider the option of using different filter combinations for extracting absorption and reduced scattering to obtain as accurate results as possible. This, however, comes at the cost of an increase in processing time.

Conclusion

In our goal to develop a new intraoperative image-guided surgery device based on the diffuse reflectance imaging method SFDI, we decided to combine the fast acquisition SSOP method of SFDI with temporal modulation of light to discriminate optical properties of tissues acquired at several wavelengths simultaneously. Legacy SSOP methods which uses rectangular filter windows at the filtering stage representing a major limitation for image quality. This filtering strategy is subject to Gibbs effect and generate unwanted artifacts essentially localized on the edges and that lead to reduced overall accuracy and to poor visual rendering of the images. To overcome this limitation, in this chapter new filtering strategies for SSOP were presented. The presented method involves the use of 2D anisotropic windows for filtering in the Fourier domain. We evaluated the combination of 39 common filters (13 DC filters and 26 AC filters) based on 2D windows for a total number of filtering combinations of 338 and they were finally tested onto N = 10 hands.

Using custom mean error percentage computation and FWHM vein measure, this study demonstrates that it is possible to significantly improve image quality while preserving accuracy in optical properties measurements using Blackman and Sine windows at high spatial frequencies (typically, 0.4 mm -1 ). Overall, using optimized filters, mean errors in predicting optical properties using SSOP remain under 8.8% in absorption and 7.5% in reduced scattering.

The next chapter will be dedicated to address the limitations of conventional multispectral acquisition method which involve sequential acquisition at each wavelength (or require a complex optical system) and their consequences: a long acquisition time, a low signal to noise ratio, an acquisition affected by motion, and the cost of the system. We will present the work we have conducted to develop a real-time multispectral acquisition system which use temporal modulation of light. Using this method, a proof of concept of a quantitative oxygenation measurement system using spatiotemporal modulation of light will be presented.

In the previous chapter we proposed an approach to solve the first limitation of our project which is the accuracy and visual quality of the optical properties measured with SSOP by designing new 2D anisotropic filters based on 2D windows such as Sine and Blackman. The measurement of optical properties at several wavelengths will enable one to provide essential physiological and structural information to surgeons such as oxygenation, lipids, melanin and water levels, as well as scattering power and scattering amplitude. However, as explained in Chapter 3, the multispectral acquisition technique used influences the quality of the measurements and the proportion of the designed system to be real-time. To overcome these limitations, we propose a real-time multispectral acquisition technique based on temporal modulation of light. Furthermore, in this chapter we will tackle the adaptation of the SSOP method to build a quantitative oxygenation imaging system using spatiotemporal modulation of light and its validation on phantoms and on in vivo experiments. The material for this chapter is our publication "Real-time,

wide-field, and quantitative oxygenation imaging using spatiotemporal modulation of light" in

Journal of Biomedical Optics, Volume 24, Issue 7.

5.

Context and motivations

Diffuse optical imaging using either temporal or spatial modulation of light has shown strong promise for quantitatively assessing functional parameters in living tissue, such as content in oxyhemoglobin, deoxy-hemoglobin, water and lipids. Such parameters can be used to determine the physiological condition of the tissue and provide critical information to healthcare professionals for making decision. For instance, temporal frequency domain monitoring of breast tissue during neoadjuvant chemotherapy treatment can help differentiating responding versus non-responding patients within 10 days of treatment [1,2]. In other studies, spatial frequency domain imaging (SFDI) has demonstrated the potential to diagnose patients with diabetic foot ulcers early, providing means to adapt treatment towards limiting morbidity and decreasing healthcare burden [3,4].

While these methods have demonstrated a strong potential to aid in decision-making and solve concrete clinical problems, their adaptation to challenging clinical environments has been very limited. Surgery is one of these environments where information must be provided quantitatively (i.e., repeatable & interpretable), in real-time (i.e., > 25 images per seconds), and over a large field-of-view (i.e., > 10 x 10 cm²) [5]. To attain this challenging objective, our group has been developing solutions based on SFDI. Because SFDI works in the frequency domain, it allows the capture of information for all the pixels of an image at once, providing a significant advantage over raster-scanning methods. The original developments of SFDI relied on the acquisition of several images at single wavelengths, providing quantitative multipixel optical properties measurements [6][7][8]. More recent developments have introduced near-real-time (e.g., 10 frames per second (fps)) and real-time (i.e., > 25 fps) imaging methods such as Single Snapshot of Optical Properties (SSOP) [9][10][11]. These methods process the acquired information directly in the frequency domain offering the opportunity to drastically reduce the acquisition time down to a single acquired image for obtaining optical properties maps of living tissues. However, these methods still acquire only one wavelength at a time, making them sub-optimal for imaging living tissue parameters (e.g., chromophores concentration) in real-time.

In order to deploy such technology in the operating room, methods capable of performing multior hyper-spectral imaging quantitatively, over a large field-of-view and in real-time must be developed. Recent work has proposed dual-wavelengths imaging for extracting blood parameters such as the concentration of oxyhemoglobin, deoxyhemoglobin or oxygenation. For instance, SSOP was used in conjunction with speckle imaging using a scientific CMOS camera by sequentially acquiring wavelengths, leading to 16 fps imaging [12]. However, this work uses a mechanical optical chopper and low power (150mW) laser diodes, making it difficult to achieve higher frame rates and limiting its scalability to a higher number of wavelengths. Another more recent work used square patterns to improve switching time using a standard digital micromirror device (DMD) to project 3 wavelengths sequentially at 3 phases (standard SFDI acquisition), and compressed single pixel imaging leading to 12.5 fps imaging [13]. However, relying on sequential frames to acquire optical properties maps naturally restricts the possibility to achieve higher frame rates. As of today, because of these instrumental limitations, no method has been developed for real-time (i.e., > 25 images per seconds), wide-field (i.e., > 10 x 10 cm²) and quantitative imaging for clinical applications.

In this work, we present a method using both spatial and temporal modulation of light to provide images interpretable by clinicians in time-constrained environments such as surgery. We propose here a method where sources are modulated sinusoidally in time using signal generators (no mechanical chopper) and projected in the form of single-frame patterns for processing using the SSOP method (no need for sequential frames). Sources are modulated in time to encode wavelengths, and combined in a multifiber bundle connected to a spatial-light modulator. Light is then projected with specific patterns onto the field of view and images captured using a fast camera. The obtained 3D cube of data (2D + time) is then processed in the temporal frequency domain to separate wavelengths, and the resulting images are processed in the spatial frequency domain using the SSOP methodology to extract optical properties for each wavelength and then calculate blood parameters (oxyhemoglobin, deoxyhemoglobin and oxygenation). In the first part of this chapter, we present the imaging method, processing and hardware. Then we validate the method in comparison with standard, non-real-time SFDI onto tissue-mimicking phantoms.

Finally, we perform a validation onto a living hand during an arm cuff occlusion. Altogether, this study contributes to the translation of endogenous diffuse optical imaging for surgical guidance.

Material and methods

Spatio-temporal modulation of light

The concept of multispectral imaging using temporal modulation of light sources for wavelength encoding is not new, nor are multispectral spatial frequency domain methods [13][14][15][16][17][18]. However, until now, imaging living tissue parameters over large fields-of-view at high resolution, quantitatively and in real-time remained a significant challenge. The principle underlying the spatio-temporal modulation of light method is shown in Fig. 5.1.. Discrete sources are modulated in time, each wavelength being at a different specifically chosen temporal frequency. The light outputs from all sources are combined using a multi-fiber bundle connected to a spatial light modulator. Patterns of light combining all temporally modulated wavelengths are projected onto the field-of-view and a fast camera is used to acquire the scene. The acquired images are first processed using a temporal Discrete Fourier Transform (DFT) to accurately separate wavelengths contributions. Finally, the resulting images are processed in the spatial frequency domain using Single Snapshot of Optical Properties (SSOP) methodology to obtain the optical properties at the chosen wavelengths.

One of the significant hurdles in implementing such method is in ensuring that wavelengths are precisely and correctly modulated, sampled and processed to unambiguously extract the information content contained in each temporally modulated source. Not performing these operations correctly would result in inaccuracies in the extracted optical properties and therefore in the assessment of tissue parameters. Two strategies can be used to perform these measurements: either a fine temporal sampling of the scene allowing to use a generic Fast Fourier Transform In this work, because we focus on being real-time, we chose to implement a custom DFT to allow ultrafast computation of the chosen sources' contributions. The discrete Fourier transform of one pixel in time is given by:

( ) = ( ). ( ) (5.1)
and the proper frequencies that should be used to modulate are exact frequencies in the discrete Fourier domain, for instance, in the case of four sources:

N F k f N F k f N F k f N F k f e e
= = = = (5.2)
with Fe being the sampling frequency, N the number of acquired images, and ki determined such that the captured signal spans an integer number of periods [19,20].

If these conditions are gathered, on top of obvious considerations such as avoiding harmonics or standard electrical current frequencies, and along with a precise synchronization of the camera and signal generators, the contribution of each source can be precisely and unambiguously determined.

Acquisition and Processing

In this work, we performed acquisitions at a sampling frequency Fe = 100 Hz, processing N = 10 images, and modulated two wavelengths to determine oxygenation, λ1 = 665 nm with k1 = 3 (f1 = 30 Hz) and λ2 = 860 nm with k2 = 2 (f2 = 20 Hz) [21]. The number of images (N=10) was chosen such that all parameters have been adjusted to satisfy the conditions for proper DFT processing (Sec. 2.1) while maintaining a proper signal-to-noise ratio (SNR) in the acquired images. We implemented a custom rolling window algorithm of the DFT method onto Matlab with the parallel computing toolbox (using a graphics processing unit, GPU) by adding the N+1 term of the DFT sum and subtracting the first term of the DFT sum in real-time. More details regarding the GPU implementation of the processing method are available elsewhere [20].

All acquisitions have been performed on the imaging system described in the next section. In order to handle the large flux of data, as well as control the hardware and perform GPU processing of the acquired images, a personal computer with the following characteristics was used: Intel i7-7800x 3.5 GHz central processing unit, 16 GB of RAM, four 1 TB solid state drives for data acquisition and one 500 GB solid state drive for system operation, and an NVIDIA GeForce GTX1080TI GPU.

Imaging system

The instrumental setup was custom built using a DMD (Vialux, Germany) for the projection of custom patterns, fiber-coupled to two laser diodes at 665nm and 860nm (LDX Optronics, Maryville, TN). Diodes were mounted in a temperature-controlled mount (TCLDM9, Thorlabs, Newton, NJ). Intensity from the diodes was controlled using current controllers (TDC240C, Thorlabs, Newton, NJ) and temperature controlled using a thermoelectric cooler controllers (TED200C, Thorlabs, Newton, NJ). Diodes were modulated using a dedicated 20 MHz dual channel function generator (4047B, BK Precision, Yorba Linda, CA). The laser diodes were coupled into a 2-to-1 fiber bundle (BFY1000HS02, Thorlabs, Newton, NJ) connected to the DMD.

The projection system projects a sine wave pattern over a 200 × 150 mm 2 surface at 45 cm working distance. Images were acquired over a 130 x 130 mm² field of view and at a resolution of 1024 x 1024 pixels using a scientific-grade monochrome sCMOS camera (PCO Edge 5.5, Kelheim, Germany). The camera and signal generators were synchronized together using a second 20 MHz dual channel function generator as a reference (4047B, BK Precision, Yorba Linda, CA).

Polarizers (PPL05C, Moxtek, Orem, UT), arranged in a crossed configuration, were used to minimize the contribution from specular reflections at the surface of the sample. A silicone-based optical phantom was used for calibration and built using titanium dioxide (TiO2) as a scattering agent and India ink as an absorbing agent [22]. Its large size (210 mm × 210 mm × 20 mm) accommodates the system's field-of-view with reduced scattering µs' = 1.08 mm -1 and absorption µa = 0.012 mm -1 at 665 nm, and µs' = 0.78 mm -1 and µa = 0.022 mm -1 at 860 nm.

Validation experiments

5.2.4.1.

Acquisition protocol SFDI and spatio-temporal modulation of light experiments were performed sequentially. SFDI was used as a reference by acquiring a total of six images for each wavelength (with no temporal modulation) at spatial frequencies of 0 mm -1 and 0.3 mm -1 with three evenly spaced phases each to extract the amplitude modulation and the calibrated diffuse reflectance. A lookup table (LUT) generated from Monte-Carlo was used to extract the optical properties for each wavelength following acquisition. Extensive details regarding the SFDI acquisition and LUT processing method can be found elsewhere [7,23,24]. Spatio-temporal acquisitions were performed at 30 Hz for 665 nm and 20 Hz for 860 nm. Acquisition rate of the camera was set at 100 fps. The SSOP projection frequency was 0.3 mm -1 and the acquired data post-processed, first in time to extract contributions from the sources, and then in space to extract the optical properties for each wavelength. More details regarding the temporal processing can be found here in Ref. [20] and details regarding the SSOP acquisition and processing method can be found in Refs. [9,11].

Phantoms experiments

A set of 10 tissue mimicking phantoms were built using titanium dioxide (TiO2) as a scattering agent and Nigrosin as an absorbing agent [25]. Phantoms were all acquired using both SFDI and spatio-temporal modulation of light, post-processed, and the mean error in percent calculated over 30 x 30 mm² region of interest taken at the center of each phantom.

In vivo experiments

A human hand was used to illustrate the capacity of the method to work properly during in vivo conditions and in real-time. First, a still hand was imaged using SFDI and spatio-temporal modulation of light to compare in vivo results using both methods. Mean error in percent calculated over 50 x 50 mm² region of interest taken on the back of the hand. Then, a simple arm cuff occlusion experiment was performed during 10 min, with 3 min of rest, 4 min of occlusion and 3 min of recovery. Images were acquired and post-processed using the spatio-temporal modulation of light method and the absorption coefficients used to determine oxygen saturation through the calculation of the concentrations of oxy-and deoxy-hemoglobin [21].

Results

Phantom experiments

Results from the phantom experiments are summarized in Fig. 5.2.. On average, optical properties obtained using spatio-temporal modulation of light show good agreement with the standard SFDI method. Mean errors were calculated for all the pixels in a 30 x 30 mm² region taken in the middle of each phantom. For 665 nm, the mean error is 1.8 % ± 1.2 % in absorption and 3.3 % ± 1 % in reduced scattering. For 860 nm, the mean error is 4.2 % ± 2.1 % in absorption and 2.

In vivo experiments

Images from the still hand are shown in Fig. 5.4. for both wavelengths, in absorption and reduced scattering. Results from standard SFDI are shown, as well as spatio-temporal modulation of light, a percentage error map comparing both methods and a line profile. Overall good visual agreement is confirmed in these images. Mean errors were calculated for all the pixels in a 50 x 50 mm² region on the back of the hand. For 665 nm, the mean error is 1.9 % ± 2.0 % in absorption and 2.0 % ± 0.7 % in reduced scattering. For 860 nm, the mean error is 3.4 % ± 1.9 % in absorption and 0.3 % ± 1.7 % in reduced scattering. These results confirm the phantom experiments with similar error values on average and allow to appreciate the image quality of the extracted optical properties maps. In particular, the vein structures are clearly visible on the absorption images confirming the potential for the spatio-temporal modulation of light method to perform adequately for clinically relevant use.

Discussion

In this work, we present the application of a novel method, spatio-temporal modulation of light, to perform real-time (100 fps), wide-field (13 x 13 cm²) and quantitative imaging of blood parameters in vivo. This method relies on the modulation of light sources in time to encode wavelengths and on the modulation of the projection of these sources in space to extract the optical properties of living tissues. Combined together, these spatio-temporal modulations allow to gather multispectral quantitative optical information in real-time (here at 100 fps and post-processed). Altogether the method enables to provide healthcare professionals with interpretable information that can be used in real-time to guide decision making in challenging clinical situations, such as surgery.

Compared to previous work, this study augments the capacity for performing multispectral diffuse optical imaging quantitatively at a high frame rate and over a large field-of-view. For instance, sequential acquisition of two wavelengths using both SSOP and Laser Speckle Imaging acquisition methods led to 16 fps imaging of simultaneous blood parameters and perfusion imaging over a wide field of view (15 x 15 cm²) [12]. In more recent work, compressed sensing was used along with square patterns projections and phase shifting leading to 12.5 fps imaging of blood parameters over a smaller field of view (3.5 x 3.5 cm²) [13].

However, this work comes with limitations. First, in this implementation, the SSOP method results in slight image degradation [9][10][11] and the 3D profile of the sample is not measured, which can cause errors in extracting optical properties onto complex, non-flat objects (~10% error per cm for absorption) [26]. Second, noise in oxygenation data, mainly caused by the poor quantum efficiency of the camera at 860 nm, as well as slight synchronization errors between the camera and the sources, resulted in the impossibility to observe fine variations such as the pulsatile nature of blood flow. Improvements in the imaging system will include stronger illumination at 860 nm and a custom light source for optimal camera-source synchronization to allow the observation of finer variations of oxygenation. Third, currently images are acquired in real-time but post-processed following acquisition. We are currently working on a hardware and software solutions enabling both real-time acquisition and real-time processing, making this technology deployable for realtime visualization of the surgical field [24,27,28]. Finally, while the concept can be expanded to virtually many more wavelengths, there is a significant instrumental challenges in acquiring properly all wavelengths with a good enough SNR, and within the limitations of the spectral quantum efficiency as well as the dynamic range of the camera. For instance, in this work, the 860 nm optical properties, necessary to optimally separate oxy-and deoxyhemoglobin, shows significantly more noise compared with 665 nm due to the decrease of the spectral quantum efficiency of silicon-based camera technology. While manageable for a few wavelengths, this certainly limits the acquisition of a large number of wavelengths concurrently.

One advantage of this method though is its capacity to work in real-time. Current CMOS technology offers cameras at a very reasonable price permitting good resolution and dynamic range imaging up to 500 fps. Our GP-GPU implementation of the temporal DFT in CUDA permits processing down to 18µs per wavelength at a resolution of 512 x 512 pixels [20]. However, the rolling window processing approach comes at the cost of a short delay between the acquired scene and the displayed information (here 1/10 th of a second, i.e., 100ms). One potential solution to decrease this delay is to increase the acquisition rate (500 fps for instance) but it requires to increase the light fluence from the sources to maintain proper SNR and prevent noise in the images.

Finally, current state-of-the-art processing methods in the spatial frequency domain include optimized look-up tables (down to 18ms for 572 x 672 pixels) [24] as well as deep learning & machine learning methods (down to 200ms for 696 x 520 pixels) [28,29], and soon, faster methods that will enable real-time processing and visualization of data will certainly be developed. Efforts are ongoing to combine rapid processing of SSOP data, with spatio-temporal modulation of light onto a clinically compatible device to test and validate the capacity for diffuse optical imaging to provide relevant data in real-time during concrete clinical applications.

Conclusion

For the design of an image-guided surgery system, we have first worked to improve the accuracy and the visual quality of the optical properties images by proposing new filters more efficient for SSOP method. In this chapter, we introduced a new multispectral acquisition method with the aim of being real-time, limiting motion-related artifacts and improving the signal-to-noise ratio. This technique used temporal modulation of light to encode wavelength and Discrete Fourier Transform at demodulation stage to recover each wavelength contribution. It was then used to design a quantitative oxygenation level measurement system using spatiotemporal modulation of light, meaning optical properties encoded in space and wavelength encoded in time. The system uses 101 two wavelengths 665nm and 860nm allowing an accurate measurement of oxygenation level. It was validated on phantom and in vivo (i.e. on hand) against the results obtained by SFDI method where a sequential multispectral acquisition was performed. The results show an error rate of less than 5% for the estimation of optical properties. The system was then used during an arm cuff occlusion experiment demonstrating the ability to detect ischemia in relevant condition such as surgery. The study ends by pointing out that, although being a real-time acquisition method it requires a real-time implementation to guarantee the real-time constraint of 25 frames per second we target for the system. In the following chapter, the real-time implementation of the oxygenation imaging system using GPU computation will be presented and validate again a simpler MATLAB implementation. We will provide more details on the equations used and how we implemented them on MATLAB and GPU in C CUDA. Real-time processing and visualization are key for the design of a surgical guidance system. This is why the filtering and multispectral acquisition techniques we have presented in the two previous chapter have been carried out taking into account the time constraint. The algorithms were designed to be simple and to require few computation operations, however this is not a guarantee to achieve real-time processing and visualization. We therefore chose to use GPUs as a computational power to reach our goal at a lower cost and this will be the content of this chapter.

Indeed, the acceleration of computational operations often requires supercomputers. However, since this solution is not always within everyone's reach, the alternative found is to use GPUs that are less expensive and take advantage of their computational power. This new way of using GPUs, technically referred to as GPGPU (General-purpose Processing on Graphics Processing Units), makes it possible to use the massively parallel architecture of GPUs particularly useful for operations on matrices, not only for its graphical functions, but also to perform calculations that would be too slow on CPU. The website of the graphic card manufacturer NVIDIA (which coincides with the brand of our choice) gives a more than exhaustive list of examples of areas of application where this is beneficial. These include Bioinformatics, Digital Fluid Mechanics, Electronic Design Automation, Computer Imaging and Vision, Medical Imaging, Digital Analysis, and weather and climate modelling to name a few. More details on the architecture and how to program these GPUs are given in the appendix. The material for this chapter is our publication "Real-time optical properties and oxygenation imaging using custom parallel processing in the spatial frequency domain" in Biomedical Optics Express, Volume 10, Issue 8.

6.

Context and motivations

The use of optical measurement methods have shown a strong potential for improving patient treatment planning, surgical guidance and post-operative monitoring [1][2][3]. Among the many methods developed, diffuse optical imaging methods have already demonstrated capabilities for the quantitative measurement of the functional and structural state of living tissues over a wide field of view [4][5][6][7]. One of these diffuse optical imaging methods termed Spatial Frequency Domain Imaging (SFDI) uses spatial modulation of light (i.e. patterns of light), and from a few images acquired (typically 6 images) allows the determination of wavelength-dependent optical properties of the entire surface of a turbid sample such as living tissue [8][9][10]. By repeating these acquisitions at several wavelengths, the method allows to measure functional parameters such as oxy-hemoglobin, deoxy-hemoglobin, water and lipids and structural parameters of the tissue [11,12]. The need to acquire several images at multiple wavelengths has been for a long time a main limitation for the use of SFDI technology in time-sensitive environment such as surgery [13][14][15].

To overcome this limitation, we first introduced Single Snapshot of Optical Properties (SSOP), a specific SFDI method allowing to extract optical properties from a single acquired image in realtime at the cost of image degradation [16,17]. More recent work optimized SSOP acquisition and processing to improve imaging quality and significantly reduced the amount of artifacts from a single phase measurement [18]. However, while only a single image was necessary for extracting optical properties at a single wavelength, multispectral acquisition were performed by using a wheel filter, a tunable filter, or a multispectral camera with sequential acquisitions to distinguish each wavelength [19][20][21]. To overcome this limitation, temporal encoding of wavelengths combined with the spatial encoding for optical properties has allowed true real-time acquisition of multispectral sets of optical properties [22][23][24].

The next logical step in providing real-time tissue status assessment is to provide real-time processing capabilities. Many contributions have been made on accelerating processing methods for SFDI, especially on interpolation techniques used to extract optical properties from linearized lookup tables in diffuse reflectance or by using machine or deep learning techniques [25][26][27].

While the state-of-the-art method allows to process large sets of data (572 x 672 pixels) in 18 milliseconds using linearized lookup tables in diffuse reflectance, this process is only applied to a single wavelength, and does not include temporal demodulation of wavelength information nor the extraction of functional parameters.

The objective of this study is to explore the possibilities offered by graphic processing unit (GPU) processing of SFDI data modulated in time (for wavelength encoding) and in space (for optical properties extraction) for fast quantitative multispectral diffuse optical imaging. More specifically, we proposed GPU implementations to perform real-time processing of spatio-temporally modulated datasets, allowing the extraction of optical properties maps and functional parameters as fast as 1.6 milliseconds for a 1 megapixel dataset. We take advantage of the parallel calculation capabilities of general purpose GPU to improve processing time. In the first part, we briefly describe the method for spatio-temporal modulation of light for the extraction of optical properties at two wavelengths and computation of tissue oxygen saturation. In the second section, we describe the specific processing implementations, both in CPU using MATLAB, in GPU using MATLAB and in GPU using a custom developed CUDA (Compute Unified Device Architecture) algorithm.

Finally, we test these implementations onto real SSOP data obtained from a human hand and compare their different processing times and accuracies.

Materials and methods

Principle of spatio-temporal modulation of light

A detailed description and validation of the spatio-temporal modulation of light method used in this study is available here [24]. Briefly, to extract the optical properties of the tissue of interest, patterns of light are projected onto the tissue to obtain optical properties (as with any SFDI acquisition), but the light source is composed of discrete wavelengths that are encoded by modulating the sources in time (temporal modulation). To perform this acquisition, as shown in Fig. 6.1., the typical system is composed of several modulated laser sources, a projection system and a high speed camera. From the acquired images, the spectral contribution of each modulated wavelength is extracted by means of discrete Fourier transforms carried out in time. The resulting images are then processed in the spatial frequency domain to quantify the optical properties of the measured tissue using SSOP. Finally, the optical properties measured at different wavelengths are used to determine the physiological and structural parameters of the studied tissue.

Wavelengths modulation and demodulation

The wavelengths multiplexing step consists of modulating in time i laser diode sources with different wavelengths λi in a sinusoidal manner at different frequencies | Chapter 6. Real-time optical properties and oxygenation imaging using GPGPU 109

In the wavelengths demodulation step, the spectral contribution of each wavelength is isolated from the acquired stack of N images by processing in time the amplitude modulation Xk at frequency Fi for each pixel pxy. This operation is performed continuously in time using discrete Fourier transforms (DFT) in a rolling window scan strategy [24,28]. At the end of the demodulation process, this scan strategy enables an output framerate equal to that of the acquisition stage framerate. Xk is given by: Finally, because each pixel pxy is independent in time from its neighbors, the method is perfectly adapted for GPGPU computing for temporal de-multiplexing.

( ) ( ) 11

Optical properties extraction

SSOP is a particular implementaiton of SFDI that is used for fast extraction of optical properties [16,17]. From a single acquired frame per wavelength, by using forward and inverse 2D fast Fourier transform (FFT), filtering operations in the spatial Fourier domain are performed with custom-made filter masks DC (Low pass) and AC (High pass) [18], allowing to extract the amplitude modulation MDC and MAC at respectively the spatial frequency 0 mm -1 and at the projected spatial frequency, e.g. fx = 0.2 mm -1 . From these measurements, the diffuse reflectance images Rd,DC and Rd,AC are obtained respectively at the spatial frequency 0 mm -1 and fx > 0 by using a calibration measurement:

( ) ( ) ( ) ( )

Physiological and structural parameters extraction

In order to determine the physiological parameters of the sample, it is necessary to decouple the contributions of the intrinsic biological chromophores. The optical absorption at a single wavelength i l is given by:

() () 1 M a i m im m c ml el = = å (6.7)
Where ( ) ai ml is the measured optical absorption coefficient, ( ) mi el the chromophore extinction coefficient at i l , m c is the chromophore concentration, and M the number of chromophores [30].

Considering the equation (6.7), the measurements of the absorption coefficient at different wavelengths i l can be expressed as a system of i equations: For example, solving a system of 4 equations (4 wavelengths measurements) leads to the determination of 4 chromophore concentrations such as oxy-hemoglobin ( 2HbO ), deoxy- hemoglobin ( HHb ), water, and lipids.

() () () () 11
In our case, we used only 2 dedicated wavelengths (665 nm and 860 nm) allowing the extraction of blood oxygenation with high accuracy, assuming a 50% water-volume fraction and a 2% lipidvolume fraction [11]. Saturated oxygenation ( 2StO ) level can be computed using the following expression:

() 2 2 2
% 100 HbO StO HbO HHb =´+ (6.9)

Further optical parameters such as scattering amplitude A and scattering power b can also be determined by fitting the reduced scattering coefficient measurements over several wavelengths [30].

Imaging system

The imaging system was custom built using a digital micro-mirror device (DMD, Vialux, Germany) for the projection of custom patterns, fiber-coupled to two laser diodes at 665 nm and 860 nm (LDX Optronics, Maryville, TN). Intensity from the diodes are modulated over time at frequency 20 Hz (860 nm) and 30 Hz (665 nm) corresponding to k = 2 and k = 3, respectively, for a chosen number of images (N = 10) and a frame rate of 100 Hz. The projection system projects a sine wave pattern at 45 cm working distance. Images were acquired over a 130 x 130 mm² field of view and at a resolution of 1024 x 1024 pixels using a scientific CMOS camera (PCO Edge 5.5, Kelheim, Germany). The camera uses a high speed frame grabber (microEnable 5 marathon AF2, SiliconSoftware, Mannheim, Germany) and was synchronized with the modulation signal generators. Polarizers (PPL05C, Moxtek, Orem, UT), arranged in a crossed configuration, were used to minimize the contribution from specular reflections at the surface of the sample. A siliconebased optical phantom (210 mm × 210 mm × 20 mm) was used for calibration and built using titanium dioxide (TiO2) as a scattering agent and India ink as an absorbing agent [31] with (

Processing

General processing workflow

The processing diagram flow of the performed implementations is illustrated in Fig. 6 element-wise product between the filter and the data, and inverse FFT. Filters used for demodulation have been described previously [18]. All the previous steps are done for both the calibration and the tissue measurement. Then, diffuse reflectance values Rd,DC and Rd,AC at wavelengths 665 nm and 860 nm are calculated with Eq. (6.5) and (6.6). Finally, the optical properties a m and ' s m are estimated with various processing methods. Two types of lookup tables (LUTs) were used to extract optical properties, one using linear optical properties LUT (standard method, used for reference) and the second using linear diffuse reflectance LUT, the fastest processing method described so far in the literature [25]. In order to study the compromise between computation time and accuracy, the latter LUT (linear diffuse reflectance) was further divided into a small-size low density table (with steps of 0.001 in Rd) and a large-size high density table (with steps of 0.0001 in Rd). Both LUTs have the same range.

In the third and last task, oxy-hemoglobin ( 2HbO ) and deoxy-hemoglobin ( HHb ) are extracted by solving a system of 2 equations (Eq. (6.8)) from the measured absorption coefficients and | Chapter 6. Real-time optical properties and oxygenation imaging using GPGPU 114 extinction coefficients at 665 nm and 860 nm [30]. The saturated oxygenation level is then computed with Eq. (6.9).

Processing implementations

The general processing workflow described above was performed using different computation environments:

-MATLAB CPU (Reference) processing was used with the linear optical properties LUT along with non-linear interpolation (griddata function in MATLAB) as a reference measurement since this is the most common implementation used today.

-MATLAB CPU processing was also used with the linear diffuse reflectance LUTs along with either bilinear interpolation (interp2 function in MATLAB) on the low density LUT, or nearest value (by rounding [25]) on both low density and high density LUTs.

-MATLAB GPU processing was used with the linear diffuse reflectance LUTs along with either bilinear interpolation (interp2 function in MATLAB) on the low density LUT, or nearest value (by rounding [25]) on both low density and high density LUTs. This particular implementation uses the MATLAB Parallel Computing Toolbox.

-C CUDA GPU processing was used with the linear diffuse reflectance LUTs along with either bilinear interpolation on the low density LUT, or nearest value (by rounding [25])

on both low density and high density LUTs. This particular implementation uses our custom-made C CUDA GPU code.

Note that all processing uses single precision with the exception of MATLAB CPU reference since the griddata function only works with double precision. A detailed description of each processing implementation is described below.

MATLAB CPU processing

The

MATLAB GPU processing

The MATLAB GPU processing uses the same code as the MATLAB CPU implementation but is specifically processed on the GPU. In fact, the previously described MATLAB CPU implementation is fully compatible with the MATLAB Parallel Computing Toolbox and the MATLAB gpuArray function is used to copy all the input variables in MATLAB workspace to the GPU at the beginning of the process.

C CUDA GPU processing

Our custom processing code in C CUDA is made with Microsoft Visual Studio and NVIDIA CUDA toolkit environment with basic operations and the cuFFT library of NVIDIA.

First, it is necessary to reserve GPU memory spaces for processing and to configure the system:

(1) a memory space is reserved for cosine and sine coefficients needed for the computation of the DFT, (2) a memory space is reserved for the pre-computed filtering masks DC and AC, (3) a ring buffer with a dimension of N+1 images is setup for the camera output images, (4) a memory space is reserved for the resulting images, (5) a cufftHandle is created and uses cufftPlanMany to configure a complex to complex (C2C) 2D FFT processing plan, (6) constant memory is used for scalar coefficients, (7) texture memory is allocated to store the LUT, and finally (8) memory space is assigned to store intermediate results. After accomplishing memory allocation, the coefficients for the DFT and the filtering masks are then computed on the CPU and they are transferred to the GPU memory. Once the initialization is finished, the processing loop can begin.

In the first task, each time an image is acquired, it is stored in the ring buffer of size N+1. The temporal demodulation kernel is then called to compute the DFT as explained above. This function directly returns an image for each demodulated wavelength. These results are stored as cufftComplex variable.

In the second task, the forward 2D FFT function cufftExecC2C is called with the created 2D plan, followed by the filtering kernel where FFT results are multiplied by the pre-computed filtering masks. Then the inverse 2D FFT is called with the same cufftExecC2C function. The MDC and MAC images are obtained by computing the complex magnitude of the previous results. Diffuse reflectance Rd,DC and Rd,AC are computed as explained above. Finally, the 3 different implementations to extract optical properties from linear diffuse reflectance LUT were performed with custom-made code. The bilinear interpolation on low density is computed using the following formula:

( )( ) ( )( ) ( )( ) ( )( ) , , , ,0 , , , ,0 , ,1 , ,0 , ,0 , ,1 , ,1 , ,1 ,, 1 ,, ,, 1 ,, ,, ,, 0 ,, 1 ,, ,, 1 ,, ,, ,, 0 , , , ,0 ,
, , ,0 Processing times depend strongly on the GPU configuration. We create here one thread for each pixel that we group into blocks of 64 threads, giving 100 % GPU occupancy, since for all the custom-made kernels the pixels pxy are independent in time and in space.

dD Cx dD C d AC x d AC d DC d AC d DC d AC d DC d AC RR R d DC d DC x d AC d AC x R R dD Cx dD C dA C dA Cx R R dD C dD Cx dA Cx dA C R R d DC x d DC d AC x d AC R RR RR R R RR RRRR RRRR mm m m m = -- + -- + -- + -- (6

Analysis

To evaluate the performance of the different implementations, we measured a hand as an in vivo complex sample positioned on a tissue-simulating phantom with a projected spatial frequency of 0.3 mm -1 for 2 seconds with a rate of 100 Hz on the system described above [24]. Processing was | Chapter 6. Real-time optical properties and oxygenation imaging using GPGPU 117 performed on the acquired images in their native format (1024 x 1024) and as well as with a smaller image size (512 x 512) obtained using software 2 x 2 binning. The resulting images (absorption, scattering and oxygenation) were compared with the reference MATLAB CPU results. These comparisons were done by measuring the mean percentage error given by the following formula: ).

() 11 ,, 00 , 1 %
In MATLAB, a stopwatch timer is used to measure performance with the "tic" "toc" function. For the GPU processing, NVIDIA Nsight environment is used to trace the CUDA code timing performance. It is important to note that we report processing time only, and not transfer time. For GPU processing, the CPU -GPU transfer time is 338 µs for 1024 x 1024 images and 86 µs for 512 x 512 images. errors compared with the reference processing (as described in equation ( 11)) are also shown in the third row. Finally cross section line profiles are shown in the fourth row, in solid blue for the reference MATLAB CPU processing and in dash red for the C CUDA GPU processing.

Results

The main results obtained from the processing methods described above are the demodulated wavelength images, the absorption and reduced scattering coefficients recovered for each wavelength and the saturated oxygenation level. First, absorption, reduced scattering and oxygenation images obtained by the reference MATLAB CPU implementation and the C CUDA GPU implementation are reported in Fig. 6.3.. In addition, the mean percentage error images and cross-section lines profiles are shown. One can notice good agreement between the C CUDA GPU results and the MATLAB CPU reference results, with maximum mean percent errors less than 0.1 %. images. In summary, compared to the reference implementation (linear optical properties LUT), all implementations with linear diffuse reflectance LUTs have a precision error of less than 1 %.

Best results are obtained with 2D interpolation on the low density LUT with a precision error less than 0.1 %, followed by the implementations using direct reading on the high density LUT with a precision error less than 0.2 % and finally, those using direct reading on the low density LUT with a precision error less than 1 %. In addition, we note globally a trend towards a positive percentage error with optical properties at 665 nm and a negative percentage error at 860 nm. However, the low level of the errors does not make this observation interpretable. Table 6.2. Processing time (in milliseconds) from wavelength demodulation to saturated oxygenation level measurement for 512×512 images.

The analysis of processing times for the different implementations is summarized in Table 6.2. for an image size of 512 x 512 pixels and Table 6.3. for an image size of 1024 x 1024 pixels. In summary, on our workstation, with the reference implementation (linear optical properties LUT)

and for an images size of 512 x 512, we can extract the absorption coefficient, the reduced scattering coefficient and determine the saturated oxygenation level in 9 seconds. It takes only a maximum of one hundred milliseconds for the CPU process when the linear diffuse reflectance LUT is used. This processing time can be divided by 5 and reach 20 milliseconds when the GPU computing power is used by MATLAB Parallel Computing Toolbox. The programming of the GPU in C CUDA that has been achieved in this work has allowed us to go even further and get | Chapter 6. Real-time optical properties and oxygenation imaging using GPGPU 120 these measurements in just 400 microseconds, 50 times faster than the MATLAB GPU process, 200 times faster than the MATLAB CPU process and 20,000 times faster than the standard method.

For 1024 x 1024 images, the processing was at best performed in only 1.6 milliseconds, 20 times faster than with the MATLAB GPU process, 200 times faster than with the MATLAB CPU process and 9,000 times faster than with the standard method. It can also be noted that for C CUDA implementations the processing time is almost linearly dependent of the image size. It should be noted that parallelized implementations on CPU are possible and would perform better than MATLAB CPU code. In this manuscript, we chose to focus our work on GPU implementations. Overall, with a C CUDA implementation using 2D interpolation on low density linear diffuse reflectance LUT the measurements can be obtained in 440 microseconds and 1.6 milliseconds respectively for images of sizes 512 x 512 and 1024 x 1024 with error rates below 0.1 %. It should be noted that acquisition time is the same for all processing methods, consisting of N = 10 frames acquired at 100 frames per second, hence an equivalent acquisition time of 100 milliseconds. Since we use a rolling window processing method for temporal demodulation, following the first 10 frames, results are calculated every frame, i.e. every 10 milliseconds. Therefore, fast processing is required to reach real-time capabilities (here processing time faster than 10 milliseconds, preferably 5 milliseconds).

Discussion

This work provides solutions to speed up processing times in SFDI. Particularly, it demonstrates that in the case of a multi-spectral acquisition (2 wavelengths in our case) using spatio-temporal modulation of light with a GPU, linearized diffuse reflectance LUT and single precision processing, we can extract optical properties and tissue oxygenation for a 1 Megapixel image in only 1.6 milliseconds. While the definition of real-time depends on the target application, the number of physiological parameters that will need to be processed will increase in the future. Our goal is to provide interpretable information to healthcare practitioners at a rate of at least 25 frames per second. Given that up to 6 physiological parameters could be extracted (oxy-hemoglobin, deoxy-hemoglobin, water, lipids, melanin, and yellow pigments [30]) both processing and acquisition time will need to be improved to allow to capture and process the necessary information.

However, this work is not without limitations. It would be expected that the processing time would increase linearly with the image size when using GPU processing. While this linearity is observed when using the C CUDA implementation, it is not the case when using the MATLAB Parallel Computing Toolbox. We believe this can be attributed to the difference in accuracy of the measurement tools between MATLAB and NVIDIA CUDA toolkit, especially when it comes to measuring very short processing times. It should also be noted that the MATLAB CPU/GPU implementation of the direct reading of the linear diffuse reflectance LUT we have written, is more suitable for CPU processing, and quite different from the one we wrote in C CUDA. Overall, with the processing times obtained for C CUDA implementation and the observed linear evolution, it can be predicted that it would be possible to process the optical properties of 6 wavelengths in less than 4.8 milliseconds and to derive 6 physiological parameters, with the acquisition time depending on the camera frame rate and the number of wavelengths used. However, scaling up to a larger number of wavelengths remains a significant challenge with the limited dynamic range of camera systems at high frame rates and with the necessity to balance the contributions of each wavelengths within the available dynamic range.

In this work we have demonstrated that a low density linear diffuse reflectance LUT (steps of 0.001) enables very good measurement accuracy by using simple bilinear interpolations. Yet, this density was chosen for the ease of handling. It should be expected that depending on the margin of error that can be tolerated, a lower density linear LUT could be used, and lead to shorter processing times. It can also be noticed that the accuracy of the measurements obtained for wavelength 665 nm is better than that obtained for the wavelength 860 nm. This difference is essentially related to a lower signal-to-noise ratio of the images obtained at 860 nm due to the lower quantum efficiency of the camera in this spectral range. Despite these limitations, we have shown by the obtained precision errors that single precision SSOP processing on linear LUT does not result from any significant loss in measurement accuracy.

Finally, for quantitative optical imaging of in vivo samples, it is necessary to take into account the sample profile since errors can reach values as high as 10% per cm. It is therefore necessary to integrate the sample profile for height correction into future implementations to have a robust, real-time measurement for clinical guidance [32].

Conclusion

In summary, we have presented in this chapter CPU and GPU implementations of the quantitative oxygenation measurement method we introduced in the previous chapter. The resulting accuracy and processing time of these implementations were compared. We demonstrated with our custommade GPU implementation an error rate lower than 1% compared to classic MATLAB implementation. We demonstrated also that we were able to extract optical properties (absorption and reduced scattering) maps for two wavelengths 665 nm and 860 nm as well as saturated oxygenation level for a 1 megapixel image in as low as 1.6 milliseconds processing time with GPU implementation directly written in C CUDA. This chapter concludes the workflow followed since Chapter 4 for the realization of a wide-field, real-time and quantitative measurement system of physiological and structural parameter for image-guided surgery. Indeed, the results obtained for the design of the oxygenation imaging system show that we can extend the system to more wavelengths and therefore measure more parameters of interest for surgeons such as lipids, water, scattering power and scattering amplitude. While the entire system workflow has been built, there are still key points such as the accuracy of the results and the visual quality of the images that can still be subject of improvements. Indeed, we have not yet considered the correction of the 3D profile of the samples; also, although the filters designed as described Chapter 4 are well performing, unpleasant artifacts on the edges remain on the final images. Addressing these two | Chapter 6. Real-time optical properties and oxygenation imaging using GPGPU 123 points is therefore essential for the use of the final system during surgical procedures. For this, in the next chapter we will discuss about the new method we developed to obtain higher visual quality and accuracy based on deep learning and 3D profile correction: a proposed method which still uses GPGPU implementation for real-time processing and visualization.

129

The previous three chapters were devoted to (1) improving the accuracy and visual quality of SSOP images by designing 2D anisotropic filters, (2) to the development of a multispectral acquisition method based on spatiotemporal light modulation that allowed us to design a system for oxygenation level measurement and (3) finally to the real-time implementation on GPU for the system guaranteeing the achievement of our real-time processing and visualization objective.

However, far from being perfect, the system still exhibits imperfections related to the filtering technique and due to not taking the 3D profile of the sample into consideration during our developments. In this chapter, which is our most advanced contribution, we will present a new method that combines convolutional neural networks and 3D profile correction to achieve highquality and accuracy of optical properties imaging. The whole being realized in real-time thanks to the custom implementation on C CUDA GPU of the entire optical properties' extraction process.

The material for this chapter is our publication "Real-time, wide-field and high-quality single snapshot imaging of optical properties with profile correction using deep learning" in Biomedical Optics Express, Volume 11, Issue 10.

7.

Context and motivations

Spatial Frequency Domain Imaging (SFDI) is a relatively inexpensive diffuse optical imaging method allowing quantitative, wide-field and rapid measurement of optical properties [1][2][3][4]. The capabilities of SFDI have been demonstrated across a large number of biomedical applications including small animal imaging, burns wound assessment, skin flap monitoring, early stage diabetic foot ulcers diagnosis, among others [5][6][7][8][9]. However, because the standard implementation of SFDI requires the sequential acquisition of several sinusoidal patterns of light, its potential for clinical applications that require real-time feedback, such as intraoperative surgical guidance, is limited.

With the aim of deploying such technology in surgical settings, Single Snapshot imaging of Optical Properties (SSOP) was introduced as a new approach for performing SFDI in real-time.

Conceptually, SSOP works by reducing the total number of necessary acquisitions to a single image -thereby increasing the speed of acquisition. First proofs of concept demonstrated the feasibility of real-time optical properties imaging (i.e. faster than 25 images per second) with high accuracy (i.e. with less than 10% error compared with SFDI). However, the approach initially

Materials and methods

Spatial Frequency Domain Imaging

State of the art in SSOP

The standard SFDI method requires the sequential acquisition of several spatially-encoded images for the retrieval of optical property maps. The acquisition is typically performed using at least at two different spatial frequencies (e.g., fx = 0 and 0.2 mm -1 ) and 3 phase shifts or more for improved accuracy and visual quality (e.g. 7 phases) [4,21]. The patterns are generated and projected using a Digital Micromirror Device (DMD) and captured with a camera as depicted in Fig. 7.

(A).

Following the acquisition, the images are used to obtain the amplitude modulation at each spatial frequency, e.g. MDC for fx = 0 mm -1 and MAC for fx = 0.2 mm -1 . From these measurements, the diffuse reflectance images Rd,DC and Rd,AC are obtained through calibration with a material of known optical properties. Finally, the measured reflectance allows for the pixel-wise estimation of the absorption and reduced scattering coefficients (µa, µs') across the image via the resolution of the inverse problem -such as with precomputed high density lookup tables (LUTs) generated using White Monte Carlo simulations [22]. In contrast to conventional SFDI, which necessitates a minimum of six images for quantification, SSOP is an alternative SFDI implementation capable of fast optical property extraction via acquisition of a single high spatial frequency image (e.g., fx = 0.2 mm -1 ) -reducing the total acquisition time substantially by a few-fold, enabling real-time acquisition [4]. Amplitude modulations are extracted by filtering directly in the Fourier domain with both a low and high pass filter [10][11][12]. The general workflow is then similar to that of SFDI (illustrated in Fig. 7.

(B)).

Of importance, one should note that demodulation is the key step with which the resulting accuracies of both OP map retrieval and visual quality depend. Details of the acquisition, demodulation and processing typically used in SFDI and SSOP are provided in [3,4], and details about creating a custom SFDI system and processing code are available on www.openSFDI.org and www.healthphotonics.org [23].

SSOP demodulation

Typically, SSOP demodulation is performed either line by line [10], or in 2D spatial Fourier domain using rectangular low and high pass filters [11]. Recent developments in 2D spatial Fourier domain include using anisotropic low pass filtering based on sine windows and anisotropic high pass filtering based on Blackman windows [12]. This approach reported accuracies under 8.8% in absorption (µa) and 7.5% in reduced scattering (µs'). While this SSOP implementation was successfully used to perform real-time oxygenation measurements [24], it is still affected by edge discontinuity, has a lower resolution than standard SFDI and does not include profile correction for the sample.

3-Dimensional profile correction applied to SFDI

Optical properties measurements of non-flat samples through SFDI conventionally suffer from inaccuracies due to intensity variations across the sample surface. To account for these errors, 3D profile correction methods for SFDI have been implemented by measuring the profile using structured-light profilometry [25,26] and correcting for its effect [27]. Several correction methods have been proposed for SFDI, all using either phase shifting profilometry or Fourier-based profilometry. Overall, they have illustrated high correction accuracy for large angles of tilt, from 40 up to 75 degrees depending on the method [11,[27][28][29]. However, the use of a single image SSOP acquisition approach leads to lower profile accuracy that significantly affects the accuracy 133 of retrieved OPs [11]. The main focus of the work herein is to leverage a novel deep learning approach to enable profile-corrected SSOP imaging in real-time and for complex non-flat samples.

Deep learning approaches to SFDI

State of the art

Several groups have focused their efforts on developing deep learning and machine learning approaches for SFDI or SSOP to improve accuracy, visual quality and/or speed. A random forest regression algorithm was developed for estimating optical properties in the spatial frequency domain, trained on Monte-Carlo simulations data [30]. The algorithm "learns" the nonlinear mapping between diffuse reflectance at two spatial frequencies, and optical properties of the tissue, with the goal of reducing the computational complexity when a high density pre-computed LUT is used. Using this method, optical properties could be obtained over a 1-megapixel image in 450 ms with errors of 0.556% in µa and 0.126% in µs'.

Following, a DNN-based framework was proposed for multi-spatial frequency processing SFDI [31]. This network, a Multi-Layered Perceptron (MLP), outperformed other classic inversion methods (iterative method and nearest search) in terms of speed with computation times of 5ms for 100x100 pixels and 200ms for 696x520 pixels, using a 5 spatial frequencies SFDI acquisition and substantially reducing measurement uncertainties compared to the widely used 2 frequency SFDI processing method. More recently, CNNs have been introduced for diffuse optical imaging owing to their remarkable performance enhancements across imaging applications compared with MLPs [32]. In the case of SFDI, a Generative Adversarial Networks (GANs) framework (named "GANPOP") has been reported for the retrieval of optical properties directly from single input images with high accuracy [33]. In this work, SFDI was used to obtain ground-truth optical properties and the model developed took approximately 40ms to process a 256x256 pixel image using an NVIDIA Tesla a P100 GPU Accelerator.

Both the random forest and MLP works have a commonality in that they both estimate optical properties from diffuse reflectance -though, GANPOP bypasses the demodulation step and retrieves optical properties directly from calibrated SSOP input. In addition, all these methods have been developed with a focus on providing a more accurate optical properties imaging ability, and though GANPOP illustrates partially profile-corrected OP retrieval, no DL-based workflow to date retrieves profilometry directly. Further, both the MLP and GANPOP's generator are comprised of a relatively large parametric size (on the order of 1x10 6 and 8x10 7 weights, respectively). In turn, these relatively large DL models can be both computationally burdensome and memory intensive, limiting their use for real-time applications. In addition, though GANs have been increasing in popularity across applications of computer vision, they are notoriously problematic with regards to hyperparameter tuning, stable convergence and computational demand for proper model training [34]. Finally, no prior DL-based approach has actually been applied in real-time.

The proposed method

Inspired by these previous developments and results, we propose herein to combine GPU accelerated computing [24] with efficient DL architectures applied for a seamless and user-friendly SSOP processing methodology -allowing for real-time image reconstruction of mega-pixel images with high accuracy and visual quality along with direct retrieval of profilometry. The method consists of the following steps:

U-NET 1: Demodulation: A first small-size U-NET (less than 2e4 parameters) is used to get demodulation images (MDC and MAC), instead of directly map raw input data to optical properties.

Since the network has dedicated weights according to the spatial frequency, both higher visual quality and more accurate optical properties can be obtained, in contrast with a Fourier filtering method with fixed spatial frequency parameters.

U-NET 2: Profilometry: A second small-size U-NET (less than 2e4 parameters), identical in structure to the first one, is used to get wrapped phase and obtain the profile, avoiding the use of Fourier domain 3D profilometry that leads to poor accuracy due to its dependency to spatial frequency. In our implementation, the real and imaginary part are extracted from the network to compute the wrapped phase by mean of atan2() function [35].

Real-time profile-corrected optical properties imaging: To reduce the inference time of both U-NETs, a custom-made GPU implementation of the networks was undertaken using the CUDA deep learning library cuDNN. Following the use of the two computationally inexpensive U-NETs described above, (used for demodulation and profilometry, respectively), the extraction of optical properties can be performed along with profile correction at high speed using a previously developed real-time optical properties GPU processing code [24]. The use of the GPU-accelerated 135 twin-U-NET framework detailed above, followed by a GPU-optimized precomputed lookup table

inversion method has a strong potential for delivering real-time and high visual quality SSOP imaging capabilities to many who would otherwise be unable to utilize more computationally burdensome approaches. The general network architecture is illustrated in Fig. 7.2. and is inspired by the U-Net architecture [36]. It is a n-dimensional network ("Line U-NET" when n = 1 and "Image U-NET" when n = 2).

It consists of a contracting path of 4 stages (left side) and an expansive path of 4 stages (right side).

The expansive path is a fork and has 2 outputs. The contracting path consists of the repetition of a 7x7 n-1 convolution layer (padded convolution), each followed by a Tanh activation layer and a 2x2

n-1 max pooling layer (with stride 2). The number of feature channels is double at each level starting by 1 feature channel at stage 1. Every stage in the expansive path consists of an up-sampling 2x2 n-1 layer followed by a 3x3 n-1 deconvolution layer that halves the number of feature channels, a concatenation with the corresponding feature map from the contracting path, and a 7x7 n-1 136 convolution layer followed by a Tanh activation layer. At the final stage, a 7x7 n-1 convolution layer with linear activation is employed to retrieve the output. In total, the twin U-NET architecture is comprised of either 3064 or 18886 weights (according to n equal to 1 or 2, respectively) for 17 convolutional layers. The python code for each architecture is available at www.healthphotonics.org in the "resources" section.

Training

The network was implemented with Tensorflow using Keras backend and is compiled with mean squared error (MSE) as loss function and Adam as optimizer [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF].

GPU C CUDA implementation

Our custom SSOP processing code (using C CUDA) is made with Microsoft Visual Studio and NVIDIA CUDA toolkit environment. Based on the previous developed code [24], our work focused on the GPU implementation of the novel demodulation via the previously described twin U-NET workflow. To perform model inference, there are three options: using the training framework itself (convention), using TensorRT [39] or writing a custom code to execute the network using cuDNN low-level libraries and math operations [START_REF]cuDNN Developer Guide[END_REF][START_REF]cuDNN API[END_REF]. We choose the latter for its high efficiency compared to other solutions. Since the network is a succession of many convolution, max pooling, up-sampling and concatenation, layers will be presented only once.

Convolution

Before one uses cuDNN functions to describe and execute a network, a handle must be created using cudnnCreate. A convolution consists of a description part and an execution part.

-Description part: (1) an input tensor is created with cudnnCreateTensorDescriptor and cudnnSetTensor4dDescriptor, this tensor memory space is then allocated; (2) the filter is described with cudnnCreateFilterDescriptor and cudnnSetFilter4dDescriptor, a memory space is allocated and the weights loaded; (3) the type of convolution is specified using cudnnCreateConvolutionDescriptor and cudnnSetConvolution2dDescriptor, to get the convolution output dimension cudnnGetConvolution2dForwardOutputDim is called; (4) a second tensor is created as previously for the convolution output; (5) the activation is set with cudnnCreateActivationDescriptor and cudnnSetActivationDescriptor; (6) a third tensor is created as previously to store the output after activation; (7) the optimized convolution algorithm is found with cudnnGetConvolutionForwardAlgorithm and the required work space for its execution with cudnnGetConvolutionForwardWorkspaceSize, this space is then allocated.

-Execution part: cudnnConvolutionForward is called first to execute the convolution, follow by a custom made add_bias_gpu kernel. Finally, to finish the execution, an activation function cudnnActivationForward is called.

Max pooling

The max pooling operation is also comprised of a description and an execution part: (1) it requires an input tensor (e.g., the previous activation tensor); (2) a description of the pooling is created via cudnnCreatePoolingDescriptor and cudnnSetPooling2dDescriptor and the pooling output dimension is defined with cudnnGetPooling2dForwardOutputDim; (3) a tensor is pre-allocated to store the output after pooling; (4) lastly, the max pooling is performed by calling cudnnPoolingForward. SSOP DL: For the two SSOP DL methods, while still affected at low spatial frequency (0.1 mm - 1 ), the image degradation is significantly reduced compared to the SSOP Filtering approach. A much better definition of the hands and organs shape and a reduction of the edge discontinuity effect can be clearly noticed (supported by SSIM measurements, Table 1

Quantitative analysis

Mean absolute percentage error of the full image for all the 24 samples SSOP error maps are computed verses 7-phase SFDI and reported in Table 7.2. (in mean ± standard deviation) for quantitative assessment.

Profile: Generally, 3D profile errors decrease when spatial frequency increases: SSOP DL Conv 7x7 demonstrated the best results with errors as low as 9.5 ±5.9% followed by SSOP Filtering that gave errors as low as 9.8 ±5.8% and SSOP DL Conv 7x1 with errors as low as 10.2 ±6.4%.

Profile-corrected absorption: Profile-corrected absorption coefficients illustrate a decreasing error trend along with spatial frequency until 0.2mm -1 and an increase at 0.3mm -1 : SSOP DL Conv 7x7 still performs best, with errors as low as 7.7 ±2.8%, followed by SSOP DL Conv 7x1 with errors as low as 8.6 ±7.0% and SSOP Filtering with errors as low as 8.7 ±2.6%. ). In this movie, one can notice the edge discontinuity effect when using SSOP Filtering and how it depends on the hand position.

In contrast, the SSOP Deep Learning workflow presented herein significantly reduces these effects and reconstructs OPs with high accuracy. Thus, for the first time, this work illustrates the potential of SSOP for use in challenging settings where motion, depth and real-time feedback are critical.

Discussion

Single snapshot imaging of optical properties was initially designed to provide fast acquisition and rapid processing but suffered from degraded visual quality compared with SFDI. Leveraging the recent developments in Deep Learning together with GPU-accelerated computation we demonstrated herein that SSOP is able to achieve reconstructive performances comparable to 7phase SFDI with concurrent sample profile correction, leading to errors as low as 7.5 ±2.7% for a complex object such as a hand and ex vivo pig organs. Such technology can readily be deployed in the intraoperative room to guide surgeons by offering real-time visualization feedback of functional and structural biomarkers.

Of note, it is important to highlight that real-time capabilities were achieved thanks to the use of a relatively simplistic CNN architecture (comprising less than 2x10 4 parameters) coupled with low-level GPU implementation. Despite the use of a small size DL architecture, our approach allowed pushing the limits of SSOP to match data accuracy of SFDI. To obtain such result, the method combines the advantages of U-NET network demodulation and rapid GPU processing. It is however interesting to note that simpler implementations such as SSOP Filtering demonstrate good results despite edge artifacts with reconstruction errors as low as 8.7 ±2.6% and 9.3 ±2.5%

for absorption and reduced scattering, respectively. Though, such results using SSOP Filtering are possible only when state-of-the-art 2D anisotropic filtering is properly implemented, as described in prior work [12].

Another interesting observation is that using the SSOP Deep Learning methods described here, accurate and good visual quality imaging of complex object can be performed at spatial frequencies as low as 0.2mm -1 . Although not demonstrated, the use of frequencies that are typically lower than 0.3mm -1 seems preferable to allow deeper penetration of the photons and better signal to noise ratio (amplitude modulation decreases significantly with spatial frequency). As such, the use of the presented SSOP DL methodology allows for a more reliable use of SFDI methods for healthcare application that necessitate real-time employability.

However, this work is not without limitations. The performance of a trained network always depends on the training dataset. In our case, we avoided the 3 phases SFDI acquisition signal to noise ratio limitation when spatial frequency increases by using 7 phases SFDI acquisition. Despite this precaution, SFDI ground truth at spatial frequencies 0. Another classical limitation regarding the datasets used for training in DL is the variety of configurations offered by the sample. Including more images of organs having different shapes and edge discontinuities would likely improve the network accuracy. Indeed, this was observed in previously published work [33] and noted as crucial for model generalizability.

This study highlights also the trade-off that exists between the choice of the spatial frequency used for acquisition and the capacity of SSOP to provide both accurate and high visual quality images. One can notice for all the SSOP Filtering images obtained at low frequency that the capability to retrieve the structure of the hand is limited by the half-period of the spatial frequency used.

Overall, this study demonstrates that with the combination of a simple (small-size, 2x10 4 parameters) DNN and low-level GPU implementation, accurate and high visual quality profilecorrected SSOP can be performed in real-time for one wavelength. Since our goal is to provide interpretable information to healthcare practitioners, additional wavelengths will need to be acquired to obtain interpretable physiological parameters such as oxy-hemoglobin, deoxyhemoglobin, lipids or water content. With a 18ms processing time using the best performing method (DL conv 7x7), at present only 2 wavelengths could reasonably be processed to meet realtime definition (i.e. frame rate faster than 25 frames per second). However, due to the rapid increase in GPU computing power and their low price, we expect that capabilities can be extended easily to process more than 2 wavelengths.

Finally, all code and data used for this work is available at www.healthphotonics.org under the resources tab. This is an open initiative for transparency and with the objective to share and accelerate developments in the field of SFDI. Sources codes, raw data as well as executable code readily usable can be downloaded. Other resources such as www.openSFDI.org can be used to fabricate a SFDI system that can be used with the codes provided here [23].

In summary, this work lays the foundation for real-time multispectral quantitative diffuse optical imaging for surgical guidance and healthcare applications.

Conclusion

As described in the previous three chapters we had succeeded in designing a multispectral system for the measurement of physiological and structural parameters of tissues. However, the system obtained was limited in accuracy and overall visual image quality. In this chapter, using a simple (less than 2x10 4 parameters) convolutional deep learning architecture and low-level GPU implementation, we have been able to push the limits of SSOP to match data quality of SFDI with a real-time capability. In particular, we have shown that by designing a computationally inexpensive twin U-Net architecture for both the demodulation stage and the sample profile extraction stage, along with a custom-made GPU-accelerated implementation using C CUDA, megapixel profile-corrected optical properties (absorption and reduced scattering) maps were extracted in 18ms with high image quality and preserving a high accuracy. The method was validated on anatomical relevant sample such as in vivo hand and ex vivo pig organs. The mean absolute percentage error relative to SFDI method is less than 10% in both profilometry and profile-corrected optical properties. This novel processing framework lays the foundation of the next version of the real-time quantitative oxygenation imaging system for surgery.

During the three years of my PhD I worked on the development of an imaging instrument for surgical guidance using diffuse optical imaging. Our objective was to develop an intraoperative image-guided surgery device based on the spatial frequency domain imaging method SFDI capable of providing physiological and structural information essential to surgeons such as oxygenation, lipids, melanin and water levels, together with scattering power and scattering amplitude. To achieve our goal, we combined the fast acquisition SFDI method named SSOP with temporal modulation of light to characterize the optical properties of tissues acquired at several wavelengths simultaneously. In addition, in order to achieve a clinical application, we had to address the issue of SSOP image degradation and loss in accuracy, design a real-time wide-field multispectral imaging acquisition method and develop innovative real-time processing techniques to guarantee real-time processing, and high quality visualization.

In Chapter 4 we discussed the first set of contributions of this work: new filtering strategies for improving the image quality of optical properties maps obtained using SSOPs. In fact, legacy SSOP methods which use rectangular DC and AC filters are limited by a degraded image quality, resulting in artifacts that diminish its potential for clinical use. Therefore, we investigated the effect of anisotropic two-dimensional filtering strategies based on other well-known filtering windows such as Sine, Hann, Blackman, and Gaussian. This study was performed for spatial frequencies ranging from 0.1 to 0.4 mm -1 onto N = 10 hands. Both results on accuracy and image quality were quantified in comparison with standard multiple image acquisitions in the spatial frequency domain. Overall, the equations enabling the design of anisotropic DC and AC filters were provided and consequently, the best optimized filtering windows combination at each spatial frequency were presented. One of the best combinations, i.e. Sine DC and Blackman AC windows, provide significantly improved image quality while preserving accuracy in optical properties measurements, in particular at high spatial frequencies (typically, 0.4 mm -1 ) with mean errors in predicting optical properties remaining under 9.5% in absorption and 8.5% in reduced scattering again respectively 11% and 9.5% for legacy filtering. This study also demonstrated the dependency of image quality and accuracy on spatial frequency.

The second set of contributions of our work described in Chapter 5, relates to the design, validation and application of the multispectral acquisition method using spatiotemporal modulation of light.

With this method, we are able to perform real-time, wide-field (13 x 13 cm²) and quantitative imaging of blood parameters in vivo. In fact, the advanced SSOP method presented in Chapter 4 allow us for the wide-field imaging of optical properties but at a single wavelength. However, for this technology to be useful for clinicians, images must be displayed in terms of metrics related to the physiological and structural state of the tissue, hence interpretable in order to guide decisionmaking. At the acquisition stage, the presented multispectral approach relies on the modulation of light sources in time to encode wavelengths and on the modulation of the projection of these sources in space to extract the optical properties of living tissues. At the processing stage instead, it is based on a temporal demodulation process using Discrete Fourier Transform (DFT) to retrieve each wavelength contribution followed by SSOP processing. The method was validated on a system build with two wavelengths (665 nm and 860 nm), and therefore capable of providing optical properties at these two wavelengths, together with saturated oxygenation level. This validation was performed by comparing the results of the method with those obtained for a SFDI sequential acquisition during two experiments, one with a phantom and the other one with an in vivo hand. The results showed a maximum error rate of 4.2%. Finally, we carried out an arm cuff occlusion experiment confirming the potential for clinically relevant use such as ischemia detection.

Chapter 6 addressed the real-time implementation step that is essential to the success of our project.

In fact, the spatiotemporal modulation technique developed is compatible with real-time acquisition (100 Fps for our system); it however requires to be supported by a high computing power to ensure real-time processing and visualization. This was performed by means of a custom low-level implementation in C CUDA (Compute Unified Device Architecture) on GPU for the oxygenation imaging system. The implementation was compared to conventional MATLAB CPU/GPU implementations and showed an error rate of less than 0.1% for the extraction of optical properties. Finally, the most important result obtained in this sense is the ability to measure in just This work has been accepted by the scientific community and has been the subject of 6 publications in international journals and 4 presentations at international conferences. The list is provided below in the scientific contribution section.

Limitations

The work we have done has surely some limitations as well. The first is the quality of the reference used for all our studies. If the SFDI method is theoretically more accurate than SSOP, it has also some limits, that can be summarized as: the need for a multi calibration for more accuracy due to the variety of optical properties present on an in vivo sample; the decline in the signal-to-noise ratio as the spatial frequency increases which has certainly put a bias in our comparisons in Chapters 4 and 7. The second frame of limitations to the quality of our work are represented by those related to the spatiotemporal modulation acquisition method. In fact, in our results we highlighted the low quality of the images obtained for the 860 nm wavelength which is due to the lower quantum efficiency of the camera at this wavelength (lower SNR). This lower efficiency at 860 nm influences the quality of the optical properties and the measurement of the oxygenation rate. Another limitation of the spatiotemporal modulation method is related to the maximum number of simultaneously adjustable wavelengths, limited by the dynamic range of the camera.

Moreover, the rolling window processing approach comes at the cost of a short delay between the acquired scene and the displayed information. Thirdly, on the new SSOP method based on deep learning, we can note the need to gather more varied acquisitions of in vivo samples to improve the precision of the CNN networks we implemented. Finally, if the deep learning solution for SSOP is a great asset, it naturally generates a processing time, 16 times slower, which reduces the number of flexible wavelengths before making our GPU architecture evolve. Rectangular bandpass

| Appendix xxvi

Appendix 1: Filtering equations and comparison mean percentage error results

( , ) = (1) × (1) for -( × ) ≤ ≤ ( × ) -( × ) ≤ ≤ ( × ) = 1.6 Exponential bandpass ( , ) = ---1 2 12 8.69 × ---1 2 12 8.69 for -( × × 1.8) ≤ ≤ ( × × 1.8) -( × ) ≤ ≤ ( × ) = 2.5
Triangular bandpass

( , ) = 1 - - -1 2 -1 2 × 1 - - -1 2 -1 2 for -( × × 3) ≤ ≤ ( × × 3) -( × ) ≤ ≤ ( × ) = 3.5 Welch bandpass ( , ) = 1 - - -1 2 -1 2 2 × 1 - - -1 2 -1 2 2 for -( × × 1.75) ≤ ≤ ( × × 1.75) -( × ) ≤ ≤ ( × ) = 2 Sine bandpass ( , ) = N -1 × m -1 for -( × × 1.88) ≤ ≤ ( × × 1.88) -( × ) ≤ ≤ ( × ) = 2.25 Lanczos bandpass ( , ) = 2 N -1 -1 × 2m -1 -1 for -( × × 2.1) ≤ ≤ ( × × 2.1) -( × ) ≤ ≤ ( × ) = 2.4
Hann bandpass

( , ) = 2 N -1 × 2 m -1 for -( × × 2.5) ≤ ≤ ( × × 2.5) -( × ) ≤ ≤ ( × ) = 2.5
Hamming bandpass Blackman-Nuttall bandpass

( , ) = 0 -1 2 N -1 × 0 -1 2 M -1 for -( × × 2.4) ≤ ≤ ( × × 2.4) -( × ) ≤ ≤ ( × ) 0 = 0.54 ; 1 = 1 -0 = 2.7 Blackman bandpass ( , ) = 0 -1 cos 2 n -1 + 2 cos 4 n -1 × 0 -1 cos 2 m -1 + 2 cos 4 m -1 for -( × × 3) ≤ ≤ ( × × 3) -( × ) ≤ ≤ ( × ) = 0.16; 0 = 1- 2 ; 1 = 1 2 ; 2 = 2 = 3 Nuttall bandpass ( , ) = 0 -1 cos 2 n -1 + 2 cos 4 n -1 -3 cos 6 n -1 × 0 -1 cos 2 m -1 + 2 cos 4 m -1 -3 cos 6 n M-1 for -( × × 3.66) ≤ ≤ ( × × 3.66) -( × ) ≤ ≤ (
( , ) = 0 -1 cos 2 n -1 + 2 cos 4 n -1 -3 cos 6 n -1 × 0 -1 cos 2 m -1 + 2 cos 4 m -1 -3 cos 6 n M-1 for -( × × 3.6) ≤ ≤ ( × × 3.6) -( × ) ≤ ≤ ( × ) 0 = 0.3635819; 1 = 0.4891775; 2 = 0.1365995; 3 = 0.0106411 = 3.4 Gaussian bandpass ( , ) = -1 2 --1 2 -1 2 2 × - 1 2 --1 2 -1 2 2 for -( × × 2.7) ≤ ≤ ( × × 2.7) -( × ) ≤ ≤ ( × ) = 0.4 = 3
, ) = 1 - - -1 2 -1 2 × 1 - - -1 2 -1 2 for (( × ) -( × × 0.8)) ≤ ≤ (( × ) + ( × × 0.8)) -( × ) ≤ ≤ ( × ) = 1.2 Welch bandpass ( , ) = 1 - - -1 2 -1 2 2 × 1 - - -1 2 -1 2 2 for (( × ) -( × × 0.8)) ≤ ≤ (( × )+ ( × × 0.8)) -( × ) ≤ ≤ ( × ) = 1.4 Sine bandpass ( , ) = N -1 × m -1 for (( × ) -( × × 0.8)) ≤ ≤ (( × ) + ( × × 0.8)) -( × ) ≤ ≤ ( × ) = 1.4 Lanczos bandpass ( , ) = 2 N -1 -1 × 2m -1 -1 for (( × ) -( × × 0.9)) ≤ ≤ (( × ) + ( × × 0.9)) -( × ) ≤ ≤ ( × ) = 1.4 (
Hann bandpass Gaussian bandpass

( , ) = 2 N -1 × 2 m -1 for (( × ) -( × × 0.9)) ≤ ≤ (( × ) + ( × × 0.9)) -( × ) ≤ ≤ ( × ) = 1.3 Hamming bandpass ( , ) = 0 -1 2 N -1 × 0 -1 2 M -1 for (( × ) -( × × 0.8)) ≤ ≤ (( × ) + ( × × 0.8)) -( × ) ≤ ≤ ( × ) 0 = 0.54 ; 1 = 1 -0 = 1.3 Blackman bandpass ( , ) = 0 -1cos 2 n -1 + 2cos 4 n -1 × 0 -1cos 2 m -1 + 2cos 4 m -1 for (( × ) -( × × 1)) ≤ ≤ (( × ) + ( × × 1)) -( × ) ≤ ≤ ( × ) = 0.16; 0 = 1- 2 ; 1 = 1 2 ; 2 = 2 = 1.3 Nuttall bandpass ( , ) = 0 -1 cos 2 n -1 + 2 cos 4 n -1 -3 cos 6 n -1 × 0 -1cos 2 m -1 + 2cos 4 m -1 -3cos 6 n M-1 for (( × )-( × × 1.1)) ≤ ≤(( × )+( × × 1.1)) -( × ) ≤ ≤ (
( , ) = - 1 2 --1 2 -1 2 2 × -1 2 --1 2 -1 2 2 for (( × )-( × × 0.8)) ≤ ≤ (( × )+ ( × × 0.8)) -( × ) ≤ ≤ ( × ) = 0.4 = 1.3 Rectangular highpass ( , ) = (1) × (1) for (( × ) -( × × 0.32)) ≤ ≤ ( × ) -( × ) ≤ ≤ ( × ) = 1.3 Exponential highpass ( , ) = ---1 2 12 8.69 × ---1 2 12 8.69 for (( × ) -( × × 0.6)) ≤ < ( × ) -( × ) ≤ ≤ ( × ) (1) × (1) for ( × ) ≤ ≤ ( × ) -( × ) ≤ ≤ ( × ) = 1.3 Triangular highpass ( , ) = ⎩ ⎪ ⎨ ⎪ ⎧ 1 - - -1 2 -1 2 × 1 - - -1 2 -1 2 for (( × ) -( × × 0.8)) ≤ < ( × ) -( × ) ≤ ≤ ( × ) (1) × (1) for ( × ) ≤ ≤ ( × ) -( × ) ≤ ≤ ( × ) = 1.2 Welch highpass ( , ) = ⎩ ⎪ ⎨ ⎪ ⎧ 1 - - -1 2 -1 2 2 × 1- - -1 2 -1 2 2 for (( × )-( × ×0.8))≤ <( × ) -( × )≤ ≤( × ) (1) × (1) for ( × ) ≤ ≤ ( × ) -( × ) ≤ ≤ ( × ) = 1.4 Sine highpass ( , ) = ⎩ ⎪ ⎨ ⎪ ⎧ N -1 × m -1 for (( × ) -( × × 0.8)) ≤ < ( × ) -( × ) ≤ ≤ ( × ) (1) × (1) for ( × ) ≤ ≤ ( × ) -( × ) ≤ ≤ ( × ) = 1.4 | Appendix xxviii Table A1.2

. Bandpass and highpass AC filters equations with their cutoff frequency and size

Lanczos highpass • 32 LD/ST (load/store) units that are in charge of reading is writing through 64KB of shared memory that also serves as level 1 cache.

( , ) = ⎩ ⎪ ⎨ ⎪ ⎧ 2 N -1 -1 × 2m -1 -1 for (( × ) -( × × 0.9)) ≤ < ( × ) -( × ) ≤ ≤ ( × ) (1) × (1) for ( × ) ≤ ≤ ( × ) -( × ) ≤ ≤ ( × ) = 1.4 Hann highpass ( , ) = ⎩ ⎪ ⎨ ⎪ ⎧ 2 N -1 × 2 m -1 for (( × ) -( × × 0.9)) ≤ < ( × ) -( × ) ≤ ≤ ( × ) (1) × (1) for ( × ) ≤ ≤ ( × ) -( × ) ≤ ≤ ( × ) = 1.3 Hamming highpass ( , ) = ⎩ ⎪ ⎨ ⎪ ⎧ 0 -1 2 N -1 × 0 -1 2 M -1 for (( × ) -( × × 0.8)) ≤ < ( × ) -( × ) ≤ ≤ ( × ) (1) × (1) for ( × )≤ ≤( × ) -( × ) ≤ ≤ ( × ) 0 = 0.54 ; 1 = 1 -0 = 1.3 Blackman highpass ( , ) = 0 -1cos 2 n -1 + 2cos 4 n -1 × 0 -1cos 2 m -1 + 2 cos 4 m -1 for (( × ) -( × × 1)) ≤ < ( × ) -( × ) ≤ ≤ ( × ) (1) × (1) for ( × )≤ ≤( × ) -( × ) ≤ ≤ ( × ) = 0.16; 0 = 1- 2 ; 1 = 1 2 ; 2 = 2 = 1.3 Nuttall highpass ( , ) = 0 -1 cos 2 n -1 + 2 cos 4 n -1 -3 cos 6 n -1 × 0 -1 cos 2 m -1 + 2 cos 4 m -1 -3 cos 6 n M -1 for (( × ) -( × × 1.1)) ≤ < ( × ) -( × )≤ ≤( × ) (1) × (1) for ( × )≤ ≤( × ) -( × ) ≤ ≤ ( × ) 0 = 0.355768; 1 = 0.487396; 2 = 0.144232; 3 = 0.012604 = 1.3 Blackman-Harris highpass ( , ) = 0 -1 cos 2 n -1 + 2 cos 4 n -1 -3 cos 6 n -1 × 0 -1 cos 2 m -1 + 2 cos 4 m -1 -3 cos 6 n M -1 for (( × ) -( × × 1.1)) ≤ < ( × ) -( × )≤ ≤( × ) (1) × (1) for ( × )≤ ≤( × ) -( × ) ≤ ≤ ( × ) 0 = 0.35875; 1 = 0.48829; 2 = 0.14128; 3 = 0.01168 = 1.3 Blackman-Nuttall highpass ( , ) = 0 -1 cos 2 n -1 + 2 cos 4 n -1 -3 cos 6 n -1 × 0 -1 cos 2 m -1 + 2 cos 4 m -1 -3 cos 6 n M -1 for (( × ) -( × × 1.1)) ≤ < ( × ) -( × )≤ ≤( × ) (1) × (1) for ( × )≤ ≤( × ) -( × ) ≤ ≤ ( × ) 0 = 0.3635819; 1 = 0.4891775; 2 = 0.1365995; 3 = 0.0106411 = 1.3 Gaussian highpass ( , ) = ⎩ ⎪ ⎨ ⎪ ⎧ -1 2 --1 2 -1 2 2 × - 1 2 --1 2 -1 2 2 for (( × ) -( × × 0.8)) ≤ < ( × ) -( × )≤ ≤( × ) (1)
• 4 x 16384 32 bit status registers or 256KB.

• Texture memories that are read-only memories that have their own cache, optimized for access to 2D data.

• 128 CUDA cores (640 cores in total for the chip) each of which incorporate an ALU and an FPU both running in a pipeline.

• And finally 32 SFU units that assist them for special operations including the functions sin, cos, tan, atan, sqrt etc… [5].

Note that there is given the architecture only 32 active threads by SMM.

A2.2. Programming Interface (API) for GPGPU

At the beginning of the era of using GPUs to perform scientific calculations, developers had to divert graphical APIs such as OpenGL (Open Graphics Library), Direct3D from their primary function to realize the programs they needed, which was very tedious. In order to overcome these difficulties, over the years he has created various APIs dedicated to GPU computing. However, it was not until 2007 that the NVIDIA graphics card manufacturer made available the first programming interface dedicated to the GPGPU: CUDA (Compute Unified Device Architecture).

This solution was followed in 2008 by OpenCL (Open Computing Language), which is an open industrial standard, provided by the industrial consortium Khronos Group allowing the execution of programs without restrictions from manufacturers [6]. Today, these two APIs remain the most popular. Since the CUDA API is specially developed for NVIDIA GPUs, it is optimized for this hardware and therefore allows more efficient implementations. In addition, NVIDIA provides quality technical support to the developer. On the other hand, OpenCL being an open standard (supported by AMD), has the major advantage of its portability [7].

A2.3. GPU programming with CUDA

In this section, we will present the basics of C CUDA programming, but there is also the possibility of programming in Fortran CUDA or Python CUDA. First, before any programming, you need the programming tools that are Visual Studio C++ for Windows and Eclipse for Linux. To this, we must add the NVIDIA CUDA Toolkit which contains the compiler and the analysis tools of the codes executed. With the development platform installed, we can address the programming part.

It should be noted that the classic C and C CUDA codes are undifferentiated. It is the compiler that distinguishes the keywords and compiles respectively the sections of code to be executed on the CPU and those for the GPU. Indeed, in the CPU system (called host) and GPU (called device) connected by the PCIe bus (Fig. A2.4.), many keywords allow the compiler to perform its task. It essentially includes :

• _global_ which tells the compiler that this function will be called by the CPU and executed on the GPU,

• _device_ for a function that will be called by the GPU to be executed on the GPU

• _host_ that we can omit which is for the CPU functions executed by the CPU. When we want to write a program that will run a function on a GPU (called Kernel and identified by _global_) its structure must include at least the following steps:

• Allocation of memory resources on the GPU.

• Initialization of some GPU variables.

• Transfer all necessary variables from the CPU to the GPU.

• Configuring threads for execution on the GPU.

• Call of calculation function for GPU.

• Transfer GPU results to CPU.

• Release of GPU memory resources and the GPU itself.

To make such an implementation, the first observation that emerges is that the functions of allocation, release, initialization and transfer are all inspired by the classic C, we have

cudaMalloc(), cudaFree(), cudaMemset() and variable transfer with cudaMemcpy(…, direction)

specifying the cudaMemcpyHostToDevice direction from CPU to GPU and cudaMemcpyDeviceToHost direction from GPU to CPU. If we continue with the implementation, we have cudaDeviceReset() which is used to release the GPU. The second observation is that the most important steps in the implementation are the threads organization and memory organization, because the performance of the function execution depends on it. The threads and memory organization are described below.

| Appendix xlv v Threads organization

As you probably know, the GPUs are SIMD (Single instruction on Multiple Data), which means that the code executed by the Kernel that we have created will have to be executed for all the data processed. This is where the organization of threads comes in. Indeed, this organization must take into account both the structure of the data, the operations carried out by the function and also the architecture of the GPU. In the basic concept, we must therefore see the set of threads called for the Kernel which runs on the GPU as a grid (Grid) which consists of several blocks (Block) each containing several Threads (Fig. A2.5.). Note that all threads sharing the same memory space but only threads in the same block can communicate with each other. For the organization of blocks in the grid and threads in blocks, it can be done in 1D, 2D or 3D. Once this organization is defined we can run the Kernel from the CPU with this syntax Namefunction <<< grid, block >>> (variables) followed by cudaDeviceSynchronize() to ensure that all results are ready before starting their transfers to the CPU memory.

Conclusion et Perspectives

Fig. 1 . 1 .Fig. 1 . 2 . 11 Fig. 1 . 3 . 12 Fig. 1 . 4 . 14 Fig. 1 . 5 . 16 Fig. 1 . 6 . 18 Fig. 1 . 7 . 19 Fig. 1 . 8 .Fig. 2 . 1 . 30 Fig. 2 . 2 . 31 Fig. 2 . 3 . 32 Fig. 2 . 4 . 32 Fig. 2 . 5 . 33 Fig. 2 . 6 . 34 Fig. 2 . 7 . 36 Fig. 2 . 8 . 39 Fig. 2 . 9 .Fig. 2 . 10 .Fig. 2 . 11 .Fig. 2 . 12 .Fig. 2 . 13 .Fig. 2 . 14 .Fig. 2 . 15 .Fig. 2 . 16 .Fig. 3 . 1 .Fig. 3 . 2 .Fig. 3 . 3 .Fig. 4 . 1 .Fig. 4 . 2 .Fig. 4 . 3 .Fig. 4 .

 111211131214141516161817191821302231233224322533263427362839292102112122132142152163132334142434 Fig. 1.1. AIRO mobile intraoperative CT scanner and a scan volume of 50 x 100cm of an entire spine ………………………………………………………………………………...………...…10 Fig. 1.2. A set of four MR sequences provides complementary information for the diagnosis of a brain tumor. (a) Proton density-weighted sequence, (b) T1-weighted native sequence, (c) contrast enhanced T1-weighted sequence, and (d) T2-weighted sequence [7] ……………………….…..11 Fig. 1.3. GE Healthcare Optima™ MR450w GEM -70cm intraoperative MR system and scan of a pediatric case ………………………………………………………………………………..….12 Fig. 1.4. BK ultrasound intraoperative system Flex Focus 800 and a few acquired images ……..14 Fig. 1.5. Biograph Vision PET/CT scanners of Siemens Healthineers and 18 F-Fluroestradiol PET/CT scan of a patient with breast carcinoma …………………………………….…………..16 Fig. 1.6. Universal OLYMPUS EVIS EXERA III imaging platform and esophagus image obtain with White Light Imaging and Narrow Band Imaging ……………………………….………….18 Fig. 1.7. Mauna Kea Technologies Cellvizio flexible microscopic imaging system and confocal laser endomicroscopy image of a healthy colon and of an adenocarcinoma-invaded colon…...….19 Fig. 1.8. Augmedics xvision TM augmented reality guidance system for surgery and an image of implantation during spine procedures ……………………………………………………....……21 Fig. 2.1. Light reflection and refraction at a surface interface ……………….………………….30 Fig. 2.2. a) Absorption in a medium. b) Intensity decay according to Beer-Lambert Law …..….31 Fig. 2.3. Scattering in a medium …………………………………………………….……….…..32 Fig. 2.4. Rayleigh and Mie scattering patterns [4] ……………………………..………………..32 Fig. 2.5. Scattering event ……………………………………………………………….………..33 Fig. 2.6. Principles of fluorescence. Jablonski energy diagram [8] ……………………………..34 Fig. 2.7. Absorption coefficient for different chromophores present in human tissue [10] ….….36 Fig. 2.8. Biological structures of various sizes for photon scattering [11] ………………….…...39 Fig. 2.9. Radiative transfer equation conservation of energy ……………………………………40 Fig. 2.10. Flowchart for the variable step size Monte Carlo technique [15] …………………….43 Fig. 2.11. Continuous wave NIR (CW NIR) [10] ……………………………………………….45 Fig. 2.12. Time domain NIR (TD NIR) [10] …………………………………………………….46 Fig. 2.13. Frequency domain NIR (FD NIR) [10] ………………………………………………47

Fig. 4 . 7 .Fig. 4 . 8 .Fig. 4 . 9 .Fig. 5 . 1 .Fig. 5 . 2 .Fig. 5 . 3 .Fig. 5 . 4 .Fig. 5 . 5 .Fig. 5 . 6 .Fig. 6 . 1 .Fig. 6 . 2 .Fig. 6 . 3 .Fig. 7 . 1 .Fig. 7 . 2 .Fig. 7 . 3 .Fig. 7 . 4 .Fig. 7 . 5 .Fig. 7 . 6 .Fig. 7 . 7 .Fig. A2. 1 .

 474849515253545556616263717273747576771 Fig. 4.7. In vivo measurements: absorption images obtained at four spatial frequencies with (a) SFDI, (b) rectangular filters DC1 and AC1 (SSOP std), (c) filters DC1 and AC3 (SSOP opt1), and (d) filters DC3 and AC3 (SSOP opt2) ………………………………….………………………..77 Fig. 4.8. In vivo measurements: reduced scattering images obtained at four spatial frequencies with (a) SFDI, (b) rectangular filters DC1 and AC1 (SSOP std), (c) filters DC1 and AC3 (SSOP opt1), and (d) filters DC3 and AC3 (SSOP opt2) ………………………………………………………78 Fig. 4.9. In vivo absorption and reducer scattering measurements: (a) standard DC1 and AC1 filter combination and (b) DC3 and AC3 filters combination (video 1, MPEG, 58MB) …….…………79

  Photon incident direction

1. 3 . 1 . 4 .

 314 Advantages and disadvantages ................................................................................ 1.3.1. Systems comparison ..................................................................................................... 1.3.2. Impact for surgeons ...................................................................................................... Conclusion ................................................................................................................

Fig. 1 . 1 .

 11 Fig. 1.1. AIRO mobile intraoperative CT scanner and a scan volume of 50 x 100cm of an entire spine.

( 1 -

 1 Fig. 1.2. A set of four MR sequences provides complementary information for the diagnosis of a brain tumor. (a) Proton density-weighted sequence, (b) T1-weighted native sequence, (c) contrast enhanced T1-weighted sequence, and (d) T2-weighted sequence [7].

Fig. 1 . 3 .

 13 Fig. 1.3. GE Healthcare Optima™ MR450w GEM -70cm intraoperative MR system and scan of a pediatric case.

Fig. 1 . 5 .

 15 Fig. 1.5. Biograph Vision PET/CT scanners of Siemens Healthineers and 18 F-Fluroestradiol PET/CT scan of a patient with breast carcinoma.

Figure 1 . 6 .

 16 shows an example of a universal imaging platform for endoscopic and surgical procedures across multiple healthcare specialties.

Fig. 1 . 6 .Fig. 1 . 7 .

 1617 Fig. 1.6. Universal OLYMPUS EVIS EXERA III imaging platform and esophagus image obtained with White Light Imaging and Narrow Band Imaging

Fig. 1 . 8 .

 18 Fig. 1.8. Augmedics xvision TM augmented reality guidance system for surgery and an image of implantation during spine procedures.

Fig. 2 . 1 .

 21 Fig. 2.1. Light reflection and refraction at a surface interface.

Fig. 2 .

 2 Fig. 2.2. a) Absorption in a medium. b) Intensity decay according to Beer-Lambert Law.

Fig. 2 . 4 .

 24 Fig. 2.4. Rayleigh and Mie scattering patterns [4].

Fig. 2 . 5 .

 25 Fig. 2.5. Scattering event.

Fig. 2 . 6 .

 26 Fig. 2.6. Principles of fluorescence. Jablonski energy diagram[8]

Fig. 2 . 7 .

 27 Fig. 2.7. Absorption coefficient for different chromophores present in human tissue [10].

Fig. 2 . 9 .

 29 Fig. 2.9. Radiative transfer equation conservation of energy.

  only depends on the diffusion angle θUnder all these approximations, the diffusion equation describes the fluence rate as a function of space r r and time t as follows:

rt f rr

  the photon density, i.e. the Green's function of the differential equation -the spatial distribution of the absorption coefficient in the medium -( ) Dr r the spatial distribution of the diffusion coefficient in the medium, which depends both on the diffusion approximation equation can be then deducted depending on boundary conditions (e.g. layers of tissue) and light source configurations. The simplest function, also known as the Green's function, is obtained for a short-pulsed point source in an infinite homogeneous medium[12]:

Fig. 2 .

 2 Fig. 2.10. Flowchart for the variable step size Monte Carlo technique [15].

2 . 11 .

 211 photodiode detector[20,21]. CW NIRS systems are less expensive and have faster sampling rates than time-resolved and frequency domain NIRS. The disadvantage of CW systems is that they cannot fully determine the optical properties of tissue and therefore they can present issues when determining functional or structural parameters. To address this issue using current continuous wave instruments, multiple optodes operating simultaneously are placed around the light emission head. The combination of the multi-distance measurements of optical attenuation enables the calculation of the relative concentrations in the illuminated tissue. This modification is called spatially resolved spectroscopy (SRS).

Fig. 2 .

 2 Fig.2.11. Continuous wave NIR (CW NIR)[10]

Fig. 2 . 2 . 4 . 2 . 4 . 1 .

 224241 Fig.2.13. Frequency domain NIR (FD NIR)[10]

Fig. 2 .

 2 Fig. 2.15. SFDI acquisition and processing workflow [25].

3 Scientific issues and objectives Contents 3 . 1 .

 331 General purpose of the project ............................................................................... 3.2. SFDI/SSOP method limitations .............................................................................. 3.3. Artifact reduction and accuracy improvement ...................................................... 3.4. Real-time multispectral imaging acquisition system design .................................. 3.5. Real-time processing and visualization imaging system ........................................ 3.6. Conclusion ............................................................................................................... | Chapter 3. Scientific issues and objectives 57

Figure 3 . 1 .

 31 , to the number of wavelengths required for the extraction of functional and structural parameters. However, the use of a single image causes a loss of spatial resolution and unwanted artifacts mostly localized on the edges. To improve the accuracy of optical properties measurements with SFDI/SSOP, the measurement of the object profile must be considered. If SFDI requires two sets of measurements to obtain both the profile and the optical properties, a 3D-SSOP version (Fig.3.2.) offers the use of a pattern grid to extract them with a single image[2]. Although operational, this solution is subject to the same image quality degradations due to the use of pattern grid.

Fig. 3 .

 3 Fig. 3.1. a) SFDI/SSOP acquisition and processing workflow. b) In vivo optical properties results obtained with SFDI/SSOP.

Fig. 3 .

 3 Fig. 3.2. 3D-SSOP acquisition and processing workflow

| Chapter 3 .

 3 Fig. 3.3. FPGA vs GPU [6].

  where , ( , ) and , ( , ) are measured on the calibration reference, and , , ( , ), , , ( , ) are modeled based on the known optical properties of the medium. Using the medium diffuse reflectance, a pre-computed lookup table can be used then to recover reduced scattering (µs ' ) and absorption (µa). This process is summarized in Fig. 4.2.(a).

Fig. 4 . 2 .

 42 Fig. 4.2. Acquisition and processing flowchart: (a) SFDI process and (b) SSOP process: from raw images, to demodulated images, to absorption and reduced scattering images. Note the use of a single raw image and a 2D FT for processing SSOP data.

  2 mm -1 ) to reduce the acquisition time and extract the MAC and MDC images by filtering the acquired image directly in the spatial frequency domain [Fig. 4.2.(b)]. A one-dimensional (1D)

Fig. 4 . 3 .

 43 Fig. 4.3. Schematics of the filters used in the study.

  Results quantifying accuracy for the selected filters on N = 10 hands are shown in Fig.4.4., colorcoded for ease of interpretation. Overall, the mean error expressed in percentage value for absorption varies from 25.31 % in the worst case to 8.84 % in the best case. Low mean error percentage value indicates a better quality, and we added a color-bar along the numerical values for easy interpretation. Results for all 338 filters combinations are shown in Appendix 1. More precisely, the best combinations to extract absorption, as a function of spatial frequency, are as follows:

2 ; 2 = 2 = 1 . 3 *Fig. 4 . 4 .

 221344 Fig. 4.4. Mean error percentage value for absorption at each spatial frequency using the selected filters combinations. Results are given as % mean ±% standard deviation. Color-coding is used for

Fig. 4 . 9 .

 49 Fig. 4.9.(a) the optimal DC3 and AC3 filters combination Fig. 4.9.(b).

Fig. 4 . 7 .Fig. 4 . 8 .

 4748 Fig. 4.7. In vivo measurements: absorption images obtained at four spatial frequencies with (a) SFDI, (b) rectangular filters DC1 and AC1 (SSOP std), (c) filters DC1 and AC3 (SSOP opt1), and (d) filters DC3 and AC3 (SSOP opt2).

Fig. 4 . 9 .

 49 Fig. 4.9. In vivo absorption and reducer scattering measurements: (a) standard DC1 and AC1 filter combination and (b) DC3 and AC3 filters combination (video 1, MPEG, 58MB).

Fig 4 . 8

 48 .. Because amplitude modulation and thereby signal-to-noise significantly decreases when increasing spatial frequency, SFDI demodulation are becoming more noisy when SSOP demodulation still give good results. This increase in noise partially explains why the mean error percentage values for reduced scattering are not decreasing when increasing the spatial frequency from 0.3 mm -1 to 0.4 mm -1 as one would expect.

(Fig. 5 . 1 .

 51 Fig. 5.1. Spatio-temporal modulation of light: acquisition and processing. Sources at different wavelengths are modulated in time at specifically chosen frequency. Light from these sources is combined and modulated in space using a spatial light modulator. The scene is acquired using a fast camera. The acquired images are then processed in the temporal frequency domain to separate the contribution from the different light sources and processed in the spatial frequency domain to extract the optical properties at all wavelengths.

5 % ± 1. 1 %

 51 in reduced scattering. Actual images of four phantoms are shown in Fig. 5.3. for both wavelengths, in absorption and reduced scattering. Results from standard SFDI are shown, as well as spatio-temporal modulation of light, a percentage error map comparing both methods and a line profile. Overall good visual agreement is confirmed in these images.

Fig. 5 . 2 .

 52 Fig. 5.2. Results from phantoms experiments. Absorption and reduced scattering are shown (mean and standard deviation) over a wide range of optical properties. Results are shown for both wavelengths using standard SFDI as a reference and spatio-temporal modulation of light.

Fig. 5 . 3 .

 53 Fig. 5.3. Actual images of tissue mimicking phantoms. Absorption and reduced scattering images of four phantoms are shown over a range of optical properties. Results are shown for both wavelengths using standard SFDI as a reference and spatio-temporal modulation of light. In addition, relative error maps are shown and line profiles compared (dash blue: standard SFDI, solid red: Spatio-temporal modulation of light).

Finally, results from

  the arm cuff occlusion experiment are shown in Fig 5.5.. Significant decrease in oxygen saturation is visible during the occlusion, with values going from 78 % down to 49 %, following by a hyperemia phase after the release of the occlusion with a slow decrease back to the resting state. A movie is shown in Fig. 5.6.. Note that the movie frame rate has been decreased from 100 frames per second to 25 fps to reduce the size of the file.

Fig. 5 . 4 .

 54 Fig. 5.4. Actual images of a still hand. Absorption and reduced scattering images of a still hand are shown over a range of optical properties. Results are shown for both wavelengths using standard SFDI as a reference and spatio-temporal modulation of light. In addition, relative error maps are shown and line profiles compared (dash blue: standard SFDI, solid red: Spatio-temporal modulation of light).

Fig. 5 . 5 .

 55 Fig. 5.5. Arm cuff occlusion results. An arm cuff occlusion was performed for 4 minutes following 3 minutes of resting state and followed by 3 minutes of recovery. Images were acquired using spatio-temporal modulation of light at 100 frames per second and processed for oxygenation. Top images show oxygen saturation at 3 time points during the baseline resting measurement, the occlusion, and after the release. The saturation of a point in the hand (blue square) is shown in the plot below. Time points of the above images are indicated by arrows.

6. 3 . 6 . 4 . 6 . 5 .

 36465 Results.....................................................................................................................118 Discussion ...............................................................................................................121 Conclusion ..............................................................................................................122 | Chapter 6. Real-time optical properties and oxygenation imaging using GPGPU 106

  Fi. As shown by the equation (1), the value of the chosen frequencies Fi (avoiding harmonics or power-line frequencies) are exact frequencies in the discrete Fourier domain (to avoid crosstalk) and depend on the camera framerate fps, the number of images used for the demodulation process N, and the frequency positions in the Fourier domain 1

Fig. 6 . 1 .

 61 Fig. 6.1. Schematics of the SSOP imaging system using spatio-temporal modulation of light associated with the processing workflow.

  A workstation with the following characteristics was used for the processing: Intel i7-7800x 3.5GHz central processing unit, 16 GB of RAM, and an NVIDIA GeForce GTX 1080TI graphics processing unit (3584 Cores CUDA, 11 GB of RAM).

. 2 .

 2 . The first task of the implementation executes the time demodulation algorithm by computing the amplitude modulation over time at each pixel with the DFT Eq. (6.3) and (6.4). For this purpose, the cosine and sine DFT coefficient at k = 2 and k = 3 are pre-computed and stored in memory. The acquired images are read in the acquisition order and stored in a ring buffer array of N+1 images allowing to have access to the oldest image of the stack and to the new one.

Fig. 6 . 2 .

 62 Fig. 6.2. General processing diagram flow

  first step is the computation of the temporal DFT using pre-computed DFT cosine and sine coefficients stored in MATLAB workspace. The extraction of MDC and MAC images in spatial frequency domain is then performed using the SSOP demodulation principle with the MATLAB | Chapter 6. Real-time optical properties and oxygenation imaging using GPGPU 115 functions fft2, ifft2 and abs. The optical properties a m and ' s m estimation for the MATLAB CPU reference implementation is performed with either griddata triangulation-based 2D interpolation MATLAB function or bilinear interpolation MATLAB function interp2 on the LUT. Note that the processing using the griddata function was carried out in double precision since this function was not designed to work in single precision. The other processing parts are custom-made code as described by their equations.

  neighbor value in the LUT of the measured ,, dA Cx R , and m represents either a m or ' s m extracted from the LUT. In the last task, a custom-made function is used to compute the saturated oxygenation level.

Fig. 6 . 3 .

 63 Fig. 6.3. Results obtained with the MATLAB CPU reference code and with GPU C CUDA code. Absorption and reduced scattering map recovered at 665 nm (left columns) and 860 nm (center columns) and saturated oxygenation level (right column) obtained with reference MATLAB CPU processing (first row) and the custom-made GPU C CUDA processing (second row). Percentage

Fig. 7 .

 7 Fig. 7.1. (A) Schematics of a SFDI/SSOP imaging system. (B) SFDI/SSOP processing workflow.

Fig. 7 . 2 .

 72 Fig. 7.2. Detailed architecture of the proposed deep learning network architecture

  The training starts with a learning rate of 0.001 and employs a learning rate reduction (factor of 1.11) if there is no loss improvement after 15 epochs. The model is trained for 900 epochs using NVIDIA GTX 1080Ti computation power in less than 2 hours. The training dataset consists of a total of 100 images that are 1024 x 1024 pixels in size composed of a combination of 7-phases SFDI acquisitions (for both accuracy and visual quality) of hands (N = 60 from 10 different hands placed in various orientations), silicon-based tissue mimicking phantoms (N = 20 from 18 different phantoms in various orientations) and ex vivo pig organs (N = 20 from 6 different organs). A 5-fold crossvalidated was performed on the network and no signs of overfitting was observed.

Fig. 7 . 3 .Fig. 7 . 4 .Fig. 7 . 6 .

 737476 Fig. 7.3. Profile measurements: 3D profile maps obtained at four spatial frequencies with SFDI method and three SSOP processing methods along with their corresponding error maps. From left to right: Raw input data, SFDI reference, SSOP Filtering, SSOP Deep Learning Convolution (DL conv) 7x1, SSOP Deep Learning Convolution 7x7, SSOP Filtering percentage error, SSOP Deep Learning Convolution 7x1 percentage error, SSOP Deep Learning Convolution 7x7 percentage error.

Fig. 7 . 7 .

 77 Fig. 7.7. Single-frame excerpt from video recording of a moving hand (visualization 1): Raw input data, 3D profile, absorption coefficient and reduced scattering coefficient. From top to bottom: SSOP Filtering, SSOP Deep Learning Convolution 7x1, SSOP Deep Learning Convolution 7x7.

  3 mm -1 and 0.4 mm -1 are still subject to some imperfections, which likely affect the network training and subsequent reconstruction upon inference (see Fig3-5). These artifacts could explain, in part, why the error values for absorption and reduced scattering are not decreasing when increasing the spatial frequency from 0.3 mm -1 to 0.4 mm -1 .

. 1 *

 1

3 *

 3 d: field of view width in mm xlii • 4 schedulder warp that schedule the batch of threads that were assigned to them by the Gigathread Engine.

Fig. A2. 3 .

 3 Fig. A2.3. SMM internal structure [5].

Fig. A2. 4 .

 4 Fig. A2.4. CPU-GPU system [8].

Fig. A2. 5 .

 5 Fig. A2.5. Threads organization [8].

2. Chapitre 1 :| Résumé étendu en français lviii 3 . Chapitre 2 :| Résumé étendu en français lix 4 . Chapitre 3 : 5 . 5 % 6 .| Résumé étendu en français lxi 7 . Chapitre 6 :. 8 .

 13243556768 Chirurgie guidée par l'imagerie La conception d'un système d'imagerie médicale doit être justifiée et ses avantages par rapport aux systèmes existants doivent être démontrés. En effet, le chemin menant à la réalisation du produit final est long et parsemé d'embûches. C'est dans ce cadre que dans le Chapitre 1, nous avons fait le tour d'horizon des modalités d'imagerie pour la chirurgie guidée par l'image pour situer la contribution que nous faisons dans le domaine de l'imagerie optique. Nous avons présenté brièvement les modes de fonctionnement, les avantages, les limites et aussi les champs d'application des modalités d'imagerie pour la chirurgie les plus utilisées que sont : l'imagerie optique, l'imagerie ultrason, l'imagerie à résonnance magnétique, l'imagerie nucléaire, l'imagerie à rayon X et l'imagerie basée sur la réalité augmentée. Les différentes caractéristiques des modalités d'imagerie présentées ont été résumées dans le tableau 1.1. Les techniques d'imagerie optique sont clairement celles qui présentent les caractéristiques les plus avantageuses et de la polyvalence. Lorsque pour les autres, aucune modalité n'est supérieure dans tous les domaines, et en fait, les modalités sont souvent complémentaires. Cela explique le développement de système d'imagerie multimodalité (hybride) comme SPECT/CT, PET/CT, PET/MRI, US/CT ou US/MRI ou encore par la conception de nouveaux systèmes basés sur la réalité augmentée. De notre point de vue, dans plusieurs années, les systèmes basés sur la réalité augmentée auront le meilleur de toutes les modalités d'imagerie en intégrant des images préopératoires et intra-opératoires dans un système de visualisation en temps réel immersif. Bien que pour l'heure de nombreux verrous technologiques restent à lever. Spectroscopie proche infrarouge L'imagerie optique apporte à faible coût beaucoup d'avantages pour la chirurgie guidée par image. Le Chapitre 2 présente les concepts principaux de l'imagerie optique utilisée dans le domaine de la spectroscopie proche infrarouge. Une description des phénomènes d'interaction entre les tissus et la lumière a été donnée pour l'absorption, la diffusion et la fluorescence. L'absorption a été définie comme la conséquence de la conversion partielle de l'énergie lumineuse électromagnétique en une autre forme d'énergie. La diffusion quant à elle est un phénomène qui se produit lorsque les photons rencontrent une particule avec un indice de réfraction différent de celui des entités biologiques, produisant une déviation de la direction incidente du photon. Enfin, la fluorescence est un phénomène au cours duquel une molécule émet de la lumière lorsqu'elle revient à son niveau fondamental après avoir été excitée par une source de lumière. Les lois qui régissent les relations entre les propriétés optiques (absorption et diffusion) et les paramètres fonctionnels et structurels des tissus (taux d'oxygénation, taux de lipides, taux de mélanine et d'eau) ont été présentées. Les différents modèles qui peuvent servir à modéliser le transport de l'énergie dans un milieu absorbant et diffusant ont été revus (transfert radiatif, approximation de diffusion et model Monte Carlo). Nous avons également décrit les trois principales techniques utilisées pour mesurer les propriétés optiques des tissus biologiques dans le proche infrarouge. La méthode de l'illumination continue qui a une réponse très rapide mais enregistre seulement les changements relatifs, et il n'est donc pas possible de faire des mesures absolues en utilisant cette technique. La méthode du temps de vol de photon qui nécessite un traitement approfondi des données, mais fournit des mesures plus précises. Elle permet d'explorer différentes informations fournies par la mesure. La troisième approche, qui utilise la technologie de modulation de fréquence ou de phase, et qui a une résolution plus faible que celle de la méthode du temps de vol. Et pour terminer une description des techniques d'imagerie par la fluorescence et par la réflectance diffuse, largement utilisées dans le proche infrarouge, a été donné avec des exemples de système et d'application. Problématiques et objectifs Dans ce Chapitre, nous avons présenté les limites à l'utilisation en chirurgie guidée par l'imagerie de la méthode d'imagerie par la réflectance diffuse SFDI et de sa variante à acquisition rapide SSOP. Ensuite nous avons présenté les pistes de solutions que nous avons étudiées. En résumé, pour ce qui est des limitations, le traitement SSOP engendre une perte de résolution spatiale des images et crée des artefacts indésirables essentiellement localisés sur les bords. De plus, pour améliorer la précision des mesures de propriétés optiques avec SFDI/SSOP, la mesure du profil de l'objet doit être prise en compte. Ensuite pour mesurer les paramètres physiologiques et structurels il est nécessaire d'effectuer une acquisition multispectrale, qui se fait généralement soit de manière séquentielle ou avec un système optique complexe. Selon l'approche utilisée, le temps d'acquisition peut être long, le rapport signal-bruit faible, l'acquisition peut être affectée par le mouvement et le coût du système élevé. Et pour finir la méthode nécessite une implémentation temps réel pour être utilisé en chirurgie. Les objectifs de cette thèse sont d'apporter des solutions à ces limitations. Dans les chapitres qui suivent les Chapitre 4 et 7 vont aborder l'aspect amélioration de la qualité des images et de leurs précisions. Dans le Chapitre 5 nous allons présenter la technique de modulation temporelle et spatiale de lumière pour résoudre les problèmes liés à l'acquisition multispectrale. Enfin les Chapitre 6 et 7 vont être dédié à l'implémentation temps réel. Chapitre 4 : Amélioration de la qualité des image SSOP par l'utilisation de filtre anisotropique 2D Dans ce Chapitre, nous avons présenté les moyens d'améliorer la qualité d'image et la précisions des propriétés optiques estimées par la méthode SSOP. Nous avons étudié en particulier deux paramètres : la stratégie de filtration dans le domaine des fréquences spatiales, et le choix de la fréquence spatiale. L'effet de ces paramètres sur les propriétés optiques estimées et la qualité d'image par rapport à la méthode SFDI standard ont été évalués. Pour cela, nous avons étudié des filtres 2D anisotropiques construits à partir de fenêtres autres que la fenêtre rectangulaire telle que Sine, Hann, Blackman et Gaussien. Cette étude a été réalisée pour des fréquences spatiales allant de 0,1 à 0,4 mm -1 directement sur N = 10 mains. Les images de ces 10 mains ont été acquises avec la méthode SFDI, puis utilisées avec la méthode SSOP pour extraire les propriétés optiques en | Résumé étendu en français lx utilisant les différents filtres conçus. L'exactitude des résultats et la qualité d'image des propriétés optiques ont été quantifiées en comparaison avec la méthode des acquisitions multiples SFDI. Dans l'ensemble, les équations permettant de concevoir les filtres anisotropes DC et AC ont été fournies et les meilleures combinaisons de fenêtres de filtrage à chaque fréquence spatiale ont été présentées. L'une des meilleures combinaisons, basée sur les fenêtres Sine DC et Blackman AC, permet une amélioration significative de la qualité des images tout en préservant la précision des mesures de propriétés optiques, en particulier à des fréquences spatiales élevées (typiquement, 0,4 mm-1) avec des erreurs moyennes de prédiction des propriétés optiques restant inférieures à 9,en absorption et 8,5 % en diffusion réduite contre respectivement 11 % et 9,5 % pour le filtrage classique. Cette étude a également démontré la dépendance de la qualité et de la précision de l'image avec la fréquence spatiale. Chapitre 5 : Mesure quantitative de l'oxygénation en utilisant la modulation spatiotemporelle de la lumière Dans le Chapitre 5, nous avons présenté la nouvelle méthode d'acquisition multispectrale utilisant la modulation spatiale et temporelle de la lumière pour fournir en temps réel des images interprétables par les chirurgiens. La méthode proposée consiste à utiliser des sources modulées sinusoïdalement dans le temps à l'aide de générateurs de signaux pour coder les longueurs d'onde et à projeter simultanément le motif structuré pour l'extraction des propriétés optiques par la méthode SSOP. Les images de la scène sont acquises à l'aide d'une caméra rapide et sont démodulées au moyen de Transformé de Fourier Discrète pour séparer les longueurs d'onde. Les images résultantes sont ensuite traitées par SSOP pour extraire les propriétés optiques de chaque longueur d'onde et estimer les paramètres physiologiques et structurels. La méthode a été validée sur un système construit avec deux longueurs d'onde (665 nm et 860 nm), et donc capable de fournir des propriétés optiques à ces deux longueurs d'onde, la concentration d'oxyhémoglobine, de désoxyhémoglobine et le taux d'oxygénation. Une validation a été effectuée en comparant les résultats de la méthode avec ceux obtenus pour une acquisition séquentielle SFDI au cours de deux expériences, l'une avec des fantômes et l'autre avec une main in vivo. Les résultats ont montré un taux d'erreur maximal de 4,2 %. Enfin, nous avons réalisé une expérience d'occlusion sur une main qui a confirmé le potentiel clinique du système pour la détection d'ischémie. Mesure en temps réel des propriétés optique et de l'oxygénation au moyen de processeur graphique L'objectif de ce Chapitre a été d'explorer les possibilités offertes par le traitement sur processeur graphique (GPU) des données SFDI modulées dans le temps (pour le codage des longueurs d'onde) et dans l'espace (pour l'extraction des propriétés optiques) pour l'imagerie optique diffuse multispectrale quantitative rapide. Plus précisément, nous avons proposé des implémentations GPU pour effectuer le traitement en temps réel de la démodulation temporelle à la mesure du taux d'oxygénation sur le système conçu au Chapitre 5. Les équations de chaque fonction à chaque étape du traitement ont été fournies. Nous avons aussi décrit les implémentations, à la fois CPU utilisant MATLAB, GPU utilisant MATLAB et GPU utilisant l'implémentation C CUDA que nous avons faites. Enfin, nous avons testé ces implémentations sur des données SSOP obtenues d'une main humaine et avons comparé leurs différences en temps de traitement et en termes de précision. L'implémentation C CUDA comparé aux implémentations conventionnelles MATLAB CPU/GPU a montré un taux d'erreur inférieur à 0,1% pour l'extraction des propriétés optiques. Enfin, le résultat le plus important obtenu dans ce Chapitre a été la capacité de mesurer en seulement 1,6ms le niveau d'oxygénation sur des images de 1 mégapixel garantissant une visualisation en 25 images par seconde sur le système de démonstration présent dans notre laboratoire. Chapitre 7 : Mesure SSOP de qualité supérieur en temps réel avec correction de profil 3D au moyen de réseau de neurones et de l'utilisation de processeur graphique Dans le Chapitre 4 nous avons présenté l'utilisation de filtre anisotropique pour améliorer la précision des résultats SSOP qui présentaient cependant encore des artefacts sur les bords. Dans le Chapitre 7 nous proposons une nouvelle approche qui permet d'avoir à la fois des images de plus grande qualité mais aussi d'avoir un profil 3D mieux défini nous permettant de fournir des mesures de propriétés optiques corrigées en fonction de leur profil 3D. En particulier, nous utilisons un CNN pour extraire simultanément les images MDC, MAC et le profil 3D de l'échantillon à partir d'une seule entrée d'image SSOP. Pour accomplir cette tâche deux variantes de réseau basées sur la structure U-NET ont été proposées, un réseau qui fonctionne ligne par ligne | Résumé étendu en français lxii (Line U-NET) et un réseau qui traite directement l'image (Image U-NET). Les deux réseaux proposés se composent d'un chemin de contraction de 4 étapes (côté gauche) et d'un chemin expansif de 4 étapes (côté droit). Le chemin expansif a une bifurcation qui conduit à 2 sorties. Au total, les réseaux U-NET ont respectivement 3064 poids pour le « Line U-NET » et 18886 poids pour le « Image U-NET » pour 17 couches convolutionnelles. Le réseau a été mis en oeuvre avec Tensorflow en utilisant Keras et est compilé avec l'erreur quadratique moyenne (MSE) comme fonction de perte et Adam comme optimiseur. Les réseaux ont été entrainés avec un jeu de données acquise avec la méthode SFDI 7 phases composé d'images de mains (femme, homme, caucasien, africain) et de fantômes prisent dans différentes positions. Enfin, afin de garantir le traitement est la visualisation en temps réel des résultats, toutes les implémentations ont été faites sur GPU en C CUDA. En particulièrement pour celle du CNN nous avons utilisé la bibliothèque cuDNN de NVIDIA. Nous avons validé notre implémentation en termes de précision en mesurant les propriétés optiques corrigées sur un ensemble d'images de mains et d'organes ex vivo. Les résultats présentent une qualité d'image similaire aux images SFDI et des taux d'erreurs de moins de 10 % en absorption, en diffusion réduite et pour le profil 3D sur les deux types de réseaux. Nous avons également évalué le temps d'inférence de nos implémentations SSOP qui donne au total 13.81ms avec le CNN « Line U-NET » et 18ms avec le CNN « Image U-NET ». Cette nouvelle implémentation pour la méthode SSOP jette les bases de la prochaine génération d'imagerie optique diffuse, quantitative, multispectrale et temps réel pour le guidage chirurgical et les applications médicales.

Durant ces 3

 3 ans de doctorat, j'ai travaillé sur le développement d'un instrument d'imagerie pour le guidage chirurgical utilisant l'imagerie optique. Dans ce manuscrit ont été présentées les innovations qui nous ont permis de progresser dans le projet. Pour rappel notre objectif était de développer un dispositif chirurgical intra-opératoire guidé par image basé sur la méthode d'imagerie par réflectance diffuse SFDI qui permettra de fournir des informations physiologiques et structurelles essentielles aux chirurgiens telles que le niveau d'oxygénation, niveau de lipides, niveau de mélanine, niveau d'eau, pouvoir de diffusion et amplitude de diffusion. Pour atteindre notre objectif, nous avons combinons la méthode d'acquisition rapide de SFDI nommée SSOP | Résumé étendu en français lxiii avec la modulation temporelle de la lumière pour discriminer les propriétés optiques des tissus acquis à plusieurs longueurs d'onde simultanément. Pour rendre cela possible, nous avons dû résoudre le problème de dégradation et de précision de l'image SSOP, concevoir un système d'acquisition d'imagerie multispectral à champ large en temps réel et développer des techniques de traitement en temps réel innovantes pour garantir le traitement en temps réel et la visualisation. Pour l'amélioration de la précision des résultats SSOP et de leurs rendus visuels. Nous avons développé dans un premier temps des filtres anisotropiques en remplacement des filtres rectangulaires habituelle qui ont permis de limiter les erreurs uniquement sur les bords de l'échantillon mesuré. Dans un deuxième temps nous sommes allés plus loin un proposant l'utilisation de CNN pour s'affranchir du défaut sur les bords et d'augmenter le rendu visuel des résultats. Cette amélioration de la qualité des images nous a permis d'aller jusqu'à l'extraction d'un profil 3D bien défini, que nous avons utilisé pour améliorer la précision des résultats SSOP. En plus de l'amélioration de la méthode SSOP nous avons proposé une technique d'acquisition multispectrale utilisant la modulation spatiotemporelle, les études effectuées ont démontré que l'on peut acquérir plusieurs longueurs d'onde simultanément en s'affranchissant des limites de la complexité des méthodes conventionnelles. Pour finir, après avoir développé un système d'acquisition multispectrale rapide, compensé la perte de qualité d'image due à cette dernière, nous avons abordé la contrainte temps primordial en chirurgie en développant des algorithmes innovants et en les implémentant sur processeur graphique en langage bas niveau pour tirer le maximum de puissance des GPU. A la fin de ma thèse nous avons pu développer le premier système de mesure d'oxygénation prenant en compte le profil 3D des échantillons, étant large champ (17.5 x 17.5 cm 2 ), à deux longueurs d'onde (665 nm et 860 nm) et fonctionnant à la cadence de 25 images par secondes. Limitations Le travail que nous avons fait comporte des limites comme toute recherche. La première est la qualité de la référence utilisée pour toutes nos études. En effet, si la méthode SFDI est théoriquement plus précise que le SSOP, elle a ses limites qui sont : la nécessité d'un multi calibrage pour plus de précision en raison de la variété des propriétés optiques présentes sur un échantillon, la baisse du rapport signal à bruit à mesure que la fréquence spatiale augmente, ce qui | Résumé étendu en français lxiv a certainement mis un biais certain dans nos comparaisons aux chapitres 4 et 7. Les deuxièmes limites à la qualité de notre travail sont celles liées à la méthode d'acquisition à modulation spatiotemporelle. En fait, dans nos résultats, nous avons des images obtenues pour la longueur d'onde 860 nm qui ont un rapport signal à bruit plus faible que à 665 nm due à l'efficacité quantique plus faible de la caméra à cette longueur d'onde. Cette plus faible efficacité à 860 nm influence la qualité des propriétés optiques mesurées à 860 nm et de faite la mesure du taux d'oxygénation. L'autre limite aux méthodes développées est liée au nombre maximum de longueur d'onde réglable simultanément qui est limité par la plage dynamique de la caméra. De plus, l'utilisation de fenêtres glissantes pour la démodulation temporelle se fait au prix d'un court délai entre la scène acquise et les informations affichées. Troisièmement sur la nouvelle méthode SSOP basée sur l'apprentissage profond, on peut noter la nécessité de faire des acquisitions plus variées notamment sur des échantillons in vivo pour améliorer la précision des réseaux CNN que nous avons créés. Enfin, si la solution à apprentissage profond pour SSOP est un atout, elle nécessite un temps de traitement 16 fois plus long que la méthode par filtrage ce qui réduit le nombre de longueurs d'onde que l'on peut mesurer simultanément avant de faire évoluer notre architecture GPU. Perspectives Le travail que nous avons fait n'est que la base de ce qui reste à faire. En effet, notre objectif pour ce projet est de fournir au chirurgien un produit fini approuvé sous la forme d'un chariot qui peut facilement être installé dans une salle d'opération et qui fournira en temps réel sur un écran d'affichage au moins 4 paramètres physiologiques (niveau d'oxygénation, niveau de lipides, niveau de mélanine, niveau d'eau) et 2 paramètres structurels (puissance de diffusion, amplitude de diffusion). Comme le montrent nos résultats, il nous manque 5 paramètres. Pour atteindre cet objectif, nous voyons comme premier travaille complémentaire, la nécessité d'augmenter le nombre de longueurs d'onde du système à au moins 6. Cette augmentation du nombre de longueurs d'onde signifie : (1) la conception des sources et d'un système de contrôle pour la modulation temporelle, (2) améliorer le système d'acquisition en choisissant une caméra pouvant atteindre au moins 500 FPS pour assurer l'acquisition en temps réel, (3) augmenter la puissance de calcul du GPU et améliorer certaine de nos implémentations, (4) faire un système de visualisation agréable

Table 1 . 1 .

 11 Comparison of Various Medical Imaging Modalities [38] ………………..….………22

Table 4 .1.

 4

Table 6 . 2 .

 62 Processing time (in milliseconds) from wavelength demodulation to saturated oxygenation level measurement for 512×512 images ……………………………….………….119

Table 6 . 3 .

 63 Processing time (in milliseconds) from wavelength demodulation to saturated oxygenation level measurement for 1024×1024 images …………………………….………….120

Table 7 . 1 .

 71 Structural similarity index of the 3D profile map and the optical properties map for all the validation dataset …………………………………………………...………………………144

Table 7 . 2 .

 72 Mean absolute percentage error (%) of the 3D profile map and the optical properties map for all the validation dataset ………………………………………….……………………145

Table 7 . 3 .

 73 GPU processing time in milliseconds for the three different SSOP processing method ……………………………………………………………………………………..……………146

Table A1 . 1 .

 A11 Bandpass DC filters equations with their cutoff frequency and size ……………..xxvi

Table A1 . 2 .

 A12 Bandpass and highpass AC filters equations with their cutoff frequency and size .xxvii

	Table

1 Image-guided surgery Contents 1

.1. Definition and objectives ........................................................................................... 1.2. Conventional main image-guided system field .........................................................

1.2.1. X-ray imaging and Computed tomography ...................................................................... 1.2.2. Magnetic resonance imaging ......................................................................................... 1.2.3. Ultrasound imaging ....................................................................................................... 1.2.4. Imaging in nuclear medicine .......................................................................................... 1.2.5. Optical imaging ............................................................................................................. 1.2.6. Augmented reality ........................................................................................................

1.7.). Properties

	| Chapter 1. Image-guided surgery
	imaging techniques are involved in intraoperative surgery, such as confocal microscopy, two
	photon microscopy, optical coherence tomography (OCT) among others. For example, OCT is
	employed in ophthalmology and optometry where it can be used to obtain detailed images from
	within the retina or in interventional cardiology to help diagnose coronary artery disease.
	such as coloration,
	texture, lesion borders, and size can therefore be used to form an indicative score that may prompt
	biopsy, treatment, or watch-and-wait paths. Microscopic detection capabilities can be enhanced
	through the use of contrast agents such as methylene blue or toluidine blue which are sprayed on
	or applied to the tissue. Given the constantly developing nature of the field, several microscopic
	20

1.3. Advantages and disadvantages 1.3.1. Systems comparison

	systems also present a few challenges. From a technical point of view, simple hybrid/AR systems
	require to perform multi-modality acquisition, accurate registration and segmentation of the
	measured parameters as well as to provide a useful visualization interface. These tasks require
	further engineering and computer science development. On the other hand, the final cost of this
	type of system will of course be higher than a simple one, even if the global healthcare cost in the
	long run will be lower.
	Characteristics	X-
	The different characteristics of the imaging modalities presented in the previous sections are
	summarized in Table 1.1.. Optical imaging techniques clearly have the most advantageous
	characteristics and versatility. When for the others no single modality is superior in all areas, in
	fact, the modalities are often complementary. That explains the development of multi-modality
	(hybrid) imaging scanners such as SPECT/CT, PET/CT, PET/MRI, US/CT or US/MRI and the
	current design of new systems using Augmented Reality capabilities. Eventually, medical
	Augmented Reality systems will have the most promising imaging modalities by integrating
	preoperative and intraoperative images in an immersive real-time visualization system. Such AR

ray imaging CT imaging MR imaging US imaging Nuclear imaging Optical imaging

	Spatial resolution	Excellent Excellent	Good	Good	Poor	Mixed
	Depth penetration	Excellent Excellent	Excellent	Good	Excellent	Poor
	Soft tissue contrast	Poor	Poor	Excellent	Good	Excellent	Excellent
	Functional parameters	No	Yes	Yes	No	Yes	Yes

  )

	f	,, r z t r	=	4 p	c D r ct r	3/2	e	-	m	a	r ct r	22 22 00 44 zz r zz r Drc t Drc t ee aeö aeö -+ ++ ç÷ ç÷ --ç÷ ç÷ èø èø éù êú ×-êú ëû êú rr	(2.18)

Table 4 .

 4 1. Filter equations and parameters. For each filter presented, the filter equation is given as well as the definition of the cutoff frequency fc as a function of spatial frequency fx. Note that the filters are also scaled according to the field of view width d.

  100

	Error	NM nm --== =ǻ å N M	Result	TEST nm Result Result REFERENCE REFERENCE nm nm -	(6.11)
	where N = M are the image pixel size (1024 or 512) and Result represents either ( a m , ' s m or

2

StO

Table 6 .

 6 1. Mean percentage error in absorption and reduced scattering at 665 nm and 860 nm compared to the reference MATLAB CPU implementation.

Table 6 .

 6 1. presents a detailed description of the percent error score of the extracted maps compared to the reference MATLAB CPU processing. Only the precision results obtained with 1024 x 1024 images are presented in the table as they are substantially equal to those obtained with 512 x 512

						Percentage of error (in %)
		LUT linear op					LUT linear rd
	Image 1024 x 1024	MATLAB CPU		MATLAB CPU			MATLAB GPU			C CUDA GPU
		(Reference)
		Nonlinear interpolation (griddata *)	Linear interpolation (interp2)	Low density rounding	High density rounding	Linear interpolation (interp2)	Low density rounding	High density rounding	Linear interpolation (interp2)	Low density rounding	High density rounding
	Absorption
	at	-	0.01 %	0.38 %	0.04 %	0.01 %	0.38 %	0.04 %	0.02 %	0.38 %	0.04 %
	= 665 nm
	Absorption
	at	-	0.02 %	0.65 %	0.08 %	0.02 %	0.65 %	0.08 %	0.08 %	0.67 %	0.12 %
	= 860 nm
	Reduced
	scattering at	-	0.00 %	0.49 %	0.05 %	0.00 %	0.49 %	0.05 %	0.02 %	0.48 %	0.04 %
	= 665 nm
	Reduced
	scattering at	-	0.02 %	0.88 %	0.11 %	0.02 %	0.88 %	0.11 %	0.08 %	0.93 %	0.16 %
	= 860 nm

* Double precision only

Table 6 .

 6

3

. Processing time (in milliseconds) from wavelength demodulation to saturated oxygenation level measurement for 1024×1024 images.

Table 7 .

 7 ). Moreover, the SSOP DL Conv 7x7 method mimics SFDI image edges with near-perfection in stark contrast with SSOP DL Conv 7x1, which illustrates slightly degraded results. 1. Structural similarity index of the 3D profile map and the optical properties map for all the validation dataset.

			Structural similarity index (SSIM)
	Full image	Spatial	SSOP	SSOP	SSOP
		frequency	Filtering	DL Conv 7x1	DL Conv 7x7
		0.1 mm -1	0.44 ±0.16	0.49 ±0.17	0.50 ±0.16
	3D Profile	0.2 mm -1 0.3 mm -1	0.47 ±0.15 0.47 ±0.15	0.50 ±0.16 0.51 ±0.15	0.51 ±0.16 0.52 ±0.15
		0.4 mm -1	0.48 ±0.15	0.51 ±0.15	0.52 ±0.15
		0.1 mm -1	0.96 ±0.02	0.98 ±0.03	0.98 ±0.01
	Absorption	0.2 mm -1 0.3 mm -1	0.98 ±0.01 0.97 ±0.02	0.98 ±0.03 0.97 ±0.03	0.98 ±0.01 0.98 ±0.02
		0.4 mm -1	0.96 ±0.02	0.96 ±0.03	0.97 ±0.02
		0.1 mm -1	0.75 ±0.08	0.84 ±0.07	0.86 ±0.06
	Reduced	0.2 mm -1	0.75 ±0.07	0.81 ±0.08	0.83 ±0.06
	scattering	0.3 mm -1	0.71 ±0.08	0.76 ±0.09	0.78 ±0.07
		0.4 mm -1	0.67 ±0.10	0.71 ±0.10	0.74 ±0.09
	To quantify this analysis, we have calculated the structural similarity (SSIM) index versus 7-phase
	SFDI images (ground-truth) -results of which are provided in Table 7.1.. The first observation is
	that, compared to the visual quality of the SFDI images at any spatial frequencies, we find a
	decreasing similarity index (hence a decrease in visual quality): (1) SSOP DL Conv 7x7, (2) SSOP
	DL Conv 7x1 and (3) SSOP Filtering. The similarity index varies according to the evaluated

Table 7 .

 7 2. Mean absolute percentage error (%) of the 3D profile map and the optical properties map for all the validation dataset.

  1.6 ms the saturated oxygenation level from images of 1 mega-pixel, guaranteeing a visualization in 25 images per seconds on the demonstrator present in our laboratory.The final achievement of this thesis involved solving some limitations of the SSOP implementation. This work was presented in Chapter 7. The improved SSOP method was still subject to artifacts on the edges, and we had not yet considered the real-time correction of the 3D profile. For this purpose, we were able to design a novel filtering method for SSOP based on a combination of CNN U-NET architecture and GPU processing in order to compensate speed losses. With this new processing strategy, our SSOP implementation is capable of high-quality image reconstruction, with simultaneous acquisition of 3D surface profile. Leading us to a 3D profile corrected and more accurate optical property (OP) extraction in real-time across large fields of view (17.5 x 17.5 cm²). At the validation stage, the methodology has demonstrated results comparable to that of profile-corrected SFDI, with a processing time of 18 ms and errors less than 10% in both profilometry and profile-corrected OPs. This novel processing framework lays the foundation for the next generation of real-time multispectral quantitative diffuse optical imaging for surgical guidance and healthcare applications.

Table A1 .

 A1 1. Bandpass DC filters equations with their cutoff frequency and size

			Cutoff
	DC Filter Name	DC Filter Equation	frequency
			in mm -1

2.5. Conclusion ...............................................................................................................

encouragement, and the confidence he had in hiring me on this project. Thanks to Wilfried who recommended me to Sylvain and who knew how to motivate me and track my progress along the whole thesis.

Time Domain Near-Infrared Spectroscopy

Time domain spectrometers are also called time-of-flight or time-resolved systems. They are based on the injection of a very short laser pulse (~50 ps) into an absorbing tissue and on the recording of the temporal response (time-of-flight) of the photons at some distance from the laser source (Fig. 2.12.). They generally employ a semiconductor or solid-state laser to generate ultra-short pulses. Such a system measures the attenuation by either a synchro-scan streak camera or a timecorrelated single photon counting in which a photon counting detector detects and sorts the received photons by their time of arrival [20,21]. This method uses a mathematical approximation that is based on diffusion theory to allow for the separation of effects due to light absorbance from those due to light scattering. However, such spectrometers are more difficult to adapt to clinical use because of the sampling rate, the instrument size, the instrument weight, and the cost. Extracting the DC modulation could be performed in a different way by calculating the mean of the images obtained at three phases at a high spatial frequency [Eq. (21) in Ref . 14]. However, when performing this operation, dark noise and other constant noise [Ioffset in Eq. (4.3)] are not eliminated and lead to incorrect estimation of diffuse reflectance. In this particular case, as well as when performing SSOP, the contribution of dark noise and other constant noise is eliminated by subtracting a dark image from the acquired image. When using SFDI, the acquisition of multiple images along with Eq. (4.2) allows to eliminate constant noise from the image through image subtraction. We chose the latter method for our reference SFDI data processing.

Following demodulation, the medium diffuse reflectance Rd,AC and Rd,DC are obtained using a calibration reference with known optical properties:

, ( , ) = ( , )

Globally, these results indicate that most accurate reduced scattering properties are extracted with the AC3 filter at low spatial frequency and with the AC2 at high spatial frequency. When combined with DC filters, best results are obtained using either DC2 or DC3 filters. Fig. 4.5. Mean error percentage value for reduced scattering at each spatial frequency using the selected filters combinations. Results are given as % mean ±% standard deviation. Color-coding is used for ease of interpretation (scale on the right).

Quantification of image visual quality

Assessing the visual quality of an image is not straightforward. To overcome this issue, we chose Results can be summarized as follow:

• fx = 0.1 mm -1 : Not measurable.

• fx = 0. as % mean ±% standard deviation. Color-coding is used for ease of interpretation (scale on the right).

Qualitative analysis

To appreciate the improvement in image quality, we provide here an example of the N = 10 hands imaged in the study. Qualitatively, one can first appreciate the image improvement with increasing spatial frequency and filter combination choice. Fine details such as the veins in absorption appear more clearly as suffered from both loss of image resolution and undesirable artifacts due to the single image acquisition strategy [10,11]. Over the years, the SSOP method has been increasingly improved upon to enhance the technique's reconstructive visual quality and capabilities. Though at present current state-of-the-art SSOP implementations, while capable of acceptable visual quality retrieval in real-time, still suffer from edge artifacts, exhibit lower resolution than SFDI and do not include sample profile correction [12]. These limitations can mainly be attributed to a conventional processing workflow.

To overcome these limitations, deep learning (DL) and Graphics Processing Units (GPU) computation provide promising alternatives to conventional processing for diffuse optical imaging. Deep Neural Networks (DNNs), and more particularly Convolution Neural Networks (CNNs), have become ubiquitous in applications such as image classification, image segmentation, image translation, natural language understanding or speech recognition [13][14][15]. Over the last few years, they have also greatly impacted the field of biomedical optics [16][17][18]. Though, if these new methodologies have provided excellent performances, they are typically not amenable to real-time processing due to the increased depth of the networks employed and software implementations. In this regard, the use of GPU-accelerated DL has demonstrated remarkable computational speed increases with regards to model inference in latency sensitives applications [19,20].

In this work, we present a real-time, wide-field and quantitative optical properties imaging implementation of SSOP that includes 3D profile correction and high visual quality reconstruction via GPU-accelerated CNNs. In particular, the workflow detailed herein is used to extract modulated images and the 3D profile of the sample concurrently from a single SSOP image input.

A custom-made GPU implementation of the network, using NVIDIA's cuDNN library, optimizes the network's inference speed -exhibiting speeds comparable to that of the fastest existing approaches with enhanced accuracy. All together, we demonstrate profile corrected optical properties imaging in real-time with visual quality similar to SFDI and errors less than 10% in absorption and reduced scattering.

Upsampling and concatenation

Upsampling and concatenation layers are custom made kernels since these two functions are not currently implemented on cuDNN's library. With each of these bricks, the entire network is then built according to the structure described in the previous section.

GPU configuration

Processing times depending strongly on the GPU configuration. For each custom-made kernel, one thread was created for each pixel, broke down into blocks of 64 threads allowing to have a GPU occupancy of 100%.

Imaging system used for experiments

The instrumental setup (see Fig. 7.1. (A)) was custom-built using a DMD (Vialux, Germany) for the projection of custom patterns, fiber-coupled to a 665 nm laser diode (LDX Optronics, Maryville, Tennessee). The laser diode was mounted using a temperature-controlled mount (TCLDM9, Thorlabs, Newton, New Jersey). Diode intensity was controlled using current controllers (TDC240C, Thorlabs, Newton, New Jersey) and temperature controlled using thermoelectric cooler controllers (TED200C, Thorlabs, Newton, New Jersey). The projection system projects a sine wave pattern over a 175 × 175 mm 2 field of view at 45 cm working distance.

Images of 1024 × 1024 pixels size are acquired using a scientific monochrome 16 bits sCMOS camera (PCO AG, pco.edge 5.5, Kelheim, Germany). Polarizers (PPL05C; Moxtek, Orem, Utah), arranged in a crossed configuration, were used to minimize the contribution from specular reflections at the surface of the sample. The projection system was set at an angle of θ = 4° to allow profilometry measurement [27]. A workstation with the following characteristics was used for controlling the acquisition and the processing: Intel i7-7800x 3.5GHz central processing unit, 16

GB of RAM, and an NVIDIA GeForce GTX 1080Ti GPU (3584 Cores CUDA, 11 GB of RAM). all expressed in mm -1 . A calibration phantom was made with a large size (210 mm × 210 mm × 20 mm) and with absorption of µa = 0.01 mm -1 and reduced scattering µs'=1.1 mm -1 at 665 nm. All phantoms optical properties were referenced using a time domain photon counting system.

Training Set

Validation and Performance Assessment

Our validation dataset consisted of the images of 4 different hands and 20 images from 6 different ex vivo pig organs in various configurations and orientations (intestine, colon, stomach, liver, esophagus and pancreas). None of the validation images were used for training. Images were acquired via 7-phase SFDI using 4 spatial frequencies (0.1 mm-1, 0.2 mm-1, 0.3 mm-1, 0.4 mm-1) to extract optical properties. These images were processed using four different processing modalities: SFDI, SSOP Filtering, SSOP Deep Learning (noted "SSOP DL conv") with convolutions of 7x1 and 7x7. For the quantitative performance assessment, 3D profile and optical properties maps obtained using SSOP methods at each spatial frequency were compared to that of 7-phase SFDI. These comparisons were undertaken by measuring the mean absolute percentage error given by the following Eq. (7.1) where N and M refer to image size (number of pixels in x and y directions):

() For the visual quality assessment, structural similarity index (SSIM) of 3D profile and optical properties maps obtained using SSOP methods at each spatial frequency were compared to that of 7-phase SFDI [12].

For the GPU implementation performance assessment, NVIDIA Nsight environment was used to trace the C CUDA code processing time. This time does not account for the transfer time between the CPU and the GPU memories (which is 338 µs for a 1024 x 1024 pixels size image) on our system. The evaluation of the quality of the images was performed visually to consider degradations, such as ripples and edge artifacts. Finally, a movie of a moving hand was acquired at fx = 0.2 mm -1 with rotation and translation to compare and contrast the three SSOP analytic method's real-time capabilities.

Note that the codes, for each demodulation method (SFDI 7 phases, SSOP Filtering, SSOP DL conv 7x1 and SSOP DL conv 7x7), and all data (calibration phantoms, hands, and ex vivo sample) used in this work are available at www.healthphotonics.org in the "resources" section.

Results

The 24 images of the test dataset were processed as described before with the four different modalities. To illustrate these results, we show in Fig. Profile-corrected reduced scattering Profile-corrected reduced scattering coefficient errors increase with spatial frequency. SSOP DL Conv 7x7 still performs better with errors as low as 7.5 ±2.7%, followed by SSOP DL Conv 7x1 with errors as low as 8.3 ±4.8% and SSOP Filtering with errors as low as 9.3 ±2.5%.

Overall, qualitative and quantitative analyses using both image-to-image similarity and pixelwise error illustrate the edge improvements provided by the twin U-NET deep learning workflow as well as an enhanced accuracy in optical properties extraction.

Processing time

The analysis of the GPU processing times for the three different SSOP implementations is summarized in On our workstation, the SSOP Filtering method was the fastest with 1.13 ms, followed by Deep

Learning with 7x1 convolution with 13.81 ms and finally Deep Learning with 7x7 convolution with 18.01 ms. As expected, the computation time to obtain diffuse reflectance and to extract optical properties is the same for all method (using the same GPU-optimized LUT). Thus, the SSOP deep learning methods are only limited by the network inference time. Overall, since a video rate of 25 images per seconds is limited at 40 ms, all three processing methods demonstrate capabilities to be used in real-time environments. Note that processing time for the 7 phase SFDI method is not reported since the method is not suited for real-time imaging.

Conclusion and Perspectives

Perspectives

The work we have achieved is only the basis of what remains to be done. Our goal for this project was to provide the surgeon with an approved finished product in the form of a cart system that can easily be installed in an operating room and that will provide in real-time at least 4 physiological parameters (oxygenation level, lipids level, melanin level, water level,) and 2 structural parameters (scattering power, scattering amplitude). As such, a first complementary work would be to increase the number of wavelengths to 6, which will involve: (1) the design of sources and a control system for temporal modulation, (2) improving the acquisition system by choosing a camera that can reach at least 500 Fps to ensure real-time acquisition, (3) increasing the GPU computing power and improving our implementations, (4) design and implement a comfortable viewing system that allows easy reading and interpretation. The second perspective is the implementation of a multi calibration technique to increase the accuracy of the SSOP results as well as the improvement of the 3D profile correction technique by taking into account the angle correction algorithm as well.

The third perspective will be to design the first cart system and to carry out the first in vivo experiments to assess the performances and limitations of the system. These experiments will also be an opportunity to obtain relevant feedback from surgeons in order to improve the ergonomics and visual interface of the system. spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right).

| Appendix xxxiv Table A1.8. Mean error percentage over 10 measured hands for reduced scattering at 0.2 mm -1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right).

| Appendix xxxv Table A1.9. Mean error percentage over 10 measured hands for reduced scattering at 0.3 mm -1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right).

| Appendix xxxvi Table A1.10. Mean error percentage over 10 measured hands for reduced scattering at 0.4 mm -1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right).

| Appendix xxxvii Table A1.11. Mean error percentage over 10 measured hands for the measurement of vein full width at half-maximum (FWHM) size. Measurement performed at 0.2 mm -1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right).

| Appendix xxxviii Table A1.12. Mean error percentage over 10 measured hands for the measurement of vein full width at half-maximum (FWHM) size. Measurement performed at 0.3 mm -1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right).

| Appendix xxxix

Table A.13. Mean error percentage over 10 measured hands for the measurement of vein full width at half-maximum (FWHM) size. Measurement performed at 0.4 mm -1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right).

| Appendix xl

Appendix 2: GPU architecture and programming interface

To be able to do low level programming on GPU, it is essential to know its structure but also the bases of the API that we use which is in our case the C CUDA.

A2.1. GPU architecture: NVIDIA QUADRO K2200

In

Floating-point Operation Per Second) in double precision [1]. In the classification of NVIDIA's computing capabilities, it is rated 5.0, making it possible to know all its specifications and characteristics [2]. The chip is accompanied by 4096Mb of GDDR5 memory, on a 128-bit bus and has a bandwidth of 80.2Gb/s. The interface bus of the card is PCIe 2.0 x16.

If we look at its internal structure, it is constituted of [3]:

• a 2Mb level 2 cache,

• a GPC (Graphics Processings Cluster) containing 5 SMM (Streaming Multiprocessor Maxwell) as shown in Fig. A2.2.,

• two 64-bit memory controllers,

• and the Gigathread Engine that distributes thread blocks to each SMM. The internal structure of one of the SMM can be depicted as in Fig. A2.3.. We can see:

• The Polymorph Engine which is the heart of tessallation (consists of breaking polygons into smaller pieces) [4] which is of more interest to graphic developers.

• The instruction cache which is linked to 4 instruction buffers. Finally, there is another way to use CUDA without having to worry about the hardware architecture of the GPU and the thread configuration. Indeed, NVIDIA and other institutions provide libraries in various areas that allow for the rapid development of complex applications. Among these libraries we can mention: cuFFT, cuDNN, cuBLAS, cuSPARSE, cuRAND, NPP, etc. In our development, we used cuFFT and cuDNN libraries. • Enagnon Aguénounon, Foudil Dadouche, Wilfried Uhring and Sylvain Gioux, "Real-time optical properties and oxygenation imaging using custom parallel processing in the spatial frequency domain," Biomed. Opt. Express 10(8), 3916-3928 (2019).

A2.4. Example of simple programming in C CUDA

• Enagnon Aguénounon, Foudil Dadouche, Wilfried Uhring, Nicolas Ducros and Sylvain Gioux, "Single snapshot imaging of optical properties using a single-pixel camera: a simulation study," J. Biomed. Opt., 24 (7), 071612 (2019).

• Enagnon Aguénounon, Foudil Dadouche, Wilfried Uhring and Sylvain Gioux, "Single snapshot of optical properties image quality improvement using anisotropic 2D windows filtering," J. Biomed. Opt., 24 (7), 071611 (2019).

• Manon Schmidt, Enagnon Aguénounon, Amir Nahas, Murielle Torregrossa, Bruce J.

Tromberg, Wilfried Uhring and Sylvain Gioux, "Real-time, wide-field, and quantitative oxygenation imaging using spatiotemporal modulation of light," J. Biomed. Opt., 24 (7), liii Proceedings:

• Enagnon Aguénounon, Foudil Dadouche, Wilfried Uhring and Sylvain Gioux, "Real-

International conferences SPIE Photonics West BIOS 2020:

• Real-time processing and visualization of functional parameters in living tissue with 3D profile correction.

Enagnon Aguénounon, Foudil Dadouche, Wilfried Uhring and Sylvain Gioux.

• Real-time, quantitative, and wide-field oxygenation imaging platform for surgery.

Enagnon Aguénounon, Silvère Ségaud, Henrique Waxin, Lucile Zorn, Julien Lamy, Murielle Torregrossa, Joseph Angelo, Sylvain Gioux.

European Conferences on Biomedical Optics -ECBO 2019 (23-27 June 2019):

• Real-time multispectral optical imaging using GPGPU processing.

Enagnon Aguénounon, Foudil Dadouche, Wilfried Uhring and Sylvain Gioux.

•

Abstract

Diffuse optical imaging has a significant potential for image-guided surgery. One of its methods called SFDI for Spatial Frequency Domain Imaging allows to extract the optical properties of biological tissue from structured light and using these measurements acquired at several wavelengths we can retrieve physiological and structural information interpretable by the surgeon. This thesis on the development of a real-time imaging instrument for surgery is based on SSOP "Single Snapshot imaging of Optical Properties" the rapid acquisition implementation of SFDI. To achieve our goals, we designed a system that uses spatiotemporal light modulation, where spatially we encode SSOP and temporally we encode wavelengths. Among our developments, we have proposed two methods, one in Fourier domain and one with deep learning for improving the visual quality and accuracy of SSOP results. A two-wavelength system (665nm and 860nm) was designed and we then validated it in vivo to measure the oxygen saturation rate and demonstrate that it can detect the ischemia of a tissue. Finally, to ensure real-time processing and visualization, all algorithms developed were implemented on GPU (Graphics Processing Unit). In conclusion, today we have the first imaging system capable of visualizing saturated oxygenation level in real-time and over a wide field of view. As perspectives, we plan to evolve the system to measure lipids, melanin and water levels and then validate it in pre-clinical studies.