Visualization of boundaries in CT volumetric data sets using dynamic M−|∇f| histogram

Computers in Biology and Medicine ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Q1 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

Contents lists available at ScienceDirect

Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm

Visualization of boundaries in CT volumetric data sets using dynamic M
j ∇f j histogram Lu Li a, Hu Peng b,n, Xun Chen b, Juan Cheng b, Dayong Gao c a

Electronic Science and Technology, University of Science and Technology of China, Anhui, China School of Medical Engineering, Hefei University of Technology, Anhui, China c Department of Mechanical Engineering and Department of Bioengineering University of Washington, Seattle, WA 98195, USA b

art ic l e i nf o

a b s t r a c t

Article history: Received 27 August 2015 Accepted 28 October 2015

Direct volume rendering is widely used for three-dimensional medical data visualization such as computed tomography and magnetic resonance imaging. Distinct visualization of boundaries is able to provide valuable and insightful information in many medical applications. However, it is conventionally challenging to detect boundaries reliably due to limitations of the transfer function design. Meanwhile, the interactive strategy is complicated for new users or even experts. In this paper, we build a generalized boundary model contaminated by noise and prove boundary middle value (M) has a good statistical property. Based on the model we propose a user-friendly strategy for the boundary extraction and transfer function design, using M, boundary height ðΔhÞ, and gradient magnitude ðj ∇f j Þ. In fact, it is a dynamic iterative process. First, potential boundaries are sorted orderly from high to low according to the value of their height. Then, users iteratively extract the boundary with the highest value of Δh in a newly defined domain, where different boundaries are transformed to disjoint vertical bars using M
j ∇f j histogram. In this case, the chance of misclassification among different boundaries decreases. & 2015 Published by Elsevier Ltd.

Keywords: Direct volume rendering Multi-dimensional transfer function Volume visualization 3-D CT volume Boundary extraction

1. Introduction Direct volume rendering is a widely used visualization technique to demonstrate the internal structures of volume data sets. To generate a meaningful visualization, transfer functions (TFs) should be designed to map data properties (e.g. scalar value, gradient) to optical properties (e.g. color, opacity). Color is used to generate a visual distinction between different data properties. Opacity determines the visual degree for each voxel of the volume data. Since data properties may be changing with different applications, designing an ad hoc TF is quite necessary to make the structures of interest visible and distinguishable. The most commonly used 1-dimensional (1-D) transfer function maps scalar value to color and opacity based on the histogram of scalar value occurrences. Users need to select the materials of interest and set color and opacity by trial-and-error method. The visualization of boundaries of 3-dimensional (3-D) computed tomography (CT) data has a wide range of applications and great significance for disease diagnosis and screening [1–3]. In n

Corresponding author. E-mail addresses: [email protected] (L. Li), [email protected] (H. Peng), [email protected] (X. Chen), [email protected] (J. Cheng), [email protected] (D. Gao).

order to highlight boundaries, gradient magnitude ðj ∇f j Þ has been adopted as a data property [4,5]. Adding j ∇f j as an independent dimension to the 1-D TF of scalar value, the 2-dimensional (2-D) f
j ∇f j histogram [6] has been proposed to facilitate boundary extractions. Boundaries in f
j ∇f j histogram are represented as arches, enabling users to extract different boundaries and design TFs for individual boundaries. Nevertheless, it is difficult to extract arches, especially considering the overlapping among arches. Low/high (LH) histogram [7], therefore, has been introduced to overcome the drawback by using the low and high intensities of materials around the boundary as two key data properties. In LH histogram, boundaries are expressed as points (without noise), lines and regions (with noise and bias), making the implementation of boundary extractions much easier. As is known to all, boundaries represented in 2-D histograms are often extracted interactively, while the strategies are too complicated for new users or even experts to follow. For instance, in f
j ∇f j histogram users have to precisely pick out each arch by widgets. In LH histogram, users have to pick out each targeted region manually. Under the situation that the noise and overlapping is large, arches with smaller size are apt to be hidden by larger ones in f
j ∇f j histogram and different targeted regions tend to merge in LH histogram. Thus, it is a tedious or even impossible task for users via trial-and-error process.

http://dx.doi.org/10.1016/j.compbiomed.2015.10.018 0010-4825/& 2015 Published by Elsevier Ltd.

Please cite this article as: L. Li, et al., Visualization of boundaries in CT volumetric data sets using dynamic , Comput. Biol. Med. (2015), http://dx.doi.org/10.1016/j.compbiomed.2015.10.018i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96

L. Li et al. / Computers in Biology and Medicine ∎ (∎∎∎∎) ∎∎∎–∎∎∎

2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

Misclassification among different boundaries, to some extent, can be avoided by carefully choosing arches and regions in f
j ∇f j and LH histograms respectively. However, as the number of boundaries increases, the effect of misclassification will accumulate gradually. This will likely have a negative impact on the subsequent arches/regions selection. When the effect reaches a certain level, it is impossible for users to pick out the remaining boundaries with high quality. They have to extract boundaries using the original histogram from the very beginning. In this paper, on one hand, we build a generalized boundary model contaminated by noise. On the other hand, based on the model we propose a novel multidimensional TF design method using boundary middle value (M), boundary height ðΔhÞ and gradient magnitude ðj ∇f j Þ as data properties. Unlike traditional histograms presenting all boundary information once, the proposed method employs dynamic M
j ∇f j histogram and presents only one or a few boundaries at each step. A simple and iterative strategy of boundary extraction is also developed. We first sort different boundaries according to the values of their heights ðΔhÞ from high to low. Users control the Δh value until a vertical bar occurs in M
j ∇f j histogram, then pick out the boundary represented by the bar in M
j ∇f j histogram. Boundaries are extracted one by one until the Δh value reaches 0. The misclassification among different boundaries can be reduced with the help of boundary ordering and the one-by-one extraction strategy. Besides, region elimination and region growing are further adopted to enhance the quality of rendering. There are two main contributions. First, we establish a twomaterial boundary model contaminated by noise and prove that M has a good statistical property in this model. Second, we propose a novel 3-D TF in a newly defined domain, using boundary middle value (M), boundary height ðΔhÞ and gradient magnitude ðj ∇f j Þ. Based on the dynamic M
j ∇f j histogram, we propose an iterative boundary extraction method. In the following section, we describe a variety of related methods. In Section 3, we build a generalized boundary model contaminated by noise. Based on the model we describe the construction and properties of the dynamic M
j ∇f j histogram, as well as refining techniques. In Section 4, we show several examples using our method and the traditional f
j ∇f j and LH histograms, and make a comparison. In Section 5, we present the conclusion and future work.

