tmcimg.xhtml

蒙特卡罗模拟光子成像C语言版,代码简洁专业

    <tr>
	<td>nxvox</td>
	<td>required</td>
	<td>Total number of voxels in X</td>
    </tr>
    <tr>
	<td>nyvox</td>
	<td>required</td>
	<td>Total number of voxels in Y</td>
    </tr>
    <tr>
	<td>nzvox</td>
	<td>required</td>
	<td>Total number of voxels in Z</td>
    </tr>
    <tr>
	<td>image_x</td>
	<td>0..nxvox-1</td>
	<td rowspan="3">Define a region of interest</td>
    </tr>
    <tr>
	<td>image_y</td>
	<td>0..nxvox-1</td>
    </tr>
    <tr>
	<td>image_z</td>
	<td>0..nxvox-1</td>
    </tr>
    <caption>Table 1. Recognized keywords.</caption>
</table>

<h3 id="sub-segmentation_file">1.2 Segmentation File</h3>

<p class="noindent">
The segmentation file is a contiguous block of data, one unsigned
8-bit number for every voxel being modeled.  A value of '0' is used
for "air"; any photon making its way into the air is assumed to have
left the system.  A photon is recorded by a given detector if it exits
into "air" and its final position is within the specified radius of
the detector location.  Values greater than one correspond to tissues
whose optical properties were specified in the configuration file; the
first tissue defined is '1', the second tissue defined is '2', and so
on.  An error occurs if the value in any voxel of the segmentation
file exceedes the total number of defined tissue types.
</p>

<p>
<strong>Important note:</strong> by default, <code>tMCimg</code>
assumes that the segmentation data was generated by Matlab, which uses
column-major storage instead the usual C/C++ row-major storage.  If
you are generating the segmentation map from a C program, you must
either transpose the X and Y coordinates when you write out the data
or use the "-R" flag when running <code>tMCimg</code>.  See
Figures&nbsp;<a href="tMCimg.xhtml#fig-matlab_ex">2</a> and 
<a href="tMCimg.xhtml#fig-C_ex">3</a> for sample Matlab and
C/C++ code respectively.
</p>

<div class="padbox">
    <pre id="fig-matlab_ex">
