mcxyz.c

Here is an update of the mcxyz.c program, as a zip file: mcxyzOct10_2014.zip, which contains
readme.txt
maketissue.m
	makeTissueList.m
	spectralLIB.mat
	makecmap.m
	reportHmci.m
mcxyz.c
lookmcxyz.m
	makec2f.m

This Oct. 10 update also adds a boundary flag to the input file (maketissue.m and mcxyz.c are updated as well):

  1. boundaryflag = 0 --> no boundaries. Photons propagate outside the cube using the optical properties correspond to the boundary voxels of the cube. For example, if photon is at x,y,z where x is outside the cube, then the optical properties of the tissue at the boundary value of x for the y and z positions is used for photon propagation. The weight is decremented appropriately, but no energy deposition into a bin occurs.
  2. boundaryflag = 1 --> photons escape at all boundaries.
  3. boundaryflag = 2 --> photons escape at surface only, and other boundaries are infinite (as in boundaryflag = 0).

NOTE: The Oct. 10 update also fixes the isotropic launch. (thanks to Baki (Reheman Baikejiang) for noticing the error.

Here is a maketissue_iso_pt.m file that demonstrates an isotropic point source.

The readme.txt file provides brief instructions on how to compile and run the program.
The program produces the following figures:

(Left) The tissue structure T(z,x) for skin with 200-µm-dia. blood vessel at 500µm depth. The rays of beam are shown as red lines.
(Middle) Fluence rate F(z,x)@y=0 [1/cm2], looking down the axis of vessel.
(Right) Fluence rate F(z,y)@x=0 [1/cm2], showing side view of vessel.
A 10-min computation is shown. Light delivered at tissue surface at z = 0.0100 cm, with water coupling above showing escaping light flux (reflectance).

Preliminary BETA version, March 4, 2013.
Will be updated with new version already being tested (Nov. 21, 2013). Appreciate any feedback: jacquess@ohsu.edu

All boundary conditions are matched. Does not yet handle mismatched boundaries.

The program mcxyz.c is a computer simulation of the distribution of light within a complex tissue that consists of many different types of tissues, each with its own optical properties. The program uses the Monte Carlo method of sampling probabilities for the stepsize of photon movement between scattering events and for the angles (θ,ψ) of photon scattering.

mcxyz.c is written in ANSNI standard C, and is easily compiled by any standard compiler. It reads input files (myname_T.bin, myname_H.mci) created by a MATLAB program maketissue.m, using a standard format. Therefore, other programs can also prepare the input files.

The 3D Monte Carlo generates an output file of relative fluence rate, F(y,x,z) [W/cm2 per W delivered] or [1/cm2].


Figure 1: An example showing the results of a simulation that includes air, epidermis, dermis, and a blood vessel. The tissue list also includes skull, gray matter and white matter (brain) but these tissue types are not used in this example. Green 532-nm light is delivered to the air/epidermis surface as a uniform-field 0.4-cm-diameter circular collimated beam. A 0.1-cm-diameter blood vessel is located at x,z = 0,0.4 cm, and runs along the y axis. In this simulation, there are 400 bins in x,y,z, i.e., T(400,400,400) is used to generate F(400,400,400). The program was run for 10 min, launching 310,092 photons.

A.
The tissue types T(z,x)@y are shown. The red lines show the photon paths if there were no scattering or absorption. The blood vessel extends along the y axis, and its cross-section is seen.

B.
The relative fluence rate φ(z,x)@y [1/cm2] is shown. The fluence rate of light escaping out of the skin at the epidermal surface into the air is also shown.
( log10(φ [cm-2]) ).

C.
The relative fluence rate φ(z,u)@x [1/cm2] is shown, where x equals the middle of the blood vessel, which appears as the horizontal bar of black at z = 0.4 cm due to strong absorption by blood.
( log10(φ [cm-2]) ).


In the following sections, you will learn

  1. how to use makeTissueList.m (MATLAB) to create a library of tissue optical properties
  2. how to use maketissue.m (MATLAB) to create the inputfiles (myname_T.bin, myname_H.mci) for use by mcxyz.c
  3. how to use mcxyz.c
  4. how to use lookmcxyz.m (MATLAB) to look at the results.

You can download a .zip file containing the 4 programs listed above, and supporting files, and the example simulation of Figure 1: mcxyzAug10_2014.zip.

Finally, some example simulations will be presented to illustrate how to adjust simulations.


How to use makeTissueList.m to create a library of tissue optical properties

To create a library of tissue types, with the associated tissue optical properties, the file spectralLIB.mat is used. This file contains the optical absorption properties, μa [cm-1], of

as functions of wavelength from 300 to 1000 nm. You can download spectralLIB.mat with the mcxyz.zip file. You can view a text version of the file here: spectralLIB.dat.

To create the scattering properties, a generic function is used:

μs' = μs.500nm' (fRayleigh(λ/500nm)-4 + fMie(λ/500nm)-bMie)

where

Each tissue type is specified as follows, using the 4th tissue type, dermis, as an example:

    j = 4;
    tissue(j).name = 'dermis';
    B	= 0.002;	% blood volume fraction
    S	= 0.67;		% oxygen saturation of hemoglobin
    W 	= 0.65;		% water volume fraction
    F	= 0.02; 	% fat volume fraction
    M 	= 0;		% volume fraction of melanosomes
    musp500 = 42.4;  	% reduced scattering coeff. at 500 nm
    fray    = 0.62;  	% fraction of Rayleigh scattering at 500 nm
    fmie 	= 1-fray;	% fraction of Mie scattering at 500 nm
    bmie    = 1.0;  	% scatter power for mie scattering
    X       = [B*S B*(1-S) W F M]';	% [oxyBlood deoxyBlood Water Fat Melansomes]'
    musp    = musp500*(fray*(nm/500).^-4 + fmie*(nm/500).^-bmie);
    gg      = 0.90;
    tissue(j).mua = MU*X;		% absorption coefficient
    tissue(j).mus = musp/(1-gg);	% scattering coefficient
    tissue(j).g   = gg;		% anisotropy of scattering
The above calculates the absorption (tissue(j).mua), the scattering (tissue(j).mus) and the anisotropy of scattering (tissue(j).g)) at the wavelength of interest. By adjusting B,S,W,F and M, the absorption properties are specified. By adjusting musp500, fray, bmie, the scattering properties are specified.