2. Background In this section, we present a more detailed description of some important work on the design of transfer functions. As the aim of this paper is to visualize boundaries, we give a special attention to the relevant methods in literature. Transfer function, which maps data properties to optical properties, is of crucial importance to volume visualizations. According to the type of the extracted data properties, designs of TFs can be roughly divided into two categories, data-driven and graph-driven. 2.1. Data-driven TFs based on boundaries Data-driven methods mainly focus on statistical features by means of histograms. The most widely used 1-D TF aims to map scalar value to opacity properties. However, it is inappropriate to employ the 1-D TF for the boundary visualization because the scalar values of the boundary vary greatly. In order to highlight boundaries, gradient magnitude ðj ∇f j Þ is added as a data property. Levoy [4] and Drebin et al. [5] put forward a 1-D TF to display different boundaries using gradient features. Later, Kindlmann and Durkin [6] proposed a 2-D f
j ∇f j histogram to enhance

boundaries with the horizontal axis representing scalar values and the vertical axis representing gradient magnitude. The arches in the histogram represent boundaries among different materials. Although f
j ∇f j histogram improves the selection of the boundaries, there still exists a drawback. The overlapping regions caused by the intersections of the arches lead to ambiguities in the classification of boundaries. Kniss [8,9] added second derivative into f
j ∇f j histogram to select only the peaks of the arches to obtain much better rendering results. Unlike the aforementioned histograms using statistical features directly, LH histogram is based on the concept of the boundary model proposed by Nickoloff and Riley [10] and treats boundaries as the transition of two idealized homogeneous regions. Lum and Ma [11] proposed the early concept of LH. LH means the high intensity and low intensity of the homogeneous regions. Arches in f
j ∇f j histogram are transferred to straight lines in a 1-D histogram with the x-axis representing L and H for an easily recognition. Sereda et al. [7] compacted the arches to points and regions using LH histogram for a more convenient extraction of boundaries. Compared with the f
j ∇f j histogram, LH histogram represents boundaries more compactly. Serlie et al. [12] derived the relationship between f
j ∇f j and LH histogram in a strict mathematic form based on the boundary model. One drawback of LH histogram is the time to establish the histogram. Praßni et al. [13] proposed a new way to efficiently construct LH histogram. Region selecting is a key issue to achieve automatic extraction of boundaries using LH histogram. Sereda et al. [14] used hierarchical clustering to automatically classify LH histogram to regions. Nguyen et al. [15] applied mean shift clustering to LH histogram and then used hierarchical clustering to group similar voxels. 2.2. Data-driven TFs for other features Apart from visualizing boundaries using gradient information, curvature information [16,17] was also used to distinguish objects based on their shapes. Other statistical data properties such as mean value [18] and texture [19,20], were used in the designs of TFs to highlight the structures of interest. These local statistical properties may produce excellent results in specific applications. 2.3. Graph-driven TFs Graph-driven methods mainly aim to visualize the local structures, which depend on topological analysis, such as skeletonization, region growing and level-set. As a key concept for graph-driven methods, connectivity is the foundation of the skeletonization and region growing. Takahashi et al. [21] extracted topological structures of volume data sets, leading to a graph called a volume skeleton tree consisting of volumetric critical points and their connectivity. Ji et al. [22] segmented the 3-D teeth in cone beam computed tomography data by level-set. Volume data segmentation proposed by Gooya et al. [23] is based on manifold distance metrics with the features of scalar value, gradient and probabilistic measures. Our proposed method mainly employs the gradient-based features and boundary model. Compared with the construction of LH histogram and f
j ∇f j histogram, the key improvement is that we build a generalized boundary model contaminated by noise and prove M has a good statistical property. Graph-driven method is employed to optimize the rendering result using region growing and region elimination based on local connectivity to reduce the misclassification, while LH histogram only uses region growing to fulfill the extracted boundaries without considering the misclassification between different boundaries.

Please cite this article as: L. Li, et al., Visualization of boundaries in CT volumetric data sets using dynamic , Comput. Biol. Med. (2015), http://dx.doi.org/10.1016/j.compbiomed.2015.10.018i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

L. Li et al. / Computers in Biology and Medicine ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 Q3 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

2.4. User interface of TF designs Almost each TF design tends to provide a user-friendly interface for convenient interaction. Histogram methods of data-driven TF design provide various manipulation widgets for users based on objective statistical features [24–26]. Though these sophisticated widgets have been proposed, these widgets often do not provide simple strategies. To achieve precise operation, users should develop professional knowledge about the widgets and volume rendering. Recently, TFs with semantics [27,28] and sketch-based user interface [29,30] have been proposed to integrate intuitions into the TF design. In this paper, we propose an interactive TF design method employing dynamic M
j ∇f j histogram. Users extract the vertical bars one by one with Δh from high to low. This is a convenience interactive strategy even for non-experts with little domain knowledge.

3. Methods 3.1. Classic boundary model The classic two-material boundary (f) with L representing the low intensity and H representing the high intensity is modeled by a scaled step function ðuðx; L; HÞÞ convolved with 1-D Gaussian edge-spread-function (g) [12]. The edge locates at x ¼0. The model is shown to be excellent to describe the two-material boundary in CT data sets [31]. We introduce two new random variables M 9 ðL þ HÞ=2, Δh 9 H
L. In the classic boundary model, M and Δh are constants. For a given boundary, M represents the middle value of the boundary and Δh represents the height of the boundary (Fig. 1): uðx; L; HÞ ¼ L þ ðH
LÞ  uðxÞ fðx; σ Þ ¼ uðx; L; HÞngðx; σ Þ ¼ L þ ðH
LÞ

ð1Þ

1  pffiffiffiffiffiffi σ 2π with uðxÞ ¼

