2.1.

Diffusion Equation and the MPDP

Perelman et al. treated the propagation of light in a turbid medium as a free random walk.^{9 10 11} They used the path integral with a quasiparticle Lagrangian to find the most probable trajectory, or classical path, over which photons can be found. Their theory can also be expressed as follows: for the given initial and final states and transport time, the probability density of finding a photon at the most probable path is maximized. As we see here, the most probable trajectory was apparently described in the time domain.

In reality, compared with the time-domain (TD) measurements, the frequency-domain (FD) method is simpler and more reliable in terms of data interpretation and immunity from noise. Medical FD equipment is more economical and portable.¹² More attention is currently being paid to the FD method. Therefore, in our work, we defined the MPDP in the FD.

In the framework of wave optics, energy transfer is along the gradient direction of the density of electromagnetic energy. However, energy transfer can also be described by the classical movement of particles in the framework of geometrical optics. Therefore, we use the same two frameworks to describe our defined quasiparticle, i.e., the particle, which represents the average effect of wave progress at the macroscale.

The diffusion equation is described in detail later, as the basis of our theory.

Light propagation is usually described by diffusion approximation in a highly scattering medium, when the source frequency is¹³ less than 1 GHz. In the FD, according to the diffusion approximation, the flux vector j(r,ω) can be expressed¹⁴ as the gradient of fluence rate Φ(r,ω):

Equation (1) represents a proportional relationship between the macroscopic flow of the scattering particles and the spatial distribution of the particle density, which is analogous to the flow conservation law of fluid dynamics. Substituting Eq. (1) into the radiative transfer equation⁷ (RTE), the RTE can be simplified as the diffusion equation, written as

The DPDW described by Eq. (4) is a new type of wave used to describe the random migration of diffuse photons in turbid media. The properties of the DPDW are different from those of general waves, such as sound and light waves. For general waves, the square of their amplitude is of the order of the magnitude of energy density, while the amplitude of DPDW itself is of the order of the magnitude of energy density. The field variables of the general wave are represented by complex numbers in the FD, but the energy density and flux are real numbers, while the energy density and flux of DPDW are complex numbers. What we are interested in is the propagation of information carried by the alternating component of energy current.

For convenience of calculation, we use j _info (r,ω) to take the place of j(r,ω) in Eq. (1). We consider the amplitude of Φ(r,ω), thus j _info (r,ω) becomes

Usually in FD method, light sources are modulated with sine waves of a single frequency. In this paper, we define the ac components of the fluence rate and flux with respect to a fixed angular frequency. Note that a dc component always exists accompanying the ac component (for the light source intensity, fluence rate, detected flux, and so on). However, we are only interested in the ac component (more specifically, the amplitude of the ac component) in this paper.

As we mentioned earlier, Perelman et al. first introduced the concept of the most probable trajectory (or classical path) in their work. Here in our work, we propose the most probable path in the steady state. In the steady state, for the given distribution of the fluence rate, and the injection point and exit point (or start and end point) the MPDP is defined as the path on which a scattered photon is most likely to occur. For each point on the MPDP, the direction of MPDP is the direction of the flux vector of this point. The MPDP can be considered as a random walk process from the injection point to the exit point. The step length is the transport mean free path, determined by l _tr =3D. The moving direction at point r is determined by j _info (r,ω) in Eq. (5).

Our definition of the most probable path in the FD does not come from photon trajectories, and is different from the TD approach. Instead, we start directly from the diffusion equation and try to localize the most important energy (information) flow carried by the diffusive photon density waves across the turbid media. Such localization is represented by the trajectory of a conceptual quasiparticle, which is not a real particle at all.

Now we must calculate the MPDP in a transmission mode for a medium with finite thickness, since the transmission mode is frequently adopted for breast cancer diagnosis.

Before giving our model, we recall the well-accepted semi-infinite model first. For a semi-infinite turbid medium, the perpendicularly incident light can be replaced by an effective source and an image source.^{16 17} The effective source is an isotropic point source located an APD (average penetration depth) into the medium. The image source is the image of the effective source by the extrapolated boundary. If we let the residual strength (RS) of the effective source be 1, then the RS of the image source will be −1.