The function call tissue = makeTissueList(nm), where nm specifies the wavelength of interest in [nm], will yield a data structure called tissue with the structure:

where j selects the jth tissue type in the library. The parameter tissue(j).nm is the same single wavelength for all tissues, but is included for convenient recall.

In the example makeTissueList(nm) presented here, 9 tissue types are specified:

  1. air
  2. water
  3. blood
  4. dermis
  5. epidermis
  6. fat or lipid
  7. skull
  8. gray matter (brain)
  9. white matter (brain)

For both air and water, a moderate scattering coefficient is specified, but the value of g is set to 1.0. Also a very very low absorption coefficient is specified (eg., μa = 0.0001 cm-1). Hence, the photon will step through air or water, and deposit a very small amount of photon weight in the voxels, thereby specifying its pathlength spent in the air or water. But the energy lost in the air or water is negligible, so it does not significantly influence the distribution of light. At the end of the program (discussed in How to use mcxy.c, below), the fluence rate φ is calculated by dividing the depositive energy by the very small μa which recovers the φ in the air and water.

View listing of makeTissueList.m: makeTissueList.m


How to use maketissue.m to create a tissue

The program maketissue.m uses the tissue library tissue created by makeTissueList.m. The program creates a complex tissue, according to the user's adjustments of the following parameters in the program:

    %%% USER CHOICES %%%%%%%% <-------- You must set these parameters ------
    SAVEON      = 1;        % 1 = save myname_T.bin, myname_H.mci 
                            % 0 = don't save. Just check the program.

    myname      = 'example';% name for files: myname_T.bin, myname_H.mci  
    nm          = 532;   	% desired wavelength of simulation
    time_min    = 10;      	% time duration [min] of the simulation
    Nbins       = 400;    	% # of bins in each dimension of cube 
    binsize     = 0.0020; 	% size of each bin, eg. [cm] or [mm]

    % Set Monte Carlo launch flags
    mcflag      = 0;     	% launch: 0 = uniform beam, 1 = Gaussian, 2 = isotropic pt. 
    launchflag  = 0;        % 0 = let mcxyz.c calculate launch trajectory
                            % 1 = manually set launch vector.

    % Sets position of photon launch
    xs          = 0;      	% x of launch
    ys          = 0;        % y of launch
    zs          = 0.10;  	% z of launch

    % Set position of focus, so mcxyz can calculate launch trajectory
    xfocus      = 0;        % set x,position of focus
    yfocus      = 0;        % set y,position of focus
    zfocus      = inf;    	% set z,position of focus (=inf for collimated beam)

    % only used if mcflag == 0 or 1 (uniform or Gaussian beam, not isotropic pt.)
    radius      = 0.2;      % 1/e radius of beam at launch position
    waist       = 0.010;  	% 1/e radius of beam at focus

    % only used if launchflag == 1 (manually set launch trajectory):
    ux0         = 0.7;      % trajectory projected onto x axis
    uy0         = 0.4;      % trajectory projected onto y axis
    uz0         = sqrt(1 - ux0^2 - uy0^2); % such that ux^2 + uy^2 + uz^2 = 1
    %%%%%%%%%%%%%%%%%%%%%%%%%
