I. INTRODUCTION
Computational phantoms have been used extensively in radiation protection for Monte Carlo (MC)-based radiation transport and simulation [1]. Accuracy of the phantom plays a critical role in dose calculation [2, 3], and computational phantoms based on color photographic images and containing detailed information of human body have been developed in many countries [4, 6]. However, it is a big problem to deal with the massive voxels accurately [5]. Lots of refinement methods were developed not directly focusing on the original voxels, and sometimes voxels were sampled with interpolation. This may generate some ambiguity that may cause serious problem in whole-body phantom construction. Due to limitation of computer performance, the large organs (e.g. skeleton system) from the image slices cannot be constructed as a whole accurately.
Based on high-resolution sectioned images of a Chinese Visible Human (CVH) dataset, a voxel phantom called Rad-HUMAN (Accurate whole-body computational phantom of Chinese adult female) was established by the FDS Team (www.fds.org.cn).
This CVH data set, which contains 3641 slices with 3286×1586 pixel resolution, was obtained from a Chinese female cadaver. Each voxel size of the data set is 0.15×0.15×0.25 mm for the head and neck and 0.15×0.15×0.5 mm for the rest of the body. The total number of the voxels is about 16.8 billion. Segmentation [6] (Fig. 1) was processed manually and calibrated by MCAM (Multi-Physics Coupling Analysis Modeling Program) developed by the FDS Team [7, 8, 9]. In addition to the fact that the voxel phantom was difficult for Monte Carlo dose calculation, demonstrating its geometry was also laborious [6].
In this work, a flexible sampling method was investigated on the modeling of Rad-HUMAN voxel phantom. Based on the method, the visualization was compared and computational geometries for Monte Carlo dose calculation was discussed.
II. VOXEL SAMPLING METHOD
A. Nearest neighbor sampling
After segmentation, the organs had different RGB colors to be distinguished. In this paper, the data set of voxels assigned with 3D coordinates can be written as
where w, l and h are the width, length and height of the CVH dataset, respectively; and i, j and k are the position subscripts of the voxels.
It is no use obtaining all the voxels of a certain organ to arduously construct its model. In the nearest neighbor sampling, the segmented voxels are appropriately enlarged and sampled from their nearest neighbors.
Change the voxel array from w × l × n to w’ × l’ × n’, the sampling enlargement factor can be expressed as
To be in consistence with the later sampling, the voxel coordinate of an organ can be described as a sparse matrix α
where ∃F(di,j,k)(or
In the nearest neighbor sampling, the sampled organ f(x,y,z) is acquired from an original organ F(X,Y,Z) described as
where symbol "[ ]" stands for rounding off the number for the result. So the resultant coordinate matrix of the sampled organ can be described as
To decrease the demand for computer memory, in present study, the nonzero coordinate (γi,j,k) was stored on disk as file format corresponded with exclusive address for compressed storage.
B. Weighted sampling
As accuracy loss is serious by the nearest sampling, a weighted sampling in the present study was proposed. With the same enlargement factor defined in Eq. (2), a sampling unit can be enclosed by a lattice with a size of Sx × Sy × Sz. Obviously the voxels in the lattice cannot always belong to one organ, hence a weighted factor to judge whether the voxels of the sampling unit can be regard as a voxel after sampling.
After sampling, the relative position of the voxel was changed corresponding with the sampled voxel array (w’ × l’ × n’). That is, the coordinates must be zoomed by a size of (1/Sx) × (1/Sy) × (1/Sz) for the consistence. In addition, described as the coordinates with exclusive addresses, the original voxels of the lattice can be projected to a new address of the sampled storage. This brought the effectiveness for numbering the amounts of the voxels which belong to one organ being projected. So according to a defined weight of the voxels in the lattice, the voxels of a certain organ can be sampled fast. The detailed steps for the weighted sampling are as follows.
Step 1: Obtain coordinate matrix α of an organ referred as Eq. (3).
Step 2: Map the matrix α to a sampling matrix β. Supposing (xα,yα,zα) is an element of α, the element of β becomes
where the mapped address in β can be described as θ with the element obtained from
ω and τ in Eq. (7) are factors to multiply every element of α.
Step 3: Record the weights of certain coordinates which are mapped into the same locations.
Let (x,y,z) be a coordinate after sampling, which are mapped from σ (a region of sampling lattice).
In this step, the coordinates in σ is stored at the same location with δ(x, y, z) recording their amounts.
Step 4: Filter the coordinates from β, with the weights of no less than a weight factor.
Supposing that the weights were defined as the proportional amounts of the sampling unit which is 0.5, the filtered coordinate matrix γ can be selected from β with the condition of
The resultant coordinate matrix γ can also be stored on disk as file format in sequential order.
C. Hybrid voxel sampling
The nearest sampling is a fast way to reduce the data size. However, it is not recommended to use this sampling method on some smaller data because of serious accuracy losses. In this paper, we combine it with the weighted sampling as a two-step hybrid sampling algorithm. In the hybrid sampling, the sampling unit is divided into smaller lattices for the nearest sampling, and then the smaller lattice can be assighed with weights for the weighted sampling.
To that end, the enlargement factor of (Sx, Sy, Sz) was divided into (Nx, Ny, Nz) and (Px, Py, Pz) for the nearest sampling and the proportional sampling, respectively.
Using an enlargement factor of (Nx, Ny, Nz) for the nearest sampling, the original coordinate matrix can be obtained according to Eq. (5) where
In order to combine the nearest sampling with the weighted sampling for the hybrid sampling, the factor in Eqs. (6) and (7) for the mapping and locating becomes
and τ = (1/Px, w/(Py × Sx), w × l / (Pz × Sx × Sy)), respectively.
The detailed steps to implement the hybrid sampling algorithm are shown in Fig. 2.
III. RESULTS AND DISCUSSION
A. Application in construction of polygon-surface phantom
The organs or tissues of interest (e.g., lungs, liver, skin etc.) from the original tomographic photography were identified by assigning every pixel with an identification numbers. All these numbers can be stored sequentially as a pixel matrix that can be used to extract equivalent matrix, from which the polygon-surface model can be constructed by Marching Cubes [10].
The equivalent matrix can be acquired by the last step of the sampling method. For example in weighted sampling, on purpose of generating the matrix γ, the coordinates satisfying the condition of δ(x,y,z)>0.5 × Sx × Sy × Sz are set as an identification number, while the coordinates not satisfying that condition are set as another identification number. The changed matrix γ can then be used as the equivalent matrix using VTK toolkit to generate the polygon-surface model.
Preparatorily, these sampling methods were processed with the just one CPU of the 3.10 GHz Intel CoreTM 2Quad Processor i5-2400 64 bit operating system. To compare with two special cases of the hybrid sampling, the construction of heart organ containing 265 slices was taken for example. Defining the proportional amounts as the weights, Table 1 presents relevant sampling parameters between the nearest and the weighted sampling methods which are compared by the visualization in Figs. 3(a) and 3(b).
| Comparison | Nearest | Weighted |
|---|---|---|
| Enlargement factor proportional | 2.0 2.0 2.0 | 2.0 2.0 2.0 |
| Weight | 0.5 | |
| Voxel size (mm3) | 0.3×0.3×1.0 | 0.3×0.3×1.0 |
| Times (s) | 19 | 268 |
| Number of triangular faces | 589 005 | 571 133 |
| Holes | 1516 | 1274 |
| Model size (MB) | 29.2 | 28.5 |
When the voxel increased with a factor of 2.0 in x, y and z directions, the time expenditure by the nearest sampling was about 7.09% of that by the proportional sampling. The number of the holes and model size produced by the weighted sampling was 84% and 95% of the one produced by the nearest sampling, respectively. By the weighted sampling, the incorrect model incurred by segmentation can be repaired (Fig. 3(d)).
Parameters between the two methods (Table 1) are mainly dependent on computer performance and data set resolution etc. For precision adjustment, the proportional weight factor can be used.
B. Improvement on the sampling and analysis
The nearest and the weighted sampling method in present paper are implied as two special cases of the hybrid sampling. To discuss their combined criteria, extra experiment was taken with a large organ (e.g. skeleton system, include marrow, pelvis, cartilage etc.). With the same enlargement size of 4.0 × 4.0 × 4.0 and a proportional weight of 0.5 for the proportional and the hybrid method, modeling parameters between the three samplings are given in Table 2 and the visualization comparison is illustrated in Fig. 4.
| Comparison | Nearest | Weighted | Hybrid |
|---|---|---|---|
| Enlargement factor | 4.0 4.0 4.0 | 4.0 4.0 4.0 | 4.0 4.0 4.0 |
| Proportional weight | 0.5 | 0.5 | |
| Voxel size (mm3) | 0.6 × 0.6 × 2.0 | 0.6 × 0.6 × 2.0 | 0.6 × 0.6 × 2.0 |
| Times (s) | 142 | 2374 | 431 |
| Number of triangular faces | 2 504 154 | 2 512 610 | 2 623 755 |
| Holes | 31 635 | 4245 | 5200 |
| Model size (MB) | 137.62 | 124.74 | 130.39 |
In the hybrid sampling, the enlargement size of 4.0 × 4.0 × 4.0 was divided to 2.0 × 2.0 × 2.0 for both the nearest and weighted. The time expenditure by the hybrid sampling was about 18.15% of that by proportional. In addition, the hole number of the model produced by the hybrid was about 16.43 percent of that produced by the nearest. It is advantageous to create the model accurately and quickly for visualization [11, 12].The visualization of the model constructed by the three sampling is shown in Fig. 4.
The nearest sampling may cause severe distortion (Fig. 4(a)) when the sampling enlargement factor becomes larger. The hybrid sampling, therefore, is an appropriate method to make up for such deficiency and enhance the process. It is a flexible sampling method to be evolved from the one to another with a distributed enlargement factor. The combined criterion depends on the sampling ratio, dataset size, computer performance etc. Generally, considering the visual impact, the enlargement size allocated for the nearest sampling shall not be larger than 2.0 × 2.0 × 2.0. This can largely reduce the data size for fast processing with current configuration of a personal computer. For the demanding of more accuracy of the model (Figs. 4(b) and 4(c)), the enlargement factor allocated for the weighted sampling should be larger. Because there are many organs contained in the whole phantom, the accepted time to visualize and check the geometry of the phantom can be satisfied by the hybrid sampling.
From Table 2, the percent of holes to faces of the nearest, the weighted and the hybrid is 1.26%, 0.17% and 0.20%, respectively. Compared with the visualization of Fig. 4, the percent of about 0.20% with the model can be as a reference acceptance criteria for better visualization.
By the hybrid sampling, the voxels are enlarged by a factor of 4.0 × 4.0 × 4.0, the whole skeleton system can be shown in Fig. 5.
C. Application in Monte Carlo computational geometries
For dose calculation in lattice geometry, the methods presented in this paper can be applied [13]. In this case for Monte Carlo simulation, the pixel matrix of the whole phantom should be generated.
The weights definition can be designed by the material, density or other biological property. Next, accumulate the weights from the resultant coordinate matrix (e.g. γ), locate the voxel in the pixel matrix satisfying the weighted factor, and assign it with the identification number. Finally, the pixel matrixes of body can be transferred to lattice representation geometry in MCNP codes [14, 15].
However, in addition to the lattice geometry [16, 17], the polygon-surface model, which is constructed by Marching Cubes [18], can hardly avoid surface intersections and holes in quality. For dose calculations with deformable phantoms [19, 20], repairing/improvement (e.g. filling, relaxation) is needed, so as to confirm no imperfection or interference, and to transfer it into those accepted by MC codes (e.g. MCNP, Geant4, etc.).
IV. CONCLUSION
In this paper, a hybrid voxel sampling algorithm, which merged the nearest sampling and a weighted sampling, was put forward. Through discussions of the model construction, the three sampling methods are able to deal with mass voxels. For accurate sampling, the weighted and the hybrid sampling are effective. As a flexible sampling, the hybrid sampling is able to show high-performance especially on construction of large organs. It is a valid sampling method on the visualization of the (MC) geometries of human voxel phantom.
Furthermore, based on these methods, the geometries for Monte Carlo dose calculation were then analyzed. The phantom of polygon-surface geometry or lattice representation geometry is able to be transferred to MC codes for dose calculations. Finally, by the proposed sampling, more accurate CAD (Computer Aided Design) (STEP etc.) model can be constructed in the field of medical diagnosis and other scientific researches.
Comp Graph, ACM Siggraph’87
,