For the transmission mode, the thickness of tissue is an important factor, i.e., the medium we studied must have at least one finite dimension, e.g., the Z direction. For simplicity, in the MPDP calculation, we assume that the quasiparticle travels in the X-Z plane. This is the MPDP in two dimensions. From Fig. 1, we can see that the 2-D model is infinity in the X direction, while the thickness in the Z direction is defined as H. As shown in Fig. 1, the light enters into the medium from the original point along the Z direction. There are two extrapolated boundaries in our model, which is different with the half-infinite model.

Figure 1

The 2-D tissue model and the boundary conditions.

In our model, image method¹⁶ is used to construct a fluence rate solution that satisfies Φ(x,z=−Δ)=0 and Φ(x,z=H+Δ)=0. The two extrapolated boundaries work like two mirrors. For the effective source with the unit strength at (0,APD), there are infinitely many image sources created by the two extrapolated boundaries at [x=0, z=2N(H+2Δ)±(Δ+APD)−Δ], where N is an integer. Figure 1 shows the 2-D tissues model and the boundary conditions with N=0. Different Ns give different image sources. In Fig. 1, Δ is the extrapolation distance. We chose Δ=APD=l _tr during calculation.¹⁷ Although the image sources are at infinity, only the first 10 to 20 terms are mainly used. The contribution of the other image sources is greatly diminished by the far distance.

According to Eq. (4), we can calculate the 2-D distribution of fluence rate in a 3-D medium by superposing the fluence rate of DPDWs generated by the effective source and image sources. Then the calculated fluence rate is substituted into Eq. (5) and the flux vector in the 2-D area can be obtained. First, a vertically injecting quasiparticle comes to (0,l _tr ). Second, the vertically injecting quasiparticle will randomly take an initial angle and move to the next step. The step length is always l _tr . The following points in the medium follow the direction of j _info (r,ω) given by Eq. (5). When the quasiparticle is close enough to the surface, i.e., less than l _tr to the surface, the path will end at the surface of the medium, either on top or at the bottom. The path ending at the detection point is the one we require. The MPDP is affected by several factors: μ_a, μ_s ^′, the refractive index n, H, the modulation frequency ω of the source, and the relative position between the source and the detector, defined as the lateral distance dx. From our model, the MPDP we calculate is not the analytic function of the trajectory, but a series of points on it. Therefore, for the convenience of calculation, a polynomial is created from these points to describe the MPDP. The function z=f(x,μ_a,μ_s ^′,n,H,dx,ω) is used to describe the polynomial, abbreviated as z=f(x).

Figure 2 shows one example of the transmitting MPDP. The source is at (X,Z)=(0,0), while the detector is at (5,50) (the units are millimeters). The parameters are μ_a=0.005 mm ⁻¹, μ_s ^′=1 mm ⁻¹, n=1.33, H=50 mm, and dx=5 mm. When f=49.86 MHz, the MPDP can be described by a polynomial as

Figure 2

MPDP at different modulation frequencies. Solid line: dc; +: 50 MHz; ○: 200 MHz. A dashed straight line is also plotted for comparison.

The calculated MPDP takes on an S shape. The symbols “○” and “+,” and the solid lines show the path at different source frequencies: 200 MHz, 50 MHz, and 0 Hz (the direct current condition), respectively. From Fig. 2 we can see that the higher the source frequency, the closer is the path to a straight line. The straight dotted line represents the direct line between the source and the detector.

2.2.

Localization of Inhomogeneities Using MPDP

Generally a projection method is useful to detect the inhomogeneities in tissue, e.g., a mammogram is a projection of the breast in film. The NIR light source can also be used to make a projection of breast tissue,¹⁸ referred to as optical mammography. However, optical mammography can acquire depth information concerning breast tissue, while conventional mammogram, x-ray mammography, cannot. This is because that the projection data of the NIR light are obtained by the point-to-point scanning of a source-detector pair. The projection data thus represent the information in a matrix of the sampling points. By using some reconstruction methods, the depth information concerning the breast tissue can be obtained.