The program then enters the section
    %%%%%%%%%%
    % Prepare Monte Carlo
    %%%
which sets up the parameters for the Monte Carlo simulation. (View the original code.)

The user then creates the desired complex tissue in the next section:

    %%%%%%
    % CREATE TISSUE STRUCTURE T(y,x,z)
    %   Create T(y,x,z) by specifying a tissue type (an integer)
    %   for each voxel in T.
    %
    %   Note: one need not use every tissue type in the tissue list.
    %   The tissue list is a library of possible tissue types.

    T = double(zeros(Ny,Nx,Nz)); 

    T = T + 4;      % fill background with skin (dermis)

    zsurf = 0.100;  % position of air/skin surface

    for iz=1:Nz % for every depth z(iz)

        % air
        if iz<=round(zsurf/dz)
            T(:,:,iz) = 1; 
        end

        % epidermis (60 um thick)
        if iz>round(zsurf/dz) & iz<=round((zsurf+0.0060)/dz)
            T(:,:,iz) = 5; 
        end

        % blood vessel @ xc, zc, radius, oriented along y axis
        xc      = 0;            % [cm], center of blood vessel
        zc      = Nz/2*dz;     	% [cm], center of blood vessel
        vesselradius  = 0.0500;      	% blood vessel radius [cm]
        for ix=1:Nx
                xd = x(ix) - xc;	% vessel, x distance from vessel center
                zd = z(iz) - zc;   	% vessel, z distance from vessel center                
                r  = sqrt(xd^2 + zd^2);	% r from vessel center
                if (r<=vesselradius)     	% if r is within vessel
                    T(:,ix,iz) = 3; % blood
                end

        end %ix

    end % iz
In the above, the user fills the backround with dermis, specifies the position of the air/skin surface (zsurf), then steps through each z(iz) depth position and determines if in air or epidermis, assigning the appropriate integer pointer value to all x,y elements at the z(iz) depth in the tissue array T(iy,ix,iz). Also, the presence of the blood vessel at depth z(iz) for each of the x,y positions is tested, based on the vessel radius being a set value (vesselradius).

When ready to save the results to a binary file for mcxyz.c to use, the array T(iy,ix,iz) is reshaped into a linear array of integer values, v(i), by the command

    % convert T to linear array of integer values, v(i)i = 0;
    v = uint8(reshape(T,Ny*Nx*Nz,1));
and this array is saved as a binary file, using the user's choice myname:
    %% write myname_T.bin file
    filename = sprintf('%s_T.bin',myname);
    disp(['create ' filename])
    fid = fopen(filename,'wb');
    fwrite(fid,v,'uint8');
    fclose(fid);
The program also saves myname_H.mci which lists the simulation parameters and tissue list of optical properties for use by mcxyz.c. The program also generates a figure similar to Figure 1A that shows the structure.