>> boundData = uint8(zeros(40,50,10));
>> boundData(:,:,2:9) = 1;
>> % Any other initializations would go here
>> fid = fopen('sample.bin', 'wb');
>> fwrite(fid, boundData, 'uint8);
>> fclose(fid);
    </pre>
    <div class="divcaption">
    Figure 2. Sample Matlab code for writing a simple "slab"
    boundary geometry.  Normally, source would be located in the 
    <math xmlns="&mathml;"><mi>z</mi><mo>=</mo><mn>1</mn></math>
    plane and detectors in either the 
    <math xmlns="&mathml;"><mi>z</mi><mo>=</mo><mn>1</mn></math>
    or <math xmlns="&mathml;"><mi>z</mi><mo>=</mo><mn>10</mn></math>
    plane for this type of boundary arrangement.  Because there is no
    layer of 0's around the edges in
    <math xmlns="&mathml;"><mi>X</mi></math> or 
    <math xmlns="&mathml;"><mi>Y</mi></math> in this example, photons 
    leaving through the sides of the box will be lost.
    </div>
    <hr />
    <pre id="fig-C_ex">
unsigned char boundData[40][50][10];
int i, j, k;
FILE *fp;

/* Initialize boundData.  The for-loops could probably
 *  be replaced with a memset(), but doing it this way
 *  is easier to follow */

memset(boundData, 0, sizeof(boundData));

for (k = 1; k &lt; 9; k++)
  for (j = 0; j &lt; 50; j++)
    for (i = 0; i &lt; 40; i++)
      boundData[i][j][k] = 1;

/* Any other initializations would go here */

fp = fopen("sample.bin", "wb");

for (k = 0; k &lt; 10; k++)
  for (i = 0; i &lt; 40; i++)
    for (j = 0; j &lt; 50; j++)
      fputc(boundData[i][j][k], fp);

fclose(fd);
    </pre>
    <div class="divcaption">Figure 3. Sample C/C++ code for writing
    the same segmentation file as in
    <a href="tMCimg.xhtml#fig-matlab_ex">Figure&nbsp;2</a>.
    </div>
</div>

<h2 id="sec-output_files">Output File Format</h2>

<p>
The output of the program consists of two datafiles: the history
file and the two-point file.
</p>

<h3 id="sub-history_file">History Files</h3>

<p>
The history file contains a record of every photon that successfully
traversed the volume.  There are two different file formats for the
history file, based on the settings in the <code>*.cfg</code> file.
If no detectors were specified, then the (x,y,z) location of every
photon exiting the medium is recorded together with the time and the
distance the photon traveled inside each tissue type.  If at least one
detector was specified, the number of the detector that saw the photon
is recorded, along with the time and lengths as before.
</p>
<div class="padbox">
    <math xmlns="&mathml;">
	<mtable equalcolumns="false" align="center"
		rowalign="center" columnalign="center">
	    <mtr>
		<mtd><msub>   <mi>x</mi><mn>1</mn></msub></mtd>
		<mtd><msub>   <mi>y</mi><mn>1</mn></msub></mtd>
		<mtd><msub>   <mi>z</mi><mn>1</mn></msub></mtd>
		<mtd><msub>   <mi>t</mi><mn>1</mn></msub></mtd>
		<mtd><msubsup><mi>l</mi><mn>1</mn><mo>(0)</mo></msubsup></mtd>
		<mtd><msubsup><mi>l</mi><mn>1</mn><mo>(1)</mo></msubsup></mtd>
		<mtd><mo>&ctdot;</mo></mtd>
		<mtd><msubsup><mi>l</mi><mn>1</mn><mo>(n)</mo></msubsup></mtd>
	    </mtr>
	    <mtr>
		<mtd><msub>   <mi>x</mi><mn>2</mn></msub></mtd>
		<mtd><msub>   <mi>y</mi><mn>2</mn></msub></mtd>
		<mtd><msub>   <mi>z</mi><mn>2</mn></msub></mtd>
		<mtd><msub>   <mi>t</mi><mn>2</mn></msub></mtd>
		<mtd><msubsup><mi>l</mi><mn>2</mn><mo>(0)</mo></msubsup></mtd>
		<mtd><msubsup><mi>l</mi><mn>2</mn><mo>(1)</mo></msubsup></mtd>
		<mtd><mo>&ctdot;</mo></mtd>
		<mtd><msubsup><mi>l</mi><mn>2</mn><mo>(n)</mo></msubsup></mtd>
	    </mtr>
	    <mtr rowalign="center">
		<mtd columnspan="4"><mspace></mspace></mtd>
		<mtd><mo>&vellip;</mo></mtd>
		<mtd columnspan="2"><mspace></mspace></mtd>
	    </mtr>
	    <mtr>
		<mtd><msub>   <mi>x</mi><mo>n</mo></msub></mtd>
		<mtd><msub>   <mi>y</mi><mo>n</mo></msub></mtd>
		<mtd><msub>   <mi>z</mi><mo>n</mo></msub></mtd>
		<mtd><msub>   <mi>t</mi><mo>n</mo></msub></mtd>
		<mtd><msubsup><mi>l</mi><mo>n</mo><mo>(0)</mo></msubsup></mtd>
		<mtd><msubsup><mi>l</mi><mo>n</mo><mo>(1)</mo></msubsup></mtd>
		<mtd><mo>&ctdot;</mo></mtd>
		<mtd><msubsup><mi>l</mi><mo>n</mo><mo>(n)</mo></msubsup></mtd>
	    </mtr>
	</mtable>
    </math>
    <div class="divcaption">Figure 4. Format of history file when no
    detectors are specified.  The file contains the
    <math xmlns="&mathml;">
	<mfenced><mi>x</mi> <mi>y</mi> <mi>z</mi> <mi>y</mi></mfenced>
    </math>
    position and the path length inside each of the
    <math xmlns="&mathml;"><mi>n</mi></math>
    tissues for all 
    <math xmlns="&mathml;"><mi>N</mi></math>
    of the exiting photons.  All values are stored as 32-bit
    floating point numbers.
    </div>
</div>

<div class="padbox">
    <math xmlns="&mathml;">
	<mtable equalcolumns="false" align="center"
		rowalign="center" columnalign="center">
	    <mtr>
		<mtd><msub>   <mi>d</mi><mn>1</mn></msub></mtd>
		<mtd><msub>   <mi>t</mi><mn>1</mn></msub></mtd>
		<mtd><msubsup><mi>l</mi><mn>1</mn><mo>(0)</mo></msubsup></mtd>
		<mtd><msubsup><mi>l</mi><mn>1</mn><mo>(1)</mo></msubsup></mtd>
		<mtd><mo>&ctdot;</mo></mtd>
		<mtd><msubsup><mi>l</mi><mn>1</mn><mo>(n)</mo></msubsup></mtd>
	    </mtr>
	    <mtr>
		<mtd><msub>   <mi>d</mi><mn>2</mn></msub></mtd>
		<mtd><msub>   <mi>t</mi><mn>2</mn></msub></mtd>
		<mtd><msubsup><mi>l</mi><mn>2</mn><mo>(0)</mo></msubsup></mtd>
		<mtd><msubsup><mi>l</mi><mn>2</mn><mo>(1)</mo></msubsup></mtd>
		<mtd><mo>&ctdot;</mo></mtd>
		<mtd><msubsup><mi>l</mi><mn>2</mn><mo>(n)</mo></msubsup></mtd>
	    </mtr>
	    <mtr rowalign="center">
		<mtd columnspan="2"><mspace></mspace></mtd>
		<mtd><mo>&vellip;</mo></mtd>
		<mtd columnspan="2"><mspace></mspace></mtd>
	    </mtr>
	    <mtr>
		<mtd><msub>   <mi>d</mi><mo>n</mo></msub></mtd>
		<mtd><msub>   <mi>t</mi><mo>n</mo></msub></mtd>
		<mtd><msubsup><mi>l</mi><mo>n</mo><mo>(0)</mo></msubsup></mtd>
		<mtd><msubsup><mi>l</mi><mo>n</mo><mo>(1)</mo></msubsup></mtd>
		<mtd><mo>&ctdot;</mo></mtd>
		<mtd><msubsup><mi>l</mi><mo>n</mo><mo>(n)</mo></msubsup></mtd>
	    </mtr>
	</mtable>
    </math>
    <div class="divcaption">Figure 5. Format of history file when a
    non-zero number of detectors are specified.  The file contains the
    number of the detector that recorded the photon and the path
    length inside each of the <math xmlns="&mathml;"><mi>n</mi></math>
    tissue for all <math xmlns="&mathml;"><mi>N</mi></math> of the
    exiting photons.  Note that it is possible for a single photon to
    be recorded by multiple detectors if the spacing between detectors
    is less than their diameter.  All values are stored as 32-bit
    floating point numbers.
    </div>
</div>

<p>
To compute the fluence at the detector, sum the detected photons and
scale by x.
</p>

<h3 id="sub-two_point_file">1.3 Two-point Files</h3>

<p>
The two-point file (<code>*.2pt</code>) contains the photon density at
every voxel and time gate, as computed by the Monte Carlo simulation.
The data is stored as "double's" (64-bit IEEE floating point
representation) and the voxels are orded the same way as the
segmentation file (see 
<a href="tMCimg.xhtml#sub-segmentation_file">Section&nbsp;1.2</a>).
The real part of all the time gates are written out sequentially,
followed by the imaginary parts (see 
<a href="tMCimg.xhtml#tbl-2pt">Figure&nbsp;6</a>).
</p>

<table id="tbl-2pt" class="noborder">
    <tr>
	<td class="noborder">Time Gate 1</td>
	<td class="noborder">Time Gate 2</td>
	<td class="noborder" align="center">
	    <math xmlns="&mathml;"><mo>&ctdot;</mo></math>
	</td>
	<td class="noborder">Time Gate N</td>
    </tr>
    <tr>
	<td class="noborder">Real Part </td>
	<td class="noborder">Real Part </td>
	<td class="noborder"></td>
	<td class="noborder">Real Part </td>
    </tr>
    <tr>	
	<td colspan="4" align="center" class="noborder">
	    <math xmlns="&mathml;"><mo>&vellip;</mo></math>
	</td>
    </tr>
    <tr>
	<td class="noborder">Time Gate 1 </td>
	<td class="noborder">Time Gate 2 </td>
	<td class="noborder" align="center">
	    <math xmlns="&mathml;"><mo>&ctdot;</mo></math>
	</td>
	<td class="noborder">Time Gate N </td>
    </tr>
    <tr>
	<td class="noborder"> Imaginary Part   </td>
	<td class="noborder"> Imaginary Part   </td>
	<td class="noborder"></td>
	<td class="noborder"> Imaginary Part </td>
    </tr>
    <caption>Figure 6. Ordering of the two-point file.  Each block
    consists of an 
<math xmlns="&mathml;">
    <msub>
	<mi>N</mi>
	<mo>y</mo>
    </msub>
    <mo>&times;</mo>
    <msub>
	<mi>N</mi>
	<mo>x</mo>
    </msub>
    <mo>&times;</mo>
    <msub>
	<mi>N</mi>
	<mo>z</mo>
    </msub>
</math>
    block of doubles.</caption>
</table> 

<p>
Thus, if the simulation is run with two time gates on a
<math xmlns="&mathml;">
    <mfenced>
	<msub>
	    <mi>N</mi>
	    <mo>x</mo>
	</msub>
	<msub>
	    <mi>N</mi>
	    <mo>y</mo>
	</msub>
	<msub>
	    <mi>N</mi>
	    <mo>z</mo>
	</msub>
    </mfenced>
</math>
volume, the first
<math xmlns="&mathml;">
    <msub>
	<mi>N</mi>
	<mo>y</mo>
    </msub>
    <mo>&times;</mo>
    <msub>
	<mi>N</mi>
	<mo>x</mo>
    </msub>
    <mo>&times;</mo>
    <msub>
	<mi>N</mi>
	<mo>z</mo>
    </msub>
</math>
doubles will be the real part of the first time gate, the second 
<math xmlns="&mathml;">
    <msub>
	<mi>N</mi>
	<mo>y</mo>
    </msub>
    <mo>&times;</mo>
    <msub>
	<mi>N</mi>
	<mo>x</mo>
    </msub>
    <mo>&times;</mo>
    <msub>
	<mi>N</mi>
	<mo>z</mo>
    </msub>
</math>
doubles will be the real part of the second time gate, the third 
<math xmlns="&mathml;">
    <msub>
	<mi>N</mi>
	<mo>y</mo>
    </msub>
    <mo>&times;</mo>
    <msub>
	<mi>N</mi>
	<mo>x</mo>
    </msub>
    <mo>&times;</mo>
    <msub>
	<mi>N</mi>
	<mo>z</mo>
    </msub>
</math>
doubles will be the the imaginary part of the first time gate, and the final 
<math xmlns="&mathml;">
    <msub>
	<mi>N</mi>
	<mo>y</mo>
    </msub>
    <mo>&times;</mo>
    <msub>
	<mi>N</mi>
	<mo>x</mo>
    </msub>
    <mo>&times;</mo>
    <msub>
	<mi>N</mi>
	<mo>z</mo>
    </msub>
</math>
doubles will be the the imaginary part of the second time gate.
</p>

<p>
Unscaled fluence for every voxel in the image.  Missing scale factor
is 
<math xmlns="&mathml;">
    <mo>&approx;</mo><mn>1</mn><mo>/</mo><mi>N</mi>
    <mtext>.</mtext>
</math>
</p>

<h2 id="appendix">2 Utility Programs</h2>

<h3 id="cfg2inp">2.1 cfg2inp</h3>

<p>
The <a href="http://www.mathworks.com/">Matlab</a> routines in the <a
href="http://www.nmr.mgh.harvard.edu/~jstott/PMI/">PMI Toolbox</a>
cannot load the <code>*.cfg</code> files used by tMCimg directly.
They use, rather, the older <code>*.inp</code> style of input files
(which, while computer-friendly aren't easily human readable).  The
program <code>cfg2inp</code> (included with <code>tMCimg</code>)
converts from <code>*.cfg</code> files to <code>*.inp</code> files.
If you have other ways of loading the data, there is no need to
compile cfg2inp.
</p>

<h3 id="tmcv">2.2 tmcv</h3>

<p>
<code>tmcv</code> is an OpenGL viewer for <code>*.2pt</code> files, 
distributed as part of the 
<a href="http://www.cims.nyu.edu/~barnett/software/">tfwd</a> package.
While written for the finite difference forward solver in the package,
the input and output formats are identical and it can also be used
with <code>*.2pt</code> files generated by <code>tMCimg</code>.
</p>

</body>
</html>