In our work, we propose a method with the MPDP in the steady state to locate an absorber in three dimensions, which includes determining the depth information of the absorber. Figures 3(a) and 3(b) are used to describe our method in three and two dimensions, respectively. We use L, W, and H to represent the length, width, and height of the medium, respectively; and d is the distance between the two detectors. For simplicity and to describe our method, we assume in Fig. 3 that the MPDP is a straight line between the source and the detector. In Fig. 3(b) there are one source and two detectors in the plane parallel to the X-Z plane, which means the source-detector triangle is always in the plane parallel to the X-Z plane. In this example, the source is at (x_s,y_s,0), the left detector is at (x_s−d/2, y_s,H), and the right detector is at (x_s+d/2, y_s,H). At y=0, the source-detector triangle scans the first row along the X direction. Then y is increased by an incremental, the source-detector triangle scans the second row along the X direction. This processes goes on until the last y is reached, and the whole scanning area has been scanned. When the source and detectors scan the sample, the absorber will shut out the light twice at position 1 (dashed lines) and position 2 (solid lines). At position 1 (X _PR ,Y _PR ,0), the source, the absorber, and the right detector are in one straight line, and the right detector shows an absorption peak. At position 2 (X _PL ,Y _PL ,0), the source, absorber, and the left detector are in one straight line, and the left detector shows an absorption peak.

Figure 3

Locating the absorber in three dimensions with the MPDP: (a) localization system in transmission mode and (b) MPDP shown by a straight line.

After scanning in the X-Y plane, two optical projections are obtained from the two detectors respectively; each has a shadow of the absorber. Ideally, the absorber center of the left projection is at (X _PL ,Y _PL ), and the right projection is at (X _PR ,Y _PR ). We assume that the projection coordinate system is the same as the source coordinate system. When the source is at (X _PL ,Y _PL ,0) and (X _PR ,Y _PR ,0), respectively, there are two MPDPs that pass through the absorber and end at the left detector (X _PL −d/2,Y _PL ,H) and the right detector (X _PR +d/2,Y _PR ,H), respectively. The question is how to reconstruct the real coordinates [defined as (X_r,Y_r,Z_r)] of the absorber from the absorber centers in the projection image, i.e., (X _PL ,Y _PL ) and (X _PR ,Y _PR ). According to the geometry of the source-detector triangle Y _PL should be equal to Y _PR in theory, and X_r should be in the middle of X _PL and X _PR . Considering the system error, the average of Y _PL and Y _PR will be used as Y_r:

For the S-shaped MPDP, we plug ΔX/2 into the polynomial function mentioned earlier to obtain the depth of absorber:

2.3.

Image Processing

After scanning we can acquire two data matrices, referred to as M _L and M _R. Each of the matrices has N (in the X direction, N=9)×M (in the Y direction, variable) points. The subscript L represents the detector at the “left” side of source and R represents the “right.” The notation M_L(i,j) represents the element in the line i and column j of the matrix of the left projection, and i represents the ordinal of the source along the X direction and j represents the step in the Y direction.

To improve the resolution in the X direction, we combined M _L and M _R to make one matrix M. The following equations describes how to combine the two matrices:

The image segmentation method is also used to improve the resolution of the image. Set a threshold and reset the pixels with a value that is less than the threshold. The nonzero part will be isolated by the reset pixels, so different parts of the image will be partitioned, i.e., segmented. It is important, but sometimes hard to choose the threshold.

It has been found that, in the image of diffuse light amplitude, the blurring level of the absorbers’ edge is associated with the depth of absorbers. For absorbers in the middle level of tissue, a threshold close to the absorption peak should be chosen.

The combination method and the segmentation method make the shadow of the absorbers clearer, which is helpful in finding the size and shape of the absorbers with more accuracy. The results of using these two methods were shown by the combined image and the deblurred image, respectively.

All the imaging and reconstruction functions used in this paper were developed in the Matlab 5.3 environment.

2.4.

Experiment