View listing of maketissue.m: makeTissueList.m


How to use mcxyz.c

In Unix, to compile mcxyz.c is simply

	cc -o mcxyz mcxyz.c
which saves the executable version as mcxyz.

From a UNIX command line, the program is executed by the command

	mcxyz myname
where myname is the name chosen by the user when creating the tissue using maketissue.m. For example, Figure 1 used myname = "example1", and mcxyz used the files "example1_H.mci" and "example1_T.bin" to generate "example1_F.bin".

As the program runs, it will look like this:

	name = example1
	time_min = 10.00 min
	Nx = 400, dx = 0.0020 [cm]
	Ny = 400, dy = 0.0020 [cm]
	Nz = 400, dz = 0.0020 [cm]
	xs = 0.0000 [cm]
	ys = 0.0000 [cm]
	zs = 0.1000 [cm]
	mcflag = 0 [cm]
	xfocus = 0.0000 [cm]
	yfocus = 0.0000 [cm]
	zfocus = 1.00e+12 [cm]
	ux0 = 0.7000 [cm]
	uy0 = 0.4000 [cm]
	uz0 = 0.5916 [cm]
	radius = 0.2000 [cm]
	waist  = 0.0100 [cm]
	Nt = 8
	muav[1] = 0.0001 [cm^-1]
	musv[1] = 1.0000 [cm^-1]
	  gv[1] = 1.0000 [--]
	
	muav[2] = 0.1000 [cm^-1]
	musv[2] = 10.0000 [cm^-1]
	  gv[2] = 1.0000 [--]
	
	muav[3] = 230.5427 [cm^-1]
	musv[3] = 93.9850 [cm^-1]
	  gv[3] = 0.9000 [--]
	
	muav[4] = 0.4585 [cm^-1]
	musv[4] = 356.5406 [cm^-1]
	  gv[4] = 0.9000 [--]
	
	muav[5] = 16.5724 [cm^-1]
	musv[5] = 375.9398 [cm^-1]
	  gv[5] = 0.9000 [--]
	
	muav[6] = 0.1154 [cm^-1]
	musv[6] = 281.9549 [cm^-1]
	  gv[6] = 0.9000 [--]
	
	muav[7] = 2.3057 [cm^-1]
	musv[7] = 181.5859 [cm^-1]
	  gv[7] = 0.9000 [--]
	
	muav[8] = 2.3057 [cm^-1]
	musv[8] = 181.5859 [cm^-1]
	  gv[8] = 0.9000 [--]
	
	central axial profile of tissue types:
	22222222222222222222222222222222222222222222222222555444444444444444444444444444
	44444444444444444444444444444444444444444444444444444444444444444444444444444444
	44444444444444433333333333333333333333333333333333333333333333334444444444444444
	44444444444444444444444444444444444444444444444444444444444444444444444444444444
	44444444444444444444444444444444444444444444444444444444444444444444444444444444
	
	Sun Mar  3 10:40:27 2013
	
	------------- Begin Monte Carlo -------------
	example
	requesting 10.0 min
	Nphotons = 310092 for simulation time = 10.00 min
	1% done
	2% done
	3% done
	4% done
	5% done
	6% done
	7% done
	8% done
	9% done
	10% done
	20% done
	30% done
	40% done
	50% done
	60% done
	70% done
	80% done
	90% done
	91% done
	92% done
	93% done
	94% done
	95% done
	96% done
	97% done
	98% done
	99% done
	100% done
	------------------------------------------------------
	Elapsed Time for 3.101e+05 photons = 9.703 min
	3.20e+04 photons per minute
	saving example_F.bin
	example is done.
	------------------------------------------------------
	Sun Mar  3 10:50:14 2013
First, the simulation parameters and tissue optical properties are listed. Then the tissue types along th central axis versus depth are listed as numerals, which provides a quick visual check on the tissue structure. Then the Monte Carlo simulation begins, giving the % of task completed throughout the run. When finished, the file myname_F.bin is saved.


How to use lookmcxyz.m to look at results

Now that the simulation is completed, the function lookmcxyz.m will view the results. In MATLAB, assign the prefix name of our files, eg., "example1", to the variable "myname". Then simply type the command lookmcxyz.m:

	myname = "example";
	lookmcxyz
and in this case the program will load the files example1_T.bin and example1_H.mci, and prepare the 3 figures in Figures 1A,B,C.

The program master_mcxyz.m is a convenient program that guides you in calling lookmcxyz.m .

View listing of master_mcxyz.m: master_mcxyz.m and lookmcxyz.m: lookmcxyz.m


EXAMPLES


Figure 2: An example like Figure 1, but with the epidermis and dermis replaced with a background tissue with μa = 0.1 cm-1, μs = 1 cm-1, and g = 0.90, and the light focused at a 0.2-cm depth. With the scattering reduced, the focus is clearly seen, and the blood vessel casts a sharp shadow. The program was run for only 1 min, launching 225,176 photons.

A.
The tissue types T(z,x)@y are shown. The red lines show the photon paths if there were no scattering or absorption. The blood vessel extends along the y axis, and its cross-section is seen.

B.
The relative fluence rate φ(z,x)@y [1/cm2] is shown. The fluence rate of light escaping out of the skin at the epidermal surface into the air is also shown.
( log10(φ [cm-2]) ).

C.
The relative fluence rate φ(z,u)@x [1/cm2] is shown, where x equals the middle of the blood vessel, which appears as the horizontal bar of black at z = 0.4 cm due to strong absorption by blood.
( log10(φ [cm-2]) ).



Figure 3: An example that mimics a small 3-cm-dia. head, with epidermis, dermis, skull, gray matter, white matter, and a central ventricle (cerebral spinal fluid ≈ water). The green light (532 nm) is delivered as a 4-mm-dia. circular beam to the top surface of the head. The backscatter of light into the air is shown. (A 10-min simulation.)

A.
The tissue types T(z,x)@y are shown. The red lines show the photon paths if there were no scattering or absorption.

B.
The relative fluence rate φ(z,x)@y [1/cm2] is shown. The fluence rate of light escaping out of the skin at the epidermal surface into the air is also shown.
( log10(φ [cm-2]) ).

C.
The relative fluence rate φ(z,u)@x [1/cm2] is shown.
( log10(φ [cm-2]) ).


The MATLAB code in maketissue.m that created this model head in Figure 3 is listed:

%%%%%%
% CREATE TISSUE STRUCTURE T(y,x,z)
%   Create T(y,x,z) by specifying a tissue type (an integer)
%   for each voxel in T.

T = double(zeros(Ny,Nx,Nz)) + 1; % background is air.

% location of the center of the head [cm]
xc = 0;
yc = 0;
zc = 2.5;

% x,y,z already specified.

[XX YY ZZ] = meshgrid(x,y,z);
RR = sqrt((XX-xc).^2 + (YY-yc).^2 + (ZZ-zc).^2);
clear XX YY ZZ

disp('creating epidermis')
mask = RR<1.5;
T = T.*(~mask) + mask*5; % epidermis

disp('creating dermis')
mask = RR<1.48;
T = T.*(~mask) + mask*4; % dermis

disp('creating skull')
mask = RR<1.40;
T = T.*(~mask) + mask*6; % skull

disp('creating gray matter')
mask = RR<1.20;
T = T.*(~mask) + mask*7; % gray matter

disp('creating white matter')
mask = RR<1.00;
T = T.*(~mask) + mask*8; % white matter

disp('creating ventricles (cerebral spinal fluid)')
mask = RR<0.25;
T = T.*(~mask) + mask*2; % water

clear RR


Figure 4: Same as Figure 3 but the light is delivered as a 200-μm-dia fiber inserted wihtin the white matter. (A 10-min simulation.)

A.
The tissue types T(z,x)@y are shown. The red lines show the photon paths if there were no scattering or absorption.

B.
The relative fluence rate φ(z,x)@y [1/cm2] is shown. The fluence rate of light escaping out of the skin at the epidermal surface into the air is also shown.
( log10(φ [cm-2]) ).

C.
The relative fluence rate φ(z,u)@x [1/cm2] is shown.
( log10(φ [cm-2]) ).