(

Z



x

exp
1

0

xo0

1

xZ0

3


η2 2σ 2



dη ¼ L þðH
LÞ 



  1 1 x pffiffiffi þ erf 2 2 σ 2

ð3Þ

  1
x2 gðx; σ Þ ¼ pffiffiffiffiffiffi exp 2 2σ σ 2π 2 erf ðxÞ ¼ pffiffiffiffi

π

Z

x 0

ð2Þ

ð4Þ

expð
η2 Þ dη

ð5Þ

3.2. Generalized boundary model The classic boundary model does not consider the interference caused by noise. The major component of the noise of CT is shot noise which refers to the uncertainty of the intensity of X-ray due to its quantum origin [32]. Shot noise can be modeled by a Poisson distribution, which except at low intensity levels approximates a Gaussian distribution [33]. In general the noise of CT is modeled by a Gaussian distribution [34,35]. As the two-material boundary model contains the material and the boundary, we discuss the effects on the material and the boundary respectively when the model is contaminated by Gaussian noise. 3.2.1. Gaussian noise of the materials We first consider the Gaussian noise of the materials, L and H obey normal distributions (Eq. (6)). Denote L  NðL0 ; σ l Þ, and H  NðH 0 ; σ h Þ respectively. In Appendix A(a) and (b), prove qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qwe ffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi that M  NððL0 þ H 0 Þ=2; 2 

σ 2l þ σ 2h Þ, Δh  NððH 0
L0 Þ;

1
ðx
μÞ2 pðxÞ ¼ pffiffiffiffiffiffi exp 2σ 2 σ 2π

!

where p(x) is the probability density at x, distribution, σ is its standard deviation.

σ 2l þ σ 2h Þ:

ð6Þ

μ is the mean of the

Fig. 1. (a) The boundary modeled by a scaled step function. (b) 1-D Gaussian edge-spread-functions with two different σ. (c) The blurred boundary modeled by the scaled step function convolved with the 1-D Gaussian edge-spread-function. (d) For a given boundary, M and Δh represent the middle value and the height of the boundary respectively.

Please cite this article as: L. Li, et al., Visualization of boundaries in CT volumetric data sets using dynamic , Comput. Biol. Med. (2015), http://dx.doi.org/10.1016/j.compbiomed.2015.10.018i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

L. Li et al. / Computers in Biology and Medicine ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

As M and Δh obey the normal distribution, the expectation (E) of each random variable is a constant. EðMÞ ¼

ðL0 þ H 0 Þ 2

EðΔhÞ ¼ H 0
L0

ð7Þ ð8Þ

3.2.2. Line search Line search is the backbone of many optimizations, which is defined as: ( x k þ 1 ¼ x k þ λ dk ð9Þ Minf ðxk þ λdk Þs:t: λ A ½0; þ 1Þ The goal is to find the appropriate λ to satisfy Eq. (9). In our boundary model, L and H can be regarded as the local maxima and minima along the boundary profile. Thus line search with the firstorder derivative [36] is adopted to find L and H. The line search stops at x, where ∇f ðx
εdÞn∇f ðx þ εdÞ r 0, where ε is a small real number. As L and H are not constants, it is unreliable to determine the initial point only by j ∇f j 4 0. The interference caused by noise should be eliminated. Assume in material L the intensities of two adjacent voxels are Ll1 and Ll2 and the distance between the two voxels is 1, we can get j ∇f l j ¼ j Ll1
Ll2 j . As L  NðL0 ; σ l Þ, the probability distribution of j ∇f l j can be calculated using Eq. (10) (Appendix pffiffiffi A(b) and (d)). Since the probability of j ∇f l j in the range of ½0; 2 2σ l  is 0.954, j ∇f j of the initial point should satisfy the condition in Eq. (11): ! 2
ðj ∇f l j Þ2 ð10Þ pðj ∇f l j Þ ¼ pffiffiffiffiffiffiexp j ∇f l j Z0 4σ 2l σ l 4π pffiffiffi pffiffiffi j ∇f j 4 maxfσ l 2 2; σ h 2 2g

ð11Þ

3.2.3. Gaussian noise of the boundary Then we consider the Gaussian noise to the boundary. Since line search with the first-order derivative is adopted as the searching method, we first calculate the first-order derivative of

the noise, which can be represented by the first order difference of two random variables that obey the normal distribution. The gradient magnitude of the noise ðj ∇f noise j Þ obeys the probability distribution in Eq. (12) (Appendix A(d)). If the noise takes place on the boundary at location x0, the line search is interfered in this Þ j j , the line search is not disturbed manner: If j ∇f noise j r j dðf dx x ¼ x0 Þ by the noise (Fig. 2(c) and (d)). If j ∇f noise j 4 j dðf j j , the line dx x ¼ x0 search stops at point x0 (Fig. 2(a) and (b)). Assume the stop points are x ¼ x1 when searching L and x ¼ x2 when searching H, M ¼ ðf ðx1 Þ þ2 f ðx2 ÞÞ: ! 2
ðj ∇f noise j Þ2 ð12Þ j ∇f noise j Z0 pðj ∇f noise j Þ ¼ pffiffiffiffiffiffi exp 2σ 2n σ n 2π Based on the interference manner, we assume that the line Þ j j for a straightforward calcusearch stops at j ∇f noise j ¼ j dðf dx x ¼ x0 lation: 8   < j ∇f j ¼ p1ffiffiffiffiffiðH
LÞ exp
x22s 2σ σ 2π ð13Þ : j ∇f j ¼ j ∇f noise j where xs is the stopping point along the boundary profile when searching the local maxima and minima. It is obvious that pðxs Þ is symmetrical about xs ¼ 0. When xs r 0, 2 ffiðH
LÞexpð
x2s Þ is monotonic increasing. Using the j ∇f j ¼ σ p1ffiffiffiffi 2σ 2π theorem in Appendix B, We can get:   dðj ∇f j Þ  xs r 0 ð14Þ pðxs Þ ¼ pðj ∇f j Þ dxs  where pðj ∇f j Þ is from Eq. (12). Eq. (13) is symmetrical about xs ¼ 0. Then we can get the probability distribution of the stopping point xs in Eq. (15). Fig. 3(a) shows the functional digraph of the following equation: pðxs Þ ¼ pðj ∇f j Þj

dðj ∇f j Þ j dxs

ð∇f Þ2 x2s ¼ ξ j xs j exp

2σ 2n 2σ 2

! ð15Þ

Þ Fig. 2. (a) and (b) Noise takes place at x ¼ x0 and j ∇f noise j 4 j dðf j j . The line search of H stops at x ¼ x0 when starting at T. (c) and (d) Noise takes place at x ¼ x1 and dx x ¼ x0 Þ j ∇f noise j r j dðf j j . The line search of H is not interfered when starting at T. The interference manner is the same for the line search of L. dx x ¼ x1

Please cite this article as: L. Li, et al., Visualization of boundaries in CT volumetric data sets using dynamic , Comput. Biol. Med. (2015), http://dx.doi.org/10.1016/j.compbiomed.2015.10.018i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

L. Li et al. / Computers in Biology and Medicine ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

5

Fig. 3. Let the given boundary be f ðL; HÞ. (a) The functional digraphs of xs in Eqs. (15) and (16) when ξ ¼ 1. The approximation does not affect our reasoning process. (b) The probability distribution of y. The probability distribution is symmetrical about y ¼ L þ2 H. (c) The probability distribution of M. The probability distribution is symmetrical about M ¼ L þ2 H.

where ξ ensures that the whole integral of the probability distribution equals to 1. Þ2 is symmetrical about xs ¼0 and its influence on As the term ð∇f 2σ 2n pðxs Þ is minor (Fig. 3(a)). We discard the term for a straightforward calculation:   x2 pðxs Þ ¼ ξ j xs j exp
s 2 ð16Þ 2σ Based on Eq. (16) and f ðxs Þ is monotonically increasing, we can get the probability distribution of y ¼ f ðxs Þ (Fig. 3(b)). Without loss of generality, we assume H
L ¼ 1: 8 R   2 xs > < pðyÞ ¼ d ξ j η j exp
2ησ 2 dη
1 dy ð17Þ >
1 : x ¼ σ pffiffiffi 2erf ð2y
1Þ s M can be indicated as M ¼ y1 þ2 y2 , where y1 and y2 obey the probability distribution in Eq. (17). The probability distribution of M can be calculated using Eq. (18). As the form of Eq. (18) is complicated, we use numerical integration to get the functional digraph of the probability distribution of M (Fig. 3(c)). We can conclude that M exists in the range of ½L þ2 H
ε; L þ2 H þ ε with a large potential, where ε is a small real number: Z þ1 pðMÞ ¼ pðyÞpð2M
yÞ dy ð18Þ
1

For a given boundary f ðL; HÞ, we calculate the expectation of M. As the probability distribution of M is symmetrical about M ¼ L þ2 H: EðMÞ ¼

ðL þ HÞ 2

ð19Þ

Combining Eqs. (7), (8) and (19), we can summarize under the situation that the noise of the material and the noise of the boundary are added simultaneously: ðL0 þH 0 Þ 2

ð20Þ

EðΔhÞ o H 0
L0

ð21Þ

EðMÞ ¼

3.3. Dynamic M
j ∇f j histogram For a given boundary f ðL; HÞ, we have proved that M exists in the range ½L þ2 H
ε; L þ2 H þ ε with a large potential, where ε is a small real number. However, it is not suitable to extract the boundaries though cluster approaches immediately because of the existence of misclassification among different boundaries. Let the set of voxels of a extracted boundary be B. The misclassification is constituted by two aspects: one is the false identifications of other boundary elements (set T), the other is the miss of true boundary elements (set N). The true boundary set is B
T þ N. Considering the fact that if the boundary is contaminated by noise, the height detected by the line search is less than its real value. Based on this

fact, the boundary with Δh1 interferes the boundaries with the height less than Δh1 , while the boundaries with the height less than Δh1 cannot interfere the boundary with Δh1 . That is to say, the boundary with the maximum Δh can interfere other boundaries, while cannot be interfered by others. In other words, set T of the boundary with the maximum Δh is an empty set. Thus we only need to consider set B þN for the boundary with the maximum Δh. In addition, if the boundary with the maximum Δh can be extracted precisely, we can establish an iterative boundary extraction method according to Δh from high to low without concerning set T. Since Δh is a dynamic variable which cannot be adopted as a coordinate of a histogram, it is not intuitive to determine the value of Δh and the range of M via the 1-D histogram of M. Thus we introduce j ∇f j as a coordinate besides M, because j ∇f j is an inherent property to represent a boundary. We establish a dynamic 2-D M
j ∇f j histogram with the horizontal coordinate representing M and the vertical coordinate representing j ∇f j . The form of a given boundary is a vertical bar in M
j ∇f j histogram. Δh is used to control the histogram. With Δh from high to low, boundaries emerge in M
j ∇f j histogram orderly as the height of the boundaries. Considering the misclassification of the boundary with the maximum Δh extracted by M
j ∇f j histogram, the aim is to rehabilitate set N þ B which refers to the true boundary set. As the whole boundary elements form an integral 3-D surface, we take set B as seed points and use region growing on the conjugate direction of ∇f to enable B towards N þ B. In region growing, some labeled voxels are chosen as initial seed points. From these seed voxels, the region is growing by iteratively adding the neighbors of already labeled voxels. In each iteration, a neighbor is chosen when it is similar to the already labeled voxel. In general, similarity measure is introduced to measure the similarity of the voxels. As our region growing aims to extract all the boundary elements of the current boundary, we introduce two similarity measures to enable the region growing to take place on the current boundary. These two measures depend on the assumption that the whole boundary elements form an integral 3-D surface. One similarity function measures the scalar value in Eq. (22). The other measures the gradient magnitude in Eq. (23). Fig. 4 shows an example of our region growing method. The blue boundary forms an integral after the region growing:   f ðxÞ
f ðx1 Þ  oδ ð22Þ   f ðxÞ where f(x) is the scalar value of the labeled voxel, f ðx1 Þ is the scalar value of the neighbor, δ is a small real number near to 0:    j ∇f ðxÞj
j ∇f ðx1 Þj   oε ð23Þ   j ∇f ðxÞj where j ∇f ðxÞj is the gradient magnitude of the labeled voxel, j ∇

Please cite this article as: L. Li, et al., Visualization of boundaries in CT volumetric data sets using dynamic , Comput. Biol. Med. (2015), http://dx.doi.org/10.1016/j.compbiomed.2015.10.018i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

6

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 Q5 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

L. Li et al. / Computers in Biology and Medicine ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Fig. 4. An instance of region growing on the conjugate direction of ∇f . (a) The blue boundary is incomplete after the boundary extraction. (b) The effect after region growing when δ ¼ 0:1 and ε ¼ 0:1. The boundary is still incomplete. (c) The effect after region growing when δ ¼ 0:2 and ε ¼ 0:2. The boundary is complete. (d) The effect after region growing when δ ¼ 0:3 and ε ¼ 0:3. The boundary is complete. In practice, we recommend δ o 0:2 and ε o 0:2. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

Fig. 5. An instance of region elimination on the conjugate direction of ∇f . There are 8 pixels in the neighborhood. A threshold is set to count the neighbors of each boundary element. (a) The red rectangular is contaminated by isolated regions after boundary extraction. (b) The effect after region elimination when the threshold is 3. Some isolated regions exist. (c) The effect after region elimination when the threshold is 4. Some isolated regions exist. (d) The effect after region elimination when the threshold is 5. There exists no isolated region. In practice, we recommend the threshold to be larger than 10 when there are 26 voxels in the neighborhood in the 3-D volume. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

f ðx1 Þj is the gradient magnitude of the neighbor, ε is a small real number near to 0. As N is an abstract set, there exists no approach to precisely determine B after region growing equals to B þ N. Considering the worst case: a slight difference exists between B and B þ N for all boundaries, cumulative error occurs in the form of isolate regions since iterative method is adopted. As the isolate regions can be directly eliminated by image morphology, we use region elimination on the conjugate direction of ∇f (Fig. 5) to remove the cumulative error. Fig. 6 gives an example of the boundary extraction based on the dynamic M
j ∇f j histogram. With Δh from high to low, different bars emerge in the M
j ∇f j histogram orderly as their heights. Different bars represent different boundaries. It is convenient for users to obtain the horizontal coordinate range of each vertical bar via clustering methods or even the intuition of users when the histogram with scale. M
j ∇f j histogram can be recognized as a 1-D histogram more than a 2-D one.

3.4. TF design based on the dynamic M
j ∇f j histogram We establish a 3-D TF based on Δh and M
j ∇f j histogram. While controlling Δh from high to low and extracting bars in M
j ∇f j histogram, users design TFs by assigning the extracted boundaries color and opacity. Color is used for the visual grouping. Users can interactively design the color of the boundaries by their intuitive feelings. Since the TF of the opacity is abstract for non-experts, we design the mathematical expression of TF for the opacity with a parameter controlled by users interactively (Eq. (22)). As j ∇f j obtains its maxima at the edge point in our boundary model, the expression is designed to emphasize the voxels close to the edge point: TFðαÞ ¼ where



j ∇f j j ∇f j max

λ

λZ1

ð24Þ

α is the opacity, λ is controlled by the users.

Please cite this article as: L. Li, et al., Visualization of boundaries in CT volumetric data sets using dynamic , Comput. Biol. Med. (2015), http://dx.doi.org/10.1016/j.compbiomed.2015.10.018i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

L. Li et al. / Computers in Biology and Medicine ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

7

Fig. 6. An artificial image with circle R1, circle R2 and rectangular S1, which is contaminated by Gaussian noise. (a) M
j ∇f j histogram when Δh 4120. (b) The blue circle is the result of the extraction of the bar of (a). (c) M
j ∇f j histogram when Δh 475. (d) The green circle is the result of the extraction of the bar of (c). (e) M
j ∇f j histogram when Δh 4 55. (f) The red rectangular is the result of the extraction of the bar of (e). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

4. Results We demonstrate our method in two steps. In the first step, we compare our interactive method with f
j ∇f j histogram and LH histogram. In the second step, our method is tested with medical data sets to demonstrate rendering quality. The data sets are obtained from the volume library http://schorsch.efi.fh-nuernberg. de/data/volume/. 4.1. User interface Fig. 7 shows the user interface of our interactive TF design. There are two steps to design the TF of a 3-D data set. The first step is finding the appropriate value of Δh according to M
j ∇f j histogram. Δh decreases from its maximum until one bar appears in M
j ∇f j histogram. We design a widget that can automatically determine the appropriate Δh for each boundary based on the half-interval search algorithm. The second step is extracting the boundary represented by the bar and assigning color and opacity to the boundary. The TF of all

boundaries is designed in an iterative way until Δh descends to 0, which is straightforward to be manipulated by users. The software package can be obtained at https://github.com/lilu630/Dynamicmiddle-gradient-histogram.git. 4.2. Interactive method Fig. 8 shows the five iterative steps to extract all the boundaries in the data set of the tooth sized 256n256n161. With Δh from high to low, the vertical bars are extracted one by one without mutual interference. On one hand, the extraction of the bar can be easily understood and manipulated by non-experts. On the other hand, the process of extracting all boundaries is an iteration of extracting one boundary no matter how many boundaries are in the data set, which means non-experts can design the TF of all boundaries without concerning the number of the boundaries. Fig. 9 shows the comparison of the interactive methods of f
j ∇f j , LH and M
j ∇f j histogram. In f
j ∇f j histogram (Fig. 9(a)) the boundaries are represented as arches. It is difficult to determine

Please cite this article as: L. Li, et al., Visualization of boundaries in CT volumetric data sets using dynamic , Comput. Biol. Med. (2015), http://dx.doi.org/10.1016/j.compbiomed.2015.10.018i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

8

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

L. Li et al. / Computers in Biology and Medicine ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Fig. 7. The user interface of our interactive TF design.

Fig. 8. Five iterative steps of the TF design of all the boundaries in the tooth data set using the dynamic M
j ∇f j histogram. The manipulation of each step is the same.

which boundary the intersection regions of arches belong. Users have to pick out part of an arch to avoid it. If the rendering result is unsatisfactory, users have to deal with the intersected regions via trial-and-error process. In LH histogram (Fig. 9(b)), the boundaries are represented as regions. Some of them are interconnecting; therefore, careful separation is needed. In M
j ∇f j histogram (Fig. 9 (c)) when Δh is set to 0, the boundaries are represented as vertical bars, and each bar is disjoint with another. On one hand, M
j ∇f j histogram is more concise to represent boundaries. On the other hand, there is a fixed strategy to extract all boundaries using M
j ∇f j histogram. Fig. 10 shows selections in M
j ∇f j histogram

and the corresponding coloring of arches in f
j ∇f j histogram and points and regions in LH histogram. 4.3. Medical data set example As our TF design method is a dynamic one, we record the entire process of TF design as videos. The videos can be watched at http://home.ustc.edu.cn/%7Elilu630/MF-Hist/. Fig. 11 shows the comparison of the visualizations of a CT head using f
j ∇f j histogram, LH histogram, M
j ∇f j histogram, 1-D histogram and region grow method, with the size of 256n256n256. The main difficulties of the boundary visualization lie in the spatial

Please cite this article as: L. Li, et al., Visualization of boundaries in CT volumetric data sets using dynamic , Comput. Biol. Med. (2015), http://dx.doi.org/10.1016/j.compbiomed.2015.10.018i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

L. Li et al. / Computers in Biology and Medicine ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

9

Fig. 9. The data set is the CT tooth. (a) The f
j ∇f j histogram. (b) The LH histogram. (c) The M
j ∇f j histogram when Δh ¼ 0. Boundaries are represented as arches, points, and regions, vertical bars in the three kinds of histograms respectively.

Fig. 10. (a) The rendering of the tooth using our method. (b) The corresponding coloring of arches in f
j ∇f j histogram. (c) The corresponding coloring of regions and points in LH histogram.

closeness and the amount of the boundaries. Fortunately, there exists little noise in the data set. The TF design manner of f
j ∇f j histogram and LH histogram is a trial-and-error process even with the help of widgets. Users have to fine-tune the widgets carefully to avoid misclassification without other technical supports. Thus, the trial-and-error process is tedious for non-experts that are lack of domain knowledge about the histograms. Fig. 11(d) shows the visualization using f
j ∇f j histogram. The shape of the boundaries is complete, but there are misclassification regions. The shin is rendered in blue and green. Fig. 11(f) shows the visualization using LH histogram. The shape is not complete, especially the blue skin. Fig. 11(b) shows the visualization using M
j ∇f j histogram. The boundaries are more complete and with less misclassification regions. In addition, users have to try all the regions in f
j ∇f j histogram and LH histogram to extract all boundaries. However, some boundaries may be ignored because of the limitation of the widgets and the noise. It is difficult to precisely determine the number of the boundaries by f
j ∇f j histogram and LH histogram. Compared with f
j ∇f j histogram and LH histogram, the boundary extraction process using M
j ∇f j histogram is an iterative way without concerning the number of the boundaries. 1-D histogram (Fig. 11(g)) and region grow method (Fig. 11(h)) mainly aim to extract the structures in the volume. In this instance, the main difference between our method and these two methods is that our method renders the boundary of the skins, while 1-D histogram and region grow method render the whole skins and muscles. 4.4. More complicated medical data sets Fig. 12 shows the rendering results of Chapel Hill CT Head, with the size of 256n256n113. It is more sophisticated than the example discussed above. The difficulties lie mainly in two aspects. Firstly, the head is divided into skin, bone and meningeal. Each of the three materials is close in spatial. Moreover, the ranges of j ∇f j of the skin and the meningeal are the same. Secondly, the scalar values of the two boards on the back of the head are

approximately equal to the skins and the two boards are so close that they may be classified as the same board. Fig. 12 gives the comparison between the visualizations of the Chapel Hill CT Head using f
j ∇f j histogram, LH histogram, M
j ∇f j histogram 1-D histogram and region growing method. The misclassification using f
j ∇f j histogram is shown clearly (Fig. 12(d)). The head is a mixed color of red, yellow and green because of the closeness of the skin, bone and meningeal. When compared with the rendering results generated by LH histogram (Fig. 12(f)) with M
j ∇f j histogram (Fig. 12(b)), the boundaries generated by M
j ∇f j histogram are more complete overall. The notable difference is the blue board. The blue board of (b) is pure blue, while the blue board of (f) mixes with yellow. In Fig. 12 (c) and (e), LH histogram cannot provide meaningful information, and in f
j ∇f j histogram other arches are hidden in the two obvious arches. It is difficult to determine the amount and the precise locations of the regions in the two histograms. Fig. 13 shows the rendering results of foot, with the size of 256n256n256. The scalar values of the bones and the muscles fluctuate greatly. The rendering results aiming to visualize the structures (Fig. 13(g) and (h)) show better results. In Fig. 13(c) and (e), f
j ∇f j histogram does not show any arches and LH histogram does not provide meaningful information. The rendering results are gotten by trail-and-error process without a clear purpose. Fig. 13(b) shows the rendering result using our method. In M
j ∇f j histogram, there are a clear distinction between the boundary of the bones and the boundary of the joints. The joints are rendered in a different color from the bones, which is the main difference from other methods.

5. Conclusions and future work In this paper, on one hand, we propose a generalized boundary model contaminated by noise and prove M has a good statistical property. On the other hand, we propose an interactive TF design method employing dynamic M
j ∇f j histogram. Users extract the

Please cite this article as: L. Li, et al., Visualization of boundaries in CT volumetric data sets using dynamic , Comput. Biol. Med. (2015), http://dx.doi.org/10.1016/j.compbiomed.2015.10.018i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

10

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

L. Li et al. / Computers in Biology and Medicine ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Fig. 11. (a) A slice of the CT head data. (b) The rendering result using M
j ∇f j histogram. (c) The f
j ∇f j histogram and the trapezoidal regions painted the same colors with the extracted boundaries. (d) The rendering result using f
j ∇f j histogram. (e) The LH histogram and the rectangular regions painted the same colors with the extracted boundaries. (f) The rendering result using LH histogram. (g) The rendering result using 1-D histogram. (h) The rendering result using region growing method.

Fig. 12. (a) A slice of Chapel Hill CT Head. (b) The rendering result using M
j ∇f j histogram. (c) The f
j ∇f j histogram and the regions painted the same colors with the extracted boundaries. (d) The rendering result using f
j ∇f j histogram. (e) The LH histogram and the regions painted the same colors with the extracted boundaries. (f) The rendering result using LH histogram. (g) The rendering result using 1-D histogram. (h) The rendering result using region growing method. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

vertical bars one by one with Δh from high to low. This is a convenience interactive strategy. However, there exist several limitations, which can be overcome to improve the quality of visualizations. First, we employ the C þ þ container of queue to implement the region growing and region elimination on CPU. This serial implementation leads to the non-real time interaction. In the future, the serial implementation can be replaced by a parallel implementation using CUDA. Second, although we add noise to the boundary, the visualization is not satisfactory when a strong bias occurs in the data set. In the future, two methods can be applied to improve the rendering quality.

Multi-scale local extremum search method may be a solution of the instability of the traditional line search. The noise in large scale can be eliminated easily, while the locations of the boundaries are blurred. It is a challenge to balance the pros and cons of the multiscale extremum search method. Local inherent properties employed in this work are scalar value and gradient. There are other properties such as curvature, texture, etc. How to select and combine these properties is worth a further investigation. Recently, deep learning method brings a new understanding of these properties. Farabet et al. [37] employed ConvNets to extract the properties and label scene for a 2-D image.

Please cite this article as: L. Li, et al., Visualization of boundaries in CT volumetric data sets using dynamic , Comput. Biol. Med. (2015), http://dx.doi.org/10.1016/j.compbiomed.2015.10.018i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

L. Li et al. / Computers in Biology and Medicine ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

11

Fig. 13. (a) A slice of foot. (b) The rendering result using M
j ∇f j histogram. (c) The f
j ∇f j histogram and the regions painted the same colors with the extracted boundaries. (d) The rendering result using f
j ∇f j histogram. (e) The LH histogram and the regions painted the same colors with the extracted boundaries. (f) The rendering result using LH histogram. (g) The rendering result using 1-D histogram. (h) The rendering result using region growing method.

This inspires us to employ deep learning method in the visualization of 3-D data set in our future work.

(b) We can use the same method of Z to calculate the probability distribution of W. (c) Let PðS r sÞ be the probability of S r s, where s is an arbitrary real number:

Conflict of interest statement PðS r sÞ ¼ PðX r 2sÞ ¼

None declared.

Acknowledgments This project was supported by National Natural Science Foundation of China (Grant nos. 61172037, 60871087).

Appendix A Let X  Nðμx ; σ x Þ and Y  Nðμy ; σ y Þ, X and Y are independent of each other. Denote ¼ X þffi Y, W ¼ X
Y, S ¼ X2 , R ¼ j X j : qZffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 þ σ 2 Þ. (a) Z  Nðμx þ μy ; q σffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x y ffi (b) W  Nðμx
μy ; σ 2x þ σ 2y Þ. (c) S  Nðμx ; 2σ x Þ. (d) R  2Nðμx ; σ x Þ; R Z 0. Proof. (a) Let pz(z) be the probability distribution of Z: Z þ1 pz ðzÞ ¼ px ðxÞpy ðz
xÞ dx

 1 exp
2 ðx
μx Þ2 dx 2σ x
1 2s

So we can get the probability distribution of S: ( ) dðPðS r sÞÞ 1 1  μx 2 pffiffiffiffiffiffiexp
¼ s
ps ðsÞ ¼ ds 2 2σ x 2π 2ð2σ x Þ2

ðA:4Þ

ðA:5Þ

(d) Let PðR r rÞ be the probability of R r r, where r is an arbitrary real number: PðR r rÞ ¼ Pð
r r X rrÞ ¼

1 pffiffiffiffiffiffi σ x 2π

Z

 1 exp
2 ðx
ux Þ2 dx 2σ x
r r

ðA:6Þ So we can get the probability distribution of R: pr ðrÞ ¼

 dðPðR rrÞÞ 2 1 ¼ pffiffiffiffiffiffiexp
2 ðx
ux Þ2 dr 2σ x σ x 2π

ðA:7Þ

Appendix B

where px and py are the probability distributions of X and Y respectively. ( "    #) z
x
μy 2 1 1 x
μx 2 exp
þ px ðxÞpy ðz
xÞ ¼ 2πσ x σ y 2 σx σy ðA:2Þ


ðz
μx
μy Þ2 1 pz ðzÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiexp 2ðσ 2x þ σ 2y Þ 2π ðσ 2x þ σ 2y Þ

Z

ðA:1Þ


1

Bring Eqs. (A.1) to (A.2):

1 pffiffiffiffiffiffi σ x 2π

!

Let X be a random variable and p(x) be the probability density function of X. In the range of ða; bÞð
1r a ob r þ1Þ, p(x) satisfies Pða o X obÞ ¼ 1, where P(x) is the probability function. Let Y ¼ f ðxÞ be a monotonic and continuous function in ða; bÞ. g(y) is the 0 inverse function of f(x) and g ðyÞ exists in ðα; βÞ which is the existence interval of g(y). Let ( pðgðyÞÞj g 0 ðyÞj yA ðα; β Þ ðB:1Þ qðyÞ ¼ 0 others

ðA:3Þ q(y) is the probability density function of Y.

Please cite this article as: L. Li, et al., Visualization of boundaries in CT volumetric data sets using dynamic , Comput. Biol. Med. (2015), http://dx.doi.org/10.1016/j.compbiomed.2015.10.018i

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

L. Li et al. / Computers in Biology and Medicine ∎ (∎∎∎∎) ∎∎∎–∎∎∎

12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 Q4 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

Proof. Without loss of generality, let f(x) be a monotone increasing function. For u A ða; bÞ: Z gðuÞ PðY ruÞ ¼ Pðf ðXÞ r uÞ ¼ PðX r gðuÞÞ ¼ pðxÞ dx ðB:2Þ a

Let x ¼ gðyÞ: Z Z u pðgðyÞÞj g 0 ðyÞj dy ¼ PðY ruÞ ¼ a

u

qðyÞ dy

ðB:3Þ

a

when u r α, PðY r uÞ ¼ 0, when u Z β, PðY r uÞ ¼ 1.

References [1] J. Klein, O. Friman, M. Hadwiger, B. Preim, F. Ritter, A. Vilanova, G. Zachmann, D. Bartz, Visual computing for medical diagnosis and treatment, Comput. Graph. 33 (4) (2009) 554–565. [2] R. Cheirsilp, R. Bascom, T.W. Allen, W.E. Higgins, Thoracic cavity definition for 3D PET/CT analysis and visualization, Comput. Biol. Med. 62 (2015) 222–238. [3] E.A. Rashed, M. al Shatouri, H. Kudo, Sparsity-constrained three-dimensional image reconstruction for c-arm angiography, Comput. Biol. Med. 62 (2015) 141–153. [4] M. Levoy, Display of surfaces from volume data, Comput. Graph. Appl. 8 (3) (1988) 29–37, IEEE. [5] R.A. Drebin, L. Carpenter, P. Hanrahan, Volume rendering, In: ACM Siggraph Computer Graphics, vol. 22, ACM, 1988, pp. 65–74. [6] G. Kindlmann, J.W. Durkin, Semi-automatic generation of transfer functions for direct volume rendering, In: IEEE Symposium on Volume Visualization, 1998, pp. 79–86. [7] P. Sereda, A.V. Bartroli, I.W.O. Serlie, F.A. Gerritsen, Visualization of boundaries in volumetric data sets using LH histograms, IEEE Trans. Vis. Comput. Graph. 12 (2) (2006) 208–218. [8] J. Kniss, G. Kindlmann, C. Hansen, Interactive volume rendering using multidimensional transfer functions and direct manipulation widgets, In: Proceedings of the Conference on Visualization'01, IEEE Computer Society, 2001, pp. 255–262. [9] J. Kniss, G. Kindlmann, C. Hansen, Multidimensional transfer functions for interactive volume rendering, IEEE Trans. Vis. Comput. Graph. 8 (3) (2002) 270–285. [10] E.L. Nickoloff, R. Riley, A simplified approach for modulation transfer function determinations in computed tomography, Med. Phys. 12 (4) (1985) 437–442. [11] E.B. Lum, K.-L. Ma, Lighting transfer functions using gradient aligned sampling, In: Proceedings of the Conference on Visualization'04, IEEE Computer Society, 2004, pp. 289–296. [12] I.W. Serlie, F.M. Vos, R. Truyen, F.H. Post, L.J. van Vliet, Classifying ct image data into material fractions by a scale and rotation invariant edge model, IEEE Trans. Image Process. 16 (12) (2007) 2891–2904. [13] J.-S. Praßni, T. Ropinski, K. Hinrichs, Efficient boundary detection and transfer function generation in direct volume rendering, In: 14th International Fall Workshop on Vision, Modeling, and Visualization (VMV 2009), 16–18 November 2009, Braunschweig, Germany, 2009, pp. 285–294. [14] P. Sereda, A. Vilanova, F.A. Gerritsen, Automating transfer function design for volume rendering using hierarchical clustering of material boundaries, In: EuroVis, 2006, pp. 243–250.

67 [15] B.P. Nguyen, W.-L. Tay, C.-K. Chui, S.-H. Ong, A clustering-based system to automate transfer function design for medical image visualization, Vis. Com68 put. 28 (2) (2012) 181–191. 69 [16] G. Kindlmann, R. Whitaker, T. Tasdizen, T. Moller, Curvature-based transfer 70 functions for direct volume rendering: methods and applications, In: Visualization, VIS 2003, IEEE, 2003, pp. 513–520. 71 [17] J. Hladuvka, A. König, E. Gröller, Curvature-based transfer functions for direct 72 volume rendering, In: Spring Conference on Computer Graphics, vol. 16, 73 Citeseer, 2000, pp. 58–65. [18] M. Haidacher, D. Patel, S. Bruckner, A. Kanitsar, M. Groller, Volume visualiza74 tion based on statistical transfer-function spaces, In: 2010 IEEE Pacific Visua75 lization Symposium (PacificVis), IEEE, 2010, pp. 17–24. 76 [19] E. LaMar, B. Hamann, K.I. Joy, Multiresolution techniques for interactive texture-based volume visualization, In: Electronic Imaging, International 77 Society for Optics and Photonics, 2000, pp. 365–374. 78 [20] A. Depeursinge, A. Foncubierta-Rodriguez, D. van de Ville, H. Müller, Three79 dimensional solid texture analysis in biomedical imaging: review and opportunities, Med. Image Anal. 18 (1) (2014) 176–196. 80 [21] S. Takahashi, Y. Takeshima, I. Fujishiro, Topological volume skeletonization and 81 its application to transfer function design, Graph. Models 66 (1) (2004) 24–49. 82 [22] D.X. Ji, S.H. Ong, K.W.C. Foong, A level-set based approach for anterior teeth 83 segmentation in cone beam computed tomography images, Comput. Biol. Med. 50 (2014) 116–128. 84 [23] A. Gooya, H. Liao, I. Sakuma, Generalization of geometrical flux maximizing 85 flow on Riemannian manifolds for improved volumetric blood vessel seg86 mentation, Comput. Med. Imaging Graph. 36 (6) (2012) 474–483. [24] H. Qin, B. Ye, R. He, The voxel visibility model: an efficient framework for 87 transfer function design, Comput. Med. Imaging Graph. Q2 88 [25] Y. Wang, J. Zhang, W. Chen, H. Zhang, X. Chi, Efficient opacity specification 89 based on feature visibilities in direct volume rendering, in: Computer Graphics Forum, vol. 30, Wiley Online Library, 2011,pp. 2117–2126. 90 [26] L. Zhou, C. Hansen, Transfer function design based on user selected samples 91 for intuitive multivariate volume exploration, In: 2013 IEEE Pacific Visualiza92 tion Symposium (PacificVis), IEEE, 2013, pp. 73–80. [27] C.R. Salama, M. Keller, P. Kohlmann, High-level user interfaces for transfer 93 function design with semantics, IEEE Trans. Vis. Comput. Graph. 12 (5) (2006) 94 1021–1028. 95 [28] P. Rautek, S. Bruckner, M.E. Gröller, Semantic layers for illustrative volume rendering, IEEE Trans. Vis. Comput. Graph. 13 (6) (2007) 1336–1343. 96 [29] H. Guo, X. Yuan, Local WYSIWYG volume visualization, In: 2013 IEEE Pacific 97 Visualization Symposium (PacificVis), IEEE, 2013, pp. 65–72. 98 [30] E. Shen, S. Li, X. Cai, L. Zeng, W. Wang, Sketch-based interactive visualization: 99 a survey, J. Vis. 17 (4) (2014) 275–294. [31] I. Serlie, F. Vos, H. Venema, L. van Vliet, CT imaging characteristics, Technical 100 Report QI-2006-01, Delft University of Technology. 101 [32] J. Hsieh, Computed Tomography: Principles, Design, Artifacts, and Recent 102 Advances, SPIE, Bellingham, WA, 2009. [33] S.O. Rice, Mathematical analysis of random noise, Bell Syst. Tech. J. 23 (3) 103 (1944) 282–332. 104 [34] T. Li, X. Li, J. Wang, J. Wen, H. Lu, J. Hsieh, Z. Liang, Nonlinear sinogram 105 smoothing for low-dose X-ray CT, IEEE Trans. Nucl. Sci. 51 (5) (2004) 2505–2513. 106 [35] J. Wang, T. Li, H. Lu, Z. Liang, Penalized weighted least-squares approach to 107 sinogram noise reduction and image reconstruction for low-dose X-ray 108 computed tomography, IEEE Trans. Med. Imaging 25 (10) (2006) 1272–1283. [36] J. Nocedal, S. Wright, Numerical Optimization, Springer Science & Business 109 Media, 2006. 110 [37] C. Farabet, C. Couprie, L. Najman, Y. LeCun, Learning hierarchical features for 111 scene labeling, IEEE Trans. Pattern Anal. Mach. Intell. 35 (8) (2013) 1915–1929.

Please cite this article as: L. Li, et al., Visualization of boundaries in CT volumetric data sets using dynamic , Comput. Biol. Med. (2015), http://dx.doi.org/10.1016/j.compbiomed.2015.10.018i

112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132