Based on our theory, we developed a prototype in our lab. The prototype was designed to detect breast tumor. The 780-nm NIR light generated by laser diode (Sharp Co. LT027) was modulated at 49.86 MHz and carried by an optical fiber bundle to the surface of the sample. The optical fiber bundles also conduct the signal to photomultiplier tubes (Hamamatsu R928). To implement the detection method described in Sec. 2.2, the prototype has 19 optical fiber bundles in total. As shown in Fig. 4, 9 of them in one row are the light sources, while the other 10 in another row are detectors. They work as nine groups of source-detector triangles, as shown in Fig. 3. The nine groups are arranged in one row along the X direction in the plane vertical to the Y axis. The distance of the adjacent groups is 6 mm and the lateral distance between source and detector (dx) in each group is 3 mm. An electrical switch took the place of a mechanical scan, which can save time. Nine sources lit alternately just like the source moved along the X axis in 6-mm steps. At one time, only one source was lit and two detectors (nearest to the lit source) received the signal. The inner diameter of each fiber bundle is 2 mm. Using the prototype, the amplitude information of the penetrating light was detected and used to image and reconstruct.

Figure 4

Sources and detectors.

In our experiments, Intralipid and pork were used to simulate breast tissues. Tumors often absorb more lights than normal tissues, so a black plastic ball or some fictile objects of different shapes and sizes or muscle (from pork) were used to simulate them. The imaging area was 45 mm wide (in the X direction) and 57 mm long (in the Y direction). The coordinate system used here was the same as the coordinate system described in Figs. 3 and 4. The z=0 point was on the top surface of the phantom. The initial position of the first light source (the first source from the left in Fig. 4) was used as the original point in the X-Y plane in each experiment.

During the experiments, the scan area was kept in the center of the phantom, about 10 mm away from the nearest edge. The absorbers were also kept at least 10 mm away from the edges of the phantom.

3.1.

Experiments with Intralipid

The reconstruction coordinates of the single absorber experiments are shown in Table 1. Experiment 2 has the largest error in reconstruction coordinate z since the absorber is the farthest from the surface of the solution among the three experiments. Figures 5 and 6 show the combined image and the deblurred image of experiments 2 and 3. The circle in Fig. 5 and the white rectangle in Fig. 6 show the real position and outline of the absorber. The prototype cannot only locate the absorber in Intralipid clearly, but can also estimate the size of the absorber. An absorber as small as a 5-mm-long cylinder can be detected, although the shadow does not look like the absorber.

Figure 5

Images of experiment 2: (a) combined image and (b) deblurred image.

Figure 6

Images of experiment 3: (a) combined image and (b) deblurred image.

Table 2 shows the experimental results from an absorber pair. It is difficult to distinguish the two absorbers. Our algorithm gives the coordinates of one absorption peak, which is the middle point of the absorber pair. Figures 7 and 8 show the images of experiments 4 and 5. After the images were deblurred, in Fig. 7 the two absorbers could be distinguished and their size could be estimated. The gray-scale values of the two absorbers are not the same in Fig. 7; this may be due to the mismatch of sources and detectors. In Fig. 8, the two absorbers have one merged shadow. Comparing the images of experiment 5 with experiment 4, we found that the shadows of the absorbers near to the surface are clearer than those farther from the surface.

Figure 7

Images of experiment 4: (a) combined image and (b) deblurred image.

Figure 8

Images of experiment 5: (a) combined image and (b) deblurred image.

Figure 9 shows the comparison of the two methods (the straight line versus the S-shaped MPDP). The horizontal axis represents the real depth (Z _org ) of the absorber. The vertical axis represents the reconstruction depth Z_r. The depth calculated from Eq. (9) is represented by “+” and the depth calculated from Eq. (10) is represented by “○.” The straight line is the ideal reconstruction. Marker “○” is closer to the straight line than marker “+.” This means that the MPDP in the S shape is a better description of the DPDW than the straight line.

Figure 9

Comparison of the two different methods.

3.2.

Experiment with Pork

Table 3 shows the experiments of absorbers in fat. As shown in Fig. 10, the 8-mm black ball in fat is very clear. Figure 11 shows the two pieces of muscle in fat. The reconstruction result gives one absorption peak, which points to the larger piece of muscle. After the images are deblurred, the shadows of the two pieces of muscle are distinguishable. This result is encouraging, for the absorption coefficients of muscle and fat are similar, while the distance between the edges of the two pieces of muscle is only 2 to 3 mm.

Figure 10

Images of experiment 6: (a) combined image and (b) deblurred image.

Figure 11

Images of experiment 7: (a) combined image and (b) deblurred image.