Yang Chen

Yinsheng Li

Weimin Yu

Limin Luo

Wufan Chen

Christine Toumoulin

Joint-MAP Tomographic Reconstruction with Patch Similarity Based Mixture Prior Model

Keywords: Maximum a posteriori (MAP), Markov Random Fields (MRF), Patch Similarity Mixture (PSM) prior model, Tomographic Reconstruction, Compute Unified Device Architecture (CUDA). I

Tomographic reconstruction from noisy projections do not yield adequate results. Mathematically, this tomographic reconstruction represents an ill-imposed problem due to information missing caused from the presence of noise. Maximum A Posteriori (MAP) or Bayesian reconstruction methods offer possibilities to improve the image quality as compared with analytical methods in particular by introducing a prior to guide the reconstruction and regularize the noise. With an aim to achieve robust utilization of continuity/connectivity information and overcome the heuristic weight update for other nonlocal prior methods, this paper proposed a novel Patch Similarity based Mixture (PSM) prior model for tomographic reconstruction. This prior is defined by a weighted Gaussian distance between neighborhood intensities. The weight quantifies the similarity between local neighborhoods and is computed using a maximization entropy constraint. This prior is then introduced within a joint image/weight MAP CT reconstruction algorithm. Several acceleration trials including Compute Unified Device Architecture (CUDA) parallelization is applied to alleviate the intensive patch distance computation involved in the joint algorithm. The method was tested with both synthetic phantoms and clinical CT data and compared in accuracy with five other reconstruction algorithms which are FBP and Bayesian-based. Reconstruction results show that the proposed reconstructions are able to produce high-quality images with ensured iteration convergence.

Nevertheless, this reconstruction becomes an ill-posed problem when it is carried out from a limited number of projections or has to deal with inconsistencies generated by the presence of noise in the measurements. Several applications are concerned with this problem such as low-dose computed tomography (CT), Positron Emission Tomography (PET) and Single Positron Emission Computerized Tomography (SPECT). The question of how to overcome such ill-posedness has been widely studied in the past twenty years [START_REF] Bertero | Ill posed problems in early vision[END_REF][START_REF] Lange | Convergence of EM image reconstruction algorithms with Gibbs smoothness[END_REF][START_REF] Erdogan | Monotonic algorithms for transmission tomography[END_REF][START_REF] Chlewicki | Noise reduction and convergence of Bayesian algorithms with blobs based on the Huber function and median root prior[END_REF][START_REF] Bouman | A generalized Gaussian image model for edge-preserving MAP estimation[END_REF][START_REF] Duval | The TVL1 model: A geometric point of view[END_REF][START_REF] Chen | Bayesian statistical reconstruction for low-dose X-ray computed tomography using an adaptive-weighting nonlocal prior[END_REF][START_REF] Yu | Edge-Preserving tomographic reconstruction with nonlocal regularization[END_REF]. Among the proposed solutions, a large number of work suggested solving it by imposing a priori knowledge as a regularization constraint in particular within a Bayesian reconstruction framework. The regularization prior, associated with the likelihood function, provide thus an approximate solution of the reconstruction problem as the maximum a priori (MAP) estimate that maximizes the log of the posterior probability [START_REF] Li | Markov Random Field Modeling in image Analysis[END_REF].

A large variety of prior has been proposed in the literature. Quadratic priors such as Gaussian-Markov prior or the convex quadratic membrane (QM) prior [START_REF] Lange | Convergence of EM image reconstruction algorithms with Gibbs smoothness[END_REF][START_REF] Li | Markov Random Field Modeling in image Analysis[END_REF], have the particularity to bring in a global smoothness of the solution with a tendency in blurring edges in the image. To preserve edges in regularization, some nonquadratic priors with nonquadratic energy potentials have been widely applied in tomographic reconstructions [START_REF] Chlewicki | Noise reduction and convergence of Bayesian algorithms with blobs based on the Huber function and median root prior[END_REF][START_REF] Bouman | A generalized Gaussian image model for edge-preserving MAP estimation[END_REF][START_REF] Duval | The TVL1 model: A geometric point of view[END_REF][START_REF] Chen | Bayesian statistical reconstruction for low-dose X-ray computed tomography using an adaptive-weighting nonlocal prior[END_REF][START_REF] Yu | Edge-Preserving tomographic reconstruction with nonlocal regularization[END_REF]. Some classic edge-preserving nonquadratic priors are the Huber potential function [START_REF] Lange | Convergence of EM image reconstruction algorithms with Gibbs smoothness[END_REF], Median Root Prior (MRP) [START_REF] Chlewicki | Noise reduction and convergence of Bayesian algorithms with blobs based on the Huber function and median root prior[END_REF] and TV prior [START_REF] Duval | The TVL1 model: A geometric point of view[END_REF]. The main drawbacks of these edge-preserving priors are the unguaranteed convexity, and the easily produced blocky artifacts when the measurements suffer from heavy noise [START_REF] Lange | Convergence of EM image reconstruction algorithms with Gibbs smoothness[END_REF][START_REF] Bouman | A generalized Gaussian image model for edge-preserving MAP estimation[END_REF][START_REF] Duval | The TVL1 model: A geometric point of view[END_REF][START_REF] Chen | Bayesian statistical reconstruction for low-dose X-ray computed tomography using an adaptive-weighting nonlocal prior[END_REF][START_REF] Yu | Edge-Preserving tomographic reconstruction with nonlocal regularization[END_REF].

In [START_REF] Bertero | Ill posed problems in early vision[END_REF][START_REF] Lange | Convergence of EM image reconstruction algorithms with Gibbs smoothness[END_REF][START_REF] Erdogan | Monotonic algorithms for transmission tomography[END_REF][START_REF] Chlewicki | Noise reduction and convergence of Bayesian algorithms with blobs based on the Huber function and median root prior[END_REF][START_REF] Bouman | A generalized Gaussian image model for edge-preserving MAP estimation[END_REF][START_REF] Duval | The TVL1 model: A geometric point of view[END_REF][START_REF] Chen | Bayesian statistical reconstruction for low-dose X-ray computed tomography using an adaptive-weighting nonlocal prior[END_REF][START_REF] Yu | Edge-Preserving tomographic reconstruction with nonlocal regularization[END_REF][START_REF] Li | Markov Random Field Modeling in image Analysis[END_REF], prior computation usually relies on intensity differences between individual pixels and a simple spatial weighting strategy between neighboring pixels. Thus to take in more information on the image (geometry, statistics distribution of the structures of interest, ), the idea has been to compute this difference between neighboring patches and apply weights that reflect the similarity between pixels within the neighborhood. This idea was inspired from Buades et al.'s work [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] on the nonlocal mean filter designed for image denoising. Our previous work [START_REF] Chen | Bayesian statistical reconstruction for low-dose X-ray computed tomography using an adaptive-weighting nonlocal prior[END_REF] and [START_REF] Chen | Nonlocal prior Bayesian tomographic reconstruction[END_REF] proposed thereby in the prior function, an adaptive weight as a decreasing function of the Euclidean distance between two square neighboring patches. This weight expressed a similarity between the neighboring pixels and was updated using the latest estimated image in a one step late (OSL) way. This method is however highly costly in computing all the patch distances and may provide inconsistent posterior energy function with some problem of convergence in some situations.

This paper aims at solving these problems by proposing a novel Patch Similarity based Mixture (PSM) prior model within a new convergent image/weight joint reconstruction algorithm. This PSM prior model provides as in [START_REF] Chen | Nonlocal prior Bayesian tomographic reconstruction[END_REF], a weighted sum of the distances between patches surrounding each pair of neighboring pixels but weights are now estimated using a constrained entropy maximization, in which a normalization is introduced to preserve the balance between the regularization and image estimate. This proposed PSM prior and the statistical framework for MAP reconstruction are described in section II. To overcome the induced computational burden, a fast implementation of the most time consuming part of the algorithm is performed in Graphics programming Units (GPU) using Compute Unified Device Architecture (CUDA) framework [START_REF]NVIDIA CUDA T M Programming Guide[END_REF][START_REF]Accelerating MATLAB with CUDA Using MEX Files[END_REF][START_REF]GPU Acceleration in MATLAB[END_REF]. An comparative evaluation of this algorithm was afterwards conducted both visually and quantitatively on simulated and real CT data. This evaluation included a study on the iteration monotonicity and sensitivity of the result according to the sizes of the comparing patches and the searching neighborhood windows. The set of results are given in section III.

II. THEORY

A. The Patch Similarity Based Mixture (PSM) Prior Model

A MAP estimate of the image can be computed as the maximization of the posterior probability P (f |g), which is defined by:

P (f |g) ∝ P (g|f ) P (f ) (1)
P (f ) = Z -1 × exp (-βU (f )) (2)
U (f ) = j U j (f ) = j b∈Nj w bj v (f b -f j ) (3)
where, f j is the estimate for a pixel j in the image space, the potential function v(f ) relates to an intensity difference between the pixel j and the neighboring pixel b. N j denotes the neighborhood window centered at pixel j. In traditional prior such as the Huber one, w bj defines the weighting coefficient that is put to 1 for a neighboring pixel located in the horizontal and vertical direction and 1/ √ 2 when it refers to a neighboring pixel whose location is in the diagonal direction. They also usually express an intensity difference between individual pixels within local neighborhoods [START_REF] Lange | Convergence of EM image reconstruction algorithms with Gibbs smoothness[END_REF][START_REF] Erdogan | Monotonic algorithms for transmission tomography[END_REF][START_REF] Chlewicki | Noise reduction and convergence of Bayesian algorithms with blobs based on the Huber function and median root prior[END_REF][START_REF] Bouman | A generalized Gaussian image model for edge-preserving MAP estimation[END_REF][START_REF] Duval | The TVL1 model: A geometric point of view[END_REF][START_REF] Chen | Bayesian statistical reconstruction for low-dose X-ray computed tomography using an adaptive-weighting nonlocal prior[END_REF][START_REF] Yu | Edge-Preserving tomographic reconstruction with nonlocal regularization[END_REF][START_REF] Li | Markov Random Field Modeling in image Analysis[END_REF]. This expression is often inadequate to characterize large density distribution variations on the image. But this issue is made difficult by the fact we have to simultaneously satisfy a double constraint that is the image has to be locally smoothed while allowing sharp intensity transitions that correspond to tissue, organ or lesion boundaries. A way to take into account this double constraint is to introduce adaptive weights in the smoothing prior so that the potential function v(f ) may be set to 0 in case of dissimilarity between neighboring pixels to avoid smoothing between these pixels and preserve thus the transitions at the boundary of two regions.

We proposed in [START_REF] Chen | Nonlocal prior Bayesian tomographic reconstruction[END_REF] a spatially adaptive nonlocal (NL) prior model in which v(t) = t 2 /2 was used and w in

(3) was determined by a similarity criterion between patches surrounding each pair of neighboring pixels within a neighborhood window of size N (Fig. 1). In [START_REF] Lyu | An implicit Markov random field model for the multi-scale oriented representations of natural images[END_REF], weights w and image f were alternatively updated within a one step late (OSL) algorithm. However, this joint updated reconstruction algorithm although efficient in practice in providing a spatially adaptive regularization, could lead to inconsistent posterior energy function with no guaranty of convergence within each iteration due to the non convexity of the potential function [START_REF] Lange | Convergence of EM image reconstruction algorithms with Gibbs smoothness[END_REF].

To overcome this problem, we propose here a novel PSM prior model, in which regularization takes effect though penalizing the Gaussian weighted Euclidean distance ||G α * n b -G α * n j || between the patches n b and n j surrounding each pair of pixel j and its neighboring pixel b in neighborhood N j . Here, G α is a two-dimensional Gaussian kernel of preset standard deviation α. All the weights w in each N are considered unknown variables with nonlinear dependence on the image f . This PSM prior can be considered as a high-level MRF prior that models the high-order patterns around each pixel in neighborhood N [START_REF] Li | Markov Random Field Modeling in image Analysis[END_REF], [START_REF] Lyu | An implicit Markov random field model for the multi-scale oriented representations of natural images[END_REF][START_REF] Winkler | Image Analysis, Random Fields And Markov Chain Monte Carlo Methods[END_REF]. For each pair of pixel j and b in neighborhood N j , w bj is considered as the similarity probability between the two patches n b and n j surrounding pixels b and j [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF]. We applied then the maximum entropy principle to derive a probability distribution which provides the best consistency in modeling the statistical properties of the image. As we looked for obtaining global patterns via stochastic propagation of local interactions, we choose thus the probability distribution that has the maximum entropy [START_REF] Peyr | Non-local regularization of inverse problems[END_REF][START_REF] Zhu | FRAME: filters, random field And Maximum Entropy:Towards a Unified Theory for Texture Modeling[END_REF]. This was done by introducing in the prior, an entropy term j b∈Nj w bj lnw bj on all the weights w bj attached to each neighborhood window N j with the following normalization constraint b∈Nj w bj = 1.

A large value of w bj is associated with an elevated similarity probability between patches n b and n j and a low entropy prior energy and conversely, a small value of w bj leads to a low similarity probability between patches n b and n j and a higher entropy prior energy.

B. Joint Weight/image Reconstruction Algorithm with the application in CT

The joint-MAP reconstruction algorithm can be therefore rewritten, according to the Bayesian theorem in order to build a posterior probability P (f, w, η|g) that includes the parameters w and Lagrange-multiplier parameter η [START_REF] Li | Markov Random Field Modeling in image Analysis[END_REF][START_REF] Hsiao | Bayesian image reconstruction for transmission tomography using deterministic annealing[END_REF] such as:

P (f, w, η|g) ∝ P (g|f )P (f |w, η)P (w, η) (4)
Substituting the log prior of the Gibbs distribution in (2) produces the objective posterior energy function:

ψ β (f, w, η|g) = log P (g|f ) + log P (f |w, η) + log P (w , η) = L(g, f ) -β( j U P SM (f j , w, η| f )) = L(g, f ) -β j (( b∈Nj w bj ||G α * n b -G α * n j || + λ b∈Nj w bj ln w bj ) + η j ( b∈Nj w bj -1)) (5)
where, f denotes the current estimated image. For each pixel j, this PSM prior uses the mixture prior energy U P SM (f j , w, η| f ) as a constrained combination of a pair-wise weighted patch similarity prior energy and an entropy prior energy -

j b∈Nj
w bj lnw bj . The parameter λ controls the entropy constraint on w , and is routinely set based on the noise level in measurement data.

The proposed PSM prior is applied in statistical CT image reconstruction. We model the scanned photon counts vectors as Poisson random variables. The likelihood function P L (g/f ) for a given CT image f , is described by the independent Poisson distribution [START_REF] Elbakri | Efficient and accurate likelihood for iterative image reconstruction in X-ray computed tomography[END_REF]:

P L (g|f ) = I i=1 [(ḡ i (f ) gi /g i !) exp (-ḡ i (f ))] (6)
with the mean

ḡi (f ) = d i exp (-[Af ] i ) = d i exp ⎛ ⎝ - J j a ij f j ⎞ ⎠ (7)
where, g i is the measured photon counts detected by projection i can be modeled as Poisson random variables with ḡi a function of the underlying attenuation map [START_REF] Kat | Principles of Computerized Tomographic imaging[END_REF][START_REF] Zhang | Fast and robust deconvolution-based image reconstruction for photoacoustic tomography in circular geometry: experimental validation[END_REF][START_REF] Beekman | Ordered subset reconstruction for x-ray CT[END_REF][START_REF] Elbakri | Efficient and accurate likelihood for iterative image reconstruction in X-ray computed tomography[END_REF]. I is the number of projection lines, J is the number of image pixels and a ij denotes the probability that a transmission from pixel j is detected by the detector i. d i factors represents the photon number detected by detector i in the absence of the absorber. Thus the log-likelihood function for obtaining g from f is:

L(g, f ) = log P L (g|f ) = i (g i log ḡi (f ) -ḡi (f ) -log g i !) (8)
Then, the objective CT image f and the weights w can be jointly learned through maximizing the ψ β (f, w, η|g) in ( 5) with respect to each element of f and w alternatively [START_REF] Hsiao | Bayesian image reconstruction for transmission tomography using deterministic annealing[END_REF]:

1) Weight update (with fixed image f ): arg max w (-U P SM (w, η| f ))

Let f denotes the current image estimate, each weight w bj is updated by solving arg max w bj (-U P SM (w, η| f ), b ∈ N j ):

w bj = arg max w bj (- b∈Nj w bj ||G α * n b -G α * n j || -λ b∈Nj w bj ln w bj -η j ( b∈Nj w bj -1)) (9)
w bj is computed in the following way:

∂UP SM ( f,w,η) ∂w bj = 0 ⇒ G α * n b -G α * n j + λ(1 + ln w bj ) + η j = 0 ∂UP SM ( f,w,η) ∂ηj = 0 ⇒ b∈Nj w bj -1 = 0 ⎫ ⎪ ⎬ ⎪ ⎭ ⇒ ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ w bj = exp(-Gα * n b -Gα * nj /λ) b ∈N j exp(-Gα * n b -Gα * nj /λ) η j = -G α * n b -G α * n j -λ -λ ln exp(-Gα * n b -Gα * nj /λ) b ∈N j exp(-Gα * n b -Gα * nj /λ) (10)
Here, G α * n b -G α * n j can be easily calculated using current image estimate f . From [START_REF] Chlewicki | Noise reduction and convergence of Bayesian algorithms with blobs based on the Huber function and median root prior[END_REF], we can see that each w bj is determined not only from the similarity metrics between patches n b and n j , but also with the summed similarity metrics between all pairs of patches. So w bj is different from w jb .

2) Image update (with fixed weight ŵ): arg max

f ψ β (f, ŵ, η|g)
With the estimated weight denoted by ŵ and the terms with no f omitted, the ψ β (f, ŵ|g) becomes

ψ β (f, ŵ|g) = L(g, f ) -β j ( b∈Nj ŵbj ||G α * n b -G α * n j ||)) (11)
Based on the theory of pseudo-likelihood and Iterated Conditional Mode (ICM) in [START_REF] Li | Markov Random Field Modeling in image Analysis[END_REF], we can get each f j by solving the following factorized problem:

arg max fj ψ β (f j , ŵbj |g, f ) ⇒ arg max fj L (f, g) -β( b∈Nj ŵbj ||G α * n b -G α * n j ||)) ( 12
)
The discretization property of the image pixels makes it suitable to apply the Iterative Coordinate Descent (ICD) algorithm to solve (12) [START_REF] Bouman | A unified approach to statistical tomography using coordinate descent optimization[END_REF][START_REF] Huang | Iterative scaling and coordinate descent methods for maximum entropy models[END_REF]. So in this study we use the Paraboloidal surrogate coordinate ascent (PSCA)

reconstruction in [START_REF] Erdogan | Monotonic algorithms for transmission tomography[END_REF] to solve image estimate at each pixel f j . As an ICD optimization procedure, the PSCA algorithm applies quadratic paraboloidal surrogate to achieve efficient reconstruction in transmission tomography [START_REF] Erdogan | Monotonic algorithms for transmission tomography[END_REF]. This PSCA method has a close theoretical relation with the Optimization Transfer idea in [START_REF] Lange | Optimization transfer using surrogate objective functions[END_REF], and the iterative scaling based CD optimization with a factorized resulting Newton step [START_REF] Huang | Iterative scaling and coordinate descent methods for maximum entropy models[END_REF]. For the PSCA reconstruction, the quadratic paraboloidal surrogate of the log-likelihood term and the quadratic prior term allow exact Taylor expansions, which can lead to an efficient exact Newton update.

When the measured projection data g is also the observed image data f , this reconstruction problem degenerates into a image denoising problem and the above joint weight/image update algorithm leads to the iterative denoising algorithm applied in the Gaussian noise case [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with Probabilistic Patch-Based Weights[END_REF]. It is trivial to prove that the -U P SM (w, η| f ) in ( 16) is concave with respect to w. ∂ 2 L (g, f ) ∂f b ∂f j is known to negative definite [START_REF] Lange | Globally convergent algorithms for maximum a posteriori transmission tomography[END_REF][START_REF] Lange | Convergence of EM image reconstruction algorithms with Gibbs smoothness[END_REF]. Thus considering [START_REF] Bouman | A generalized Gaussian image model for edge-preserving MAP estimation[END_REF], and given the fixed computed weight ŵ, we conclude that the Hessian matrix of each ψ β (f, ŵ|g) in ( 11) is also negative definite. So each ψ β (f, w, η|g) is separately concave with respect to each f and w. Denoting f k , ŵk , ηk and f k+1 , ŵk+1 , ηk+1 , as the k th and k + 1 th iterated estimations, the joint weight/image reconstruction algorithm yields to 10)), respectively. We can therefore infer the following relation:

the convexity of ||G α * n b -G α * n j || in
ψ β ( f k , ŵk+1 , ηk+1 |g) > ψ β ( f k , ŵk , ηk |g) and ψ β ( f k+1 , ŵk+1 , ηk+1 |g) > ψ β ( f k , ŵk+1 , ηk+1 |g) (η is updated with ŵ in (
ψ β ( f k+1 , ŵk+1 , ηk+1 |g) > ψ β ( f k , ŵk+1 , ηk+1 |g) > ψ β ( f k , ŵk , ηk |g), which implies that convergence to local maximum is ensured for the iterated objective function ψ β (f, w, η|g).
Finally, the joint reconstruction algorithm is outlined as follows:

-

--------------Algorithm -------Outline ----------------- • Define: L(g, f ) = I i=1 h i (l i ), h i (l i ) = (g i log (d i exp (-l i )) -d i exp (-l i ) ) , l i = [Af ] i = J j a ij f j
Initialize: set the initial image estimate f , built projection matrix A based on system geometry, and set iteration number niter.

• For each iteration k = 1, ......, niter:

With current estimated image f , update all η and ŵbj using [START_REF] Chlewicki | Noise reduction and convergence of Bayesian algorithms with blobs based on the Huber function and median root prior[END_REF], and calculate li = A

f i = J j a ij fj .
For each pixel j in image f , based on the PSCA method in [START_REF] Erdogan | Monotonic algorithms for transmission tomography[END_REF], we update the image estimate fj :

ḣi li = ⎛ ⎝ 1 - g i d i exp -li ⎞ ⎠ d i exp -li , ḧi li = d i g i exp -li d i exp -li 2 -d i exp -li c i li = ⎧ ⎨ ⎩ 2 l2 i h i li -h i (0) -li ḣi li , li > 0 -ḧi li li = 0, UPSM (f j , ŵ, η| f ) = ∂ ∂f j U P SM (f j , ŵ, η| f ) fj = fj = -2 b∈Nj ŵbj ( fb -fj ) ||G α * n b -G α * n j || ÜPSM (f j , ŵ, η| f ) = ∂ 2 ∂f 2 j U P SM (f j , ŵ, η| f ) fj = fj = 2 b∈Nj ŵbj ||G α * n b -G α * n j ||
We can update the image estimate f j :

fj = ⎡ ⎢ ⎢ ⎣ fj + I i=1 a ij ḣi ( li ) -β UPSM (f j , ŵ, η| f ) I i=1 a 2 ij c i ( li ) + β ÜPSM (f j , ŵ, η| f ) ⎤ ⎥ ⎥ ⎦ + End End

---------------Algorithm -------Outline -----------------

C. Graphic Hardware Acceleration of the PSM Prior Computation

As illustrated in above algorithm outline, the PSM prior reconstruction involves pixelwise calculations of the UPSM (f j , ŵ, η| f ) and ÜPSM (f j , ŵ, η| f ), in which intensive weight computations are required for each pair of patches with centers in each N j over the whole image. This leads to a computational burden that is not compatible with the medical practice.

CUDA offers an unified hardware and software solution for parallel computing on NIVIA GPUs supporting the standard C/C++ programming language together with high performance computing numerical libraries [START_REF]NVIDIA CUDA T M Programming Guide[END_REF].

A CUDA computation model is broken down into connected CPU (host) and GPU (device) parts. The direct heterogeneous CPU-GPU (Graphic Processing Unit) parallelism technique can be used to develop and deploy kernel functions for computations. The host executes code sequentially for GPU initialization, communication, and making kernel calls, and the GPU executes the kernel functions. Thus we can make well use of the parallel power to achieve time optimized result. Once the data loaded into the GPU memory, multi-kernel functions allow processing the data simultaneously [START_REF]NVIDIA CUDA T M Programming Guide[END_REF][START_REF]Accelerating MATLAB with CUDA Using MEX Files[END_REF][START_REF]GPU Acceleration in MATLAB[END_REF]. The CPU calls kernel functions and each call results in the creation of threads and blocks in the GPU.

……

For all pixel j with its neighborhood N j Parallelly calculate summed weighted intensity differences and the summed weights :

P1= exp j b j b j b j b N G n G n f f G n G n
For all pixel j with its neighborhood We exploited therefore this environment to improve the performance of the proposed reconstruction, proposing a parallel implementing of the PSM prior. The structure of the CUDA implementation is illustrated in Fig. 2. The total number of blocks in the grid is given by image row size and the total number of threads in each block by image column size. Therefore, at each iteration in the PSM prior reconstruction, all the threads in this block-grid structure are executed in parallel to perform all the pixelwise operations i.e those involved in the calculations of the weights for all patch pairs and then the pixelwise functions UPSM (f j , ŵ, η| f ) and ÜPSM (f j , ŵ, η| f ) .

N j P2= exp P3 exp j j b j b j b N b j b N G n G n G n G n G n G n

Parallel calculation for only 1 time

j Parallelly calculate ˆ, , | =P1/P3 and , , | =P2/P3 PSM j PSM j U f w f U f w f
In addition, we devised a more efficient method to calculate the distance dis bj between the two comparing patches n b and n j in [START_REF] Bouman | A generalized Gaussian image model for edge-preserving MAP estimation[END_REF] in each iteration. In the following (13), we expand the weighted patch distance in [START_REF] Bouman | A generalized Gaussian image model for edge-preserving MAP estimation[END_REF] into two summations of the G α weighted squared intensities in n b and n j , and one summations of the G α weighted dot product between n b and n j :

||G α * n b -G α * n j || = l (α l f l∈n b -α l f l∈nj ) 2 = l α l 2 (f 2 l∈n b + f 2 l∈nj -2f l∈n b × f l∈nj ) = l α l 2 f 2 l∈n b + l α l 2 f 2 l∈nj -2 l α l 2 (f l∈n b × f l∈nj ) (13)
Here, under the square root operator, the two first sums, that relates to each patch n b and n j , can be preliminary computed before performing the pixel-wise calculation at each iteration. The third sum related to the intensity dot product between n b and n j can also be considered as one element of the convolution matrix between the central patch n j and the other patches in the neighborhood N . Such convolution can be easily obtained with efficient multiplications, FFT (Fast Fourier Transform) and IFFT (Inverse FFT) operators. In addition, considering that the

distance ||G α * n b -G α * n j || is equal to the distance ||G α * n j -G α * n b ||,
we can save time by computing one distance on two.

III. EXPERIMENTATION AND ANALYSIS

In order to assess the algorithm effectiveness and speed, we present results on simulated and real data successively.

We compare these results with two kinds of reconstruction methods, which are the conventional FBP and Bayesianbased algorithms respectively. For each of them, we considered different filters and priors. The list below gives the characteristics of each method:

1) FBP1 rec.: Filtered Backprojection reconstruction considering a ramp filter and a cutoff frequency equal to the Nyquist frequency;

2) FBP2 rec.: Filtered Backprojection reconstruction considering a hamming window and a cutoff frequency equal to 80% of the Nyquist frequency;

3) TV prior rec. : MAP reconstruction Algorithm based on a Paraboloidal Surrogates optimization procedure [START_REF] Erdogan | Monotonic algorithms for transmission tomography[END_REF] and using a TV prior:

U (f ) = j Δ h j f 2 + Δ v j f 2 + 2 ;
4) Huber prior rec. : MAP reconstruction Algorithm based on a Paraboloidal Surrogates optimization procedure [START_REF] Erdogan | Monotonic algorithms for transmission tomography[END_REF] and using a Huber prior :

U (f ) = j U j (f ) = j b∈Nj w bj v (f b -f j ) with v(t) = ⎧ ⎨ ⎩ t 2 2, |t| ≤ γ γ |t| -γ 2 2, |t| > γ
(γ is the bulit-in parameter controlling the sensitivity to noise); In all the reconstructions, the negative intensities in reconstructed images are set to 0 to preserve intensity positivity. The stopping criterion, for reconstructions using TV, Huber and MRP prior, is given by the condition (k denotes the iteration number):

ψ β (f k+1 |g) -ψ β (f k |g) < 0.99(ψ β (f k |g) -ψ β (f k-1 |g)) ( 14
)
In the case of the PSM prior based reconstruction, this criterion relies on the current estimated image intensities and weights (k denotes the iteration number):

ψ β (f k+1 , w k+1 |g) -ψ β (f k , w k |g) < 0.99(ψ β (f k , w k |g) -ψ β (f k-1 , w k-1 |g)) (15)
In practice, between 800 and 1000 iterations are required for all theses methods to reach the stopping conditions. For the proposed PSM prior reconstruction, we use a 1-value unitary image as the initial condition because the streak artifacts in FBP initial condition tends to mess up the whole iteration process. We will illustrate this in the later part of the paper. It is noted in [START_REF] Zbijewski | Suppression of intensity transition artifacts in statistical X-ray computer tomography through Radon inversion initialization[END_REF] that, though showing a relative low convergence in the early iterations, reconstructions with a uniform initial condition are ultimately able to converge to results at least no worse than the ones with a FBP initial condition. For equal comparison, other Bayesian reconstructions also use 1-value unitary images as the initial conditions. Routinely, larger λ is set to suppress more noise when the noise level in the measured sinogram is higher, so we set λ as δ * var s with sinogram variance var s . It is hard to derive a deterministic relation between sinogram variances and the λ that will lead to the best reconstructions in terms of SNR or visual criterions. So we choose suitable λ by modulating the δ in δ * var s . In practical experiment, we found the suitable range for the parameter is from 8 to 14. This parameter determination strategy was also used in the Nonlocal means denoising in [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF].

The algorithms have been written in hybird MATLAB and C++ language using the MATLAB release R2009b

and NVIDIA CUDA libraries. The Operating System is Linux, the CPU core is INTEL E6600 and the GPU card is a NVIDIA GTX465 (Core Clock: 625MHz, Shader Clock: 1250MHz, Stream Processors: PCI Express 2.0 × 16

Memory: Effective Memory Clock: 3240MHz Memory Size: 1GB GDDR5 Memory). We build the projection matrix based on ASPIRE software system of Professor fessler's group was used as the programming platform [START_REF] Fessler | Aspire 3.0 user's guide: A sparse reconstruction library[END_REF].

A. Reconstruction Results on Simulated Data

Two synthetic phantoms have been considered that were both gray scale images whose intensity ranged within the [0-255] interval with a resolution of 8 bit 256 × 256 pixels and a pixel size of 1mm × 1mm. The first phantom was made of an ellipsoid of uniform background in which a cold square and circle was introduced, the circle having blurred edges and including a smaller hot square for simulating a lesion. The second phantom was a modified shepp-logan phantom in which 4 straight lines of different width was added. We modeled a monochromatic single slice detector CT scanner with a parallel beam geometry. We simulated a CT scan acquisitions with 367 radial bins and 360 angles uniformly spaced to produce a sinogram comprising 360 angular projections, each containing 367 linear samples. For the first phantom, we simulated a low-dose sinogram with a total count 0.7 × 10 8 (further named Phantom1-SL), and considered sinograms with two levels of total counts for the second phantom :1.2 × 10 8 and 4 × 10 8 (further named Phantom2-SL and Phantom1-SH). All the sinograms were generated by adding a 10% zero mean Gaussian noise (with a variance of 0.5) to take into account the electronic and quantum noises. The d i factor in [START_REF] Bertero | Ill posed problems in early vision[END_REF] were generated using pseudo-random log normal variates with a standard deviation of 0.3 to take into account the characteristics of the energy-integrating sensor in the X-Ray CT detector. The two original phantoms are displayed in Fig. 3.

Parameter setting related to the Bayesian based reconstructions are given Table I for each simulated phantoms Phantom1-SL, Phantom2-SL and Phantom2-SH. We used the following SNR equation to modulate parameter settings and evaluate the reconstruction qualities:

⎛ ⎜ ⎜ ⎜ ⎝ J j f phantom-j -fphantom 2 J j (f phantom-j -f j ) 2 ⎞ ⎟ ⎟ ⎟ ⎠ (16)
f j , f phantom-j denote the intensity of the pixel j in the reconstructed image and the original phantom image respectively. fphantom provides the mean of all the pixels in the original phantom image.

Fig. 4 shows the reconstruction results for the two sinograms simulated with a low dose scan (Phantom1-SL and Phantom2-SL). The visual comparison suggests that FBP1 method (with Ramp filter) suffers from severe noise effects and streak artifacts. Although FBP2 method (with Hanning filter) outperforms FBP1 method in noise suppression, many artifacts remains. Huber prior and MRP reconstructions tend to produce blocky piecewise smooth regions in the background and staircase effects around the edges. The straight lines introduced in the Shepp-logan phantom appear smoothed especially with the MRP method, the thinner lines being drowned out in the background.

With TV prior, the results are improved, edges are better preserved even if we can observe a light staircase effect around the edges. Finally the best results are given with the PSM prior for which much less oversmoothing and staircase effects are observable. with a higher level of count. We can observe the same effects in a lower proportion since the magnitude of the noise is less important in the original measurements. FBP1 and FBP2 still produce some streak artifacts, Huber and MRP priors provide a smoothed image excepted that edges are better preserved in the image reconstructed with MRP prior, yielding thus a higher level contrast than in the former one. Smaller objects such as the straight lines appears nevertheless relatively smoothed. TV prior performs better in preserving edges and PSM prior still produce the best result both in term of smoothing and edge preservation. In TABLE II, we list the calculated SNR of all the reconstructions with the simulated data. Obviously, the SNR comparisons indicate that Bayesian reconstructions can produce images with higher SNR than analytic FBP reconstructions, and the reconstructed images using the proposed prior can obtain the highest SNR among all the reconstruction methods.

(a1) (b1) (c1) (d1) (e1) (f1) (a2) (b2) (c2) (d2) (e2) (f2)

B. Reconstruction Results on Clinical Data

Real data came from a pelvis scan of a 56 year aged patient. This dataset was acquired on a monoenergetic Siemens Somatom volume zoom 1 detector scanner using the following protocol: 100 mA, 120 KV, pixel size:1mm × 1mm, slice thickness: 1 mm, detector cell spacing :

Methods

Parameter setting

TV prior rec.

β = 5 × 10 -3 , = 1 × 10 -10
Huber prior rec. normalized, corrected for attenuation and scatter and rebinned to 2D using Fourier rebinning. The d i factor in [START_REF] Bertero | Ill posed problems in early vision[END_REF] was set to 1. As for the simulated data, we performed a set of experiments to select the optimal parameters for each method and asked a radiologist to evaluate the reconstruction quality. Listed in TABLE III, the parameters were thus chosen according to a visual inspection and based on the capability of the method to preserve the original anatomical structures.

β = 4 × 10 -1 , γ = 0.95 MRP rec. β = 8 PSM prior rec. β = 1.5 × 10 -2 , λ = 7.94
Fig. 6 illustrates the reconstructions on real clinical data. All the images are displayed using a window level and width respectively equal to 60 and 250 HU (Housfield units). FBP1 (with Ramp filter) produces similar effects to those met on the synthetic phantoms and FBP2 (with Hanning filter, Fig. 6(b)) provides an improved image with lower streak artifacts. For the other methods, the differences locate in the details. We can thus observe in a lower proportion, blocky piecewise smooth regions in the background with undesirable staircase effects around the edges in the reconstructed images (Fig. 6(c),(d)and (e)) from TV prior, Huber prior and MRP reconstructions. However, we can still notice blocky piecewise smooth regions with the TV prior, Huber prior and MRP reconstructions in Fig. 6(c),(d)and (e). TV prior reconstruction provides a smoothing background with a loss of resolution on the small structures while preserving the sharp transitions in the images (see transitions between bones and soft tissues).

Notably, in the proposed reconstruction of pelvis data in Fig. 6(f), compared to other results, the bones (the brightest part of the image), scrotum, muscles(the regions around the bones) and fat tissue(the outer region beneath the skin surface) can be better identified with clearer detail preservation (the parts pointed by red arrows).

To complete this qualitative analysis, each method was evaluated in term of noise suppression, artifact suppression, contrast preservation and global image quality, using a five-point scale (1 = unacceptable, 2 = substandard, 3 = acceptable, 4 = above average, 5 = excellent). This evaluation was realized by three radiologists working in the department of radiology at the Nanjing Hospital in China (X.D.Y., with 15 years of experience, X.H.Y., with 8 years of experience and Y.M.D., with 5 years of experience). We asked them to independently assess the six randomized reconstructed images of the Figure 6. Un-stationary shapes that had a passive effects on image quality and included streak and blocky artifacts, were considered as artifacts. TABLE IV reports the quality scores of the different methods under the form of a mean score calculated on the responses of the three radiologists. The global mean score computed on the four criteria, for each method, hauls up the PSM method as the best one, followed by the the TV method. If we now look at each score, we can observe that the proposed PSM prior reconstruction is always superior to other reconstructions.

(a) (b) (c) (d) (f) (e)

C. Computation Time Comparison

We measured then the acceleration factor, which is defined as the reduction factor in reconstruction time between the CUDA accelerated and the unaccelerated reconstruction algorithms. This study was carried out on all the algorithms and for the data Phantom1-SL using the same parameters as defined in above phantom experiment.

The computation time is given in TABLE V for each algorithm in each situation: accelerated vs unaccelerated. For reference, the FBP reconstruction takes 21.89 CPU seconds. Without surprise, statistical methods appears much more expensive in computation time with a CPU time of about 4300 CPU seconds for the PSM prior method, and 900 CPU seconds for other prior methods. The GPU-based CUDA acceleration significantly reduces the overall computation time and allows a speed up factor of respectively 14 (for TV prior, Huber prior and MRP reconstructions) and 30

(for the proposed PSM reconstruction). The reason that the PSM prior reconstruction can obtain higher acceleration than other Bayesian reconstructions is that the prior involved calculation has a higher proportion in the PSM prior reconstruction. This parallelization was here limited to the weight calculation, and we would try to realize global acceleration of the algorithm (including forward/backward projections, dot multiplication between large matrixes [START_REF] Fessler | Aspire 3.0 user's guide: A sparse reconstruction library[END_REF][START_REF] Scherl | Fast GPU-based CT reconstruction using the common unified device architecture (CUDA)[END_REF])

D. Initial Conditions and Monotonicity for The Joint Image/Weight Reconstruction Algorithm

We carry out the PSM prior reconstruction on the Phantom1-SL data with the same parameters as previous experiment and a total number of iterations equal to 1000. Fig. 7 that using a FBP reconstructed initial image to initialize the reconstruction process, lead to the formation of streak artifacts in the final reconstructed image while a simple uniform initial image provides a good result.

We analyzed then the monotonicity of the PSM algorithm by calculating the posterior function energy ψ β (f, w, η|g)

and the SNR with respect to the iteration numbers. Fig. 8 depicts their evolutions according to the iteration numbers.

It shows that we have a monotonic increase of the SNR and posterior function energy with the iteration progression.

We can thus conclude that the convergence of the posterior energy is ensured in the iteration of the PSM algorithm. Fig. 9 provides the iterated weight maps for the two specified pixels in the image. It shows that the algorithm can achieve stable weight updates as iteration grows.

E. Study on Size Settings of N and n

We study here the impact of the size setting of N and n on the reconstruction qualities when using the proposed PSM prior. We performed the PSM reconstruction while making vary the sizes of N and n. We limited the size of n to 7×7 because a larger n might fail to present discriminative characteristics between the local structures and noise. The above CUDA accelerated algorithm was applied to perform reconstruction with data Phantom1-SL.

Parameters β and λ were changed according to patch sizes based on the rule of SNR minimization. The best value for β was found to be 0.065, and λ was set to 8.32, 1.53 and 0.17 for n size respectively equals to 7×7, 3×3 and 1×1. Table VI provides for different couples (N ,n), the calculated SNRs of the reconstructed images and the CPU times. When n is set to 1, the patch reduces to one single pixel and we recover the conventional pairewise TV-like priors. Fig. 10 illustrates the reconstructed images with zoomed regions of interests (in the red line boxes).

We can see that when a small size of n is associated with a large size of N , the reconstruction provides blocky artifacts or staircase effects around the regions of interest. These effects decrease and the intensity transitions in the \ 1 st iter 21 st iter 41 st iter 81 st iter 101 st iter 121 st iter 141 st iter 161 st iter 181 st iter 201 st iter 221 st iter 241 st iter 261 st iter 281 st iter 301 st iter 321 st iter 341 st iter 361 st iter 381 st iter 401 st iter 421 st iter 441 st iter 461 st iter 481 st iter 501 st iter 521 st iter 541 st iter 561 st iter 581 st iter 601 st iter 621 st iter 641 st iter 661 st iter 681 st iter 701 st iter 721 st iter 741 st iter 761 st iter 781 st iter 801 st iter 821 st iter 841 st iter 861 st iter 881 st iter 901 st iter 921 st iter 941 st iter 961 st iter 981 st iter 1000 th i t e r zoomed regions are better restored when the size of n increases. For this kind of images, the best SNR is given for the highest values of N and n. As expected, we can note in TABLE VI that more CPU time costs are needed for the reconstructions using the proposed prior with larger N and n. We observe in Fig. 10

IV. CONCLUSIONS AND FUTURE WORK PLAN

We proposed in this paper a novel PSM prior which penalizes the distances between neighboring patches through a constrained mixture prior model. A convergent joint image/weight update algorithm, which estimates sequentially the image and the weights, has been also proposed to avoid the heuristically OSL weight determination. This algorithm has been applied in CT reconstructions with both synthetic phantom and real scan data. The quantitative and subjective analysis demonstrated the efficiency of the regularization in both suppressing noise and preserving detail as its supremacy over the classical edge-preserving TV prior, Huber prior and MRP algorithms. The uniform initial condition has shown to provide a better reconstruction result than the FBP initial condition and monotonic update of both image and weight has been demonstrated. The study on the sizes of the neighborhoods N and n has shown improved performances for larger values N and n.

The reconstruction using this new prior needs blockwise computations of the distances between patches over a large region, which leads to a computational burden. To overcome this problem, we suggested a graphic hardware acceleration to calculate the patch distance in a parallel way using the CUDA framework. A speed up of factor 30 has thus been achieved. This parallelization was here limited to the weight calculation. Next stage will be devoted to further accelerating the whole reconstruction by parallelizing the projection and image updates [START_REF] Scherl | Fast GPU-based CT reconstruction using the common unified device architecture (CUDA)[END_REF][START_REF] Knaup | GPU-based parallel-beam and cone-beam forward-and backprojection using CUDA[END_REF], developing more robust parameter estimations, and evaluating the proposed algorithm on more clinical CT data and other imaging modalities. is associated with the couple (N ,n) having the same index in TABLE VI. Let's note that it is not easy to distinguish the differences between the images reconstructed with a size of N equal to 11×11 and 13×13 especially when n is equal to 7×7.

Fig. 1 .

 1 Fig. 1. (a),pixel-wise relations in pixel discontinuity (difference) calculation for traditional prior models; (b), Patch-wise relations in patch discontinuity (difference) calculation in the proposed PSM prior model.

  Allocate global memory for variables in the device;For all pixel j with its neighborhood N j

Fig. 2 .

 2 Fig. 2. CUDA acceleration of the computation of the UP SM (f j , ŵ, η| f ) and ÜP SM (f j , ŵ, η| f ) in the PSM prior reconstruction: Serial vs parallel pixel-wise computation

5 )

 5 MRP rec. : Bayesian-based Median Root Prior reconstruction algorithm [10]; 6) PSM prior rec. : Above joint MAP algorithm applying the PSM prior with 7×7 patch n and 11×11 neighborhood N .

Fig. 3 . 5 PSM prior reconstruction β = 6 . 5 ×

 3565 Fig. 3. Simulated phantom images.

Fig 5

 5 Fig 5 provides the reconstruction results on the sinogram Phantom2-SL for which the sinogram has been simulated

Fig. 4 .Fig. 5 .

 45 Fig. 4. Low dose reconstructions for the simulated data Phantom1-SL (a1-f1) and Phantom2-SL (a2-f2). (a1) and (a2), FBP1 reconstruction using Ramp filter; (b1) and (b2): FBP2 reconstruction using the Hanning filter with cutoff at 80% Nyquist frequency; (c1) and (c2): TV prior reconstruction; (d1) and (d2): Huber prior reconstruction; (e1) and (e2): MRP reconstruction; (f1) and (f2): The proposed PSM prior reconstructions.

Fig. 6 .

 6 Fig. 6. FBP and Bayesian reconstructions using different priors for clinical pelvis data. (a), FBP1 reconstruction using Ramp filter; (b), FBP2 reconstruction using the Hanning filter with cutoff at 80% Nyquist frequency; (c), TV prior reconstruction; (d), Huber prior reconstruction; (e), MRP reconstruction; (f), the proposed PSM prior reconstruction.

FBP1

Fig. 7 .

 7 Fig. 7. Impact of different initial conditions on the reconstruction results. (a), FBP reconstructed initial image; (b), unitary initial image.

Fig. 8 .

 8 Fig. 8. Monotonicity Study : SNR (upper left:1-1000 iterations;upper right:101-1000 iterations) and posterior function energy (lower left:1-1000 iterations;lower right:101-1000 iterations) with respect to iteration number.

Fig. 9 .

 9 Fig. 9. Phantom1-SL reconstruction using the PSM prior algorithm -The total number of iteration was set to 1000 : weight evolution map during iterations for two points (Point 1 and Point 2).

Fig. 10 .

 10 Fig. 10. Reconstruction results according to the sizes of N and n. Zoomed regions of interest are also illustrated in red line boxes. Each image

TABLE I PARAMETER

 I

SETTING IN THE EXPERIMENTS ON SIMULATED DATA SN R = 10 log 10

TABLE III PARAMETER

 III SETTING IN THE EXPERIMENTS ON CLINICAL DATA

  (a) and (b) illustrate this situation with a 1-value unitary image and an image previously built from a FBP reconstruction with Ramp filter, respectively. Results show

		TV prior rec.	Huber prior rec.	MRP rec.	PSM prior rec.
	CPU computation time for accelerated algorithms	66.54 sec	67.09 sec	55.74 sec	1474.42 sec
	CPU computation time for Unaccelerated algorithms	909.62 sec	907.23 sec	835.57 sec	43090.83 sec
	Acceleration factor	13	13	15	30

TABLE V CPU

 V TIMES (SECONDS) FOR EACH BAYESIAN BASED RECONSTRUCTION ALGORITHM

  and TABLE VI that no significantly distinct visual and SNR differences are made between 11×11 N and larger 13×13 N . This verifies our choice of 11×11 N and 7×7 n in this work.

	N and n	SNR	CPU times	N and n	SNR	CPU times	N and n	SNR	CPU times
	(a), N 3×3 n 1×1	18.88	158.98 sec	(b), N 3×3 n 3×3	19.62	378.59 sec	(c), N 3×3 n 7×7	21.76	679.34 sec
	(d), N 7×7 n 1×1	18.95	119.12 sec	(e), N 7×7 n 3×3	21.59	437.30 sec	(f ), N 7×7 n 7×7	22.78	542.09 sec
	(g), N 11×11 n 1×1	19.23	198.59 sec	(h), N 11×11 n 3×3	21.05	632.42 sec	(i), N 11×11 n 7×7	23.75	944.42 sec
	(j), N 13×13 n 1×1	19.19	491.54 sec	(k), N 13×13 n 3×3	21.69	837.67 sec	(l), N 13×13 n 7×7	24.06	1240.31 sec

TABLE VI PHANTOM1

 VI -SL RECONSTRUCTION USING THE PSM PRIOR ALGORITHM: SNR AND CPU COMPUTATION TIME (SECONDS) ACCORDING TO THE SIZES OF THE PATCH n AND THE NEIGHBORHOOD N .

ACKNOWLEDGEMENT

This research was supported by National Basic Research Program of China under grant (2010CB732503),National Natural Science Foundation under grant (81000636), and the Project supported by Natural Science Foundation of Jiangsu Province (BK2009012).


Mathias John

Christian Tominski

Heidrun Schumann

Visual and Analytical Extensions for the Table Lens

Keywords: Information Visualization, Two-Tone Pseudo Coloring, Table Lens

Many visualization approaches teach us that ease of use is the key to effective visual data analysis. The Table Lens is an excellent example of a simple, yet expressive visual method that can help in analyzing even larger volumes of data.

, we obtain visual representations that can communicate larger volumes of data while still maintaining precision. Secondly, we propose to integrate a data analysis step that implements a hybrid clustering based on self-organizing maps and hierarchical clustering. The analysis step helps to extract and communicate complementary structural information about the data and also serves to drive interactive information drill-down.

INTRODUCTION

Today's datasets are growing in size and complexity. The challenge of visualizing larger datasets is more current than ever before. In order to extract global information from such data efficiently, a compact representation that provides overview and structural insight is needed. On the other hand, it is important to display data values very precisely, because otherwise detail information gets lost. Combining overview and detail is not always an easy task -for both visualization designers as well as users. [START_REF]Readings in Information Visualization: Using Vision to Think[END_REF] We think that simplicity is a requirement for user acceptance and a wider range of applicability. In this work, we adhere to simple, yet expressive visualization methods. In particular, we present two extensions of the well-known Table Lens, [START_REF] Rao | The Table Lens: Merging Graphical and Symbolic Representations in an Interactive Focus + Context Visualization for Tabular Information[END_REF] which is an elegant approach potentially useful in many applications domains. The first extension concerns the integration of Two-Tone Pseudo Coloring (TTPC) [START_REF] Saito | Two-Tone Pseudo Coloring: Compact Visualization for One-Dimensional Data[END_REF] into the Table Lens. While the Table Lens is particularly suited to provide overview and detail for relational datasets with a large number of data items, TTPC has been developed to provide compact overviews of larger value ranges. At the same time, TTPC increases the precision that can be achieved when reading details from a visualization. Consequently, the goal is to combine the strengths of Table Lens and Two-Tone Pseudo Coloring. The second extension concerns the integration of analytical methods into the Table Lens. With the help of an additional analysis step, we are able to rearrange table rows according to clustering information. Our hybrid clustering approach can be seen as a complementary feature to the sorting functionality of the original Table Lens to gain structural insight into the data being visualized.

In Section 2, we provide brief background information to lay the ground for extending the Table Lens by TTPC in Section 3. There, we will also describe a heuristic to improve readability of TTPC-based visualization significantly. In Section 4, we will explain how to integrate self-organizing maps (SOM) and hierarchical clustering in order to generate visual representations that convey structural information as well. Section 5 is devoted to further extensions and visual examples of the proposed approaches. A summary and an outlook on future work will be given in Section 6.

BACKGROUND

To facilitate understanding of later sections, we will provide a brief view on the methods we use in this work. We will take a look at the role of color and consider the importance of overview+detail and precision in visualization.

Two-Tone Pseudo Coloring

In his Semiology of Graphics Jacques Bertin describes color as an important visual variable. However, proper use of color is a requirement to fully exploit the potential of color coding. [START_REF] Brewer | Color Use Guidelines for Mapping and Visualization[END_REF] Still, many visual representations make use of rather basic color coding schemes, even though better approaches are available in the literature. [START_REF] Bergman | A Rule-based Tool for Assisting Colormap Selection[END_REF][START_REF] Tominski | Enhancing Visual Exploration by Appropriate Color Coding[END_REF] Presumed such methods are used, color coding is an effective means to visualize qualitative and quantitative data. However, it is not always easy for users to read data values precisely when looking at color coded visual representations. Fig. 1 shows examples of discrete and continuous color coding. In the former case, only ranges of colors, and hence, only ranges of data values can be identified. In contrast to that, continuous color coding is based on a scale of smooth color transitions. However, users tend to perceive a sequence of colors as a set of discrete ones. [START_REF] Ware | Information Visualization: Perception for Design[END_REF] So, even with very smooth transitions, values cannot be read very accurately.

Recently, Two-Tone Pseudo Coloring (TTPC) has been introduced as a technique to facilitate visualization of univariate quantitative data. [START_REF] Saito | Two-Tone Pseudo Coloring: Compact Visualization for One-Dimensional Data[END_REF] TTPC uses a discrete color scale that consists of only few colors (less than eight). Whereas classic color coding (irrespective of whether it is discrete or continuous) maps each data value to one color, TTPC uses two colors, which are adjacent in the color scale. These two colors are interpreted in two steps (see Fig. 1). First, the two colors guide users to a particular interval in the value range. If users find that interval interesting, they go into detail: The proportion of use of both colors encodes the precise data value. This two-step interpretation is the basis for overview+detail. While the first step relies on the perception of color, which facilitates overview, the second step is based on the human capabilities to judge lengths, which provides for sufficient detail. Compared to other approaches, TTPC generates very compact, yet precise visual representations. As seen in Fig. 1 line charts require much more screen space than a TTPC visualization.

However, in its original form, TTPC is suitable for univariate data only. To visualize multivariate data, other methods must be applied, as we will see in the next section.

Table Lens

Multivariate data (i.e., data with more than one variable) are usually modeled as a data table, where columns accommodate variables and rows represent data items. Even though common spreadsheet representations are very intuitive, they cannot be used to display datasets with a larger number of data items at a glance. The Table Lens 2 addresses this problem. It reduces the required amount of display space by shrinking the height of table rows: Thin bars are used to encode data values. The length or position of a bar serves as indication for the particular data value. To facilitate precise reading of data values, user can expand rows to see data values in textual form. In this sense, the Table Lens supports overview (representation as thin bars) and provides details on demand (value print out). In the original Table Lens, overview and detail concern data items (i.e., table rows). In the next section, we will extend the Table Lens by TTPC to achieve more compact representations with respect to variables (i.e., table columns). TTPC encoding results in representations that are based on both perception of color and length. With such a two-channel encoding it is easy to maintain an overview, even when table columns are narrow in the case that many attributes are being displayed.

In order to assist users in finding hidden relationships, it makes sense to support rearrangement of table rows. Users of the original Table Lens can sort the data on demand on a column-by-column basis. As we will see in Section 4, providing a hybrid clustering mechanism complementary to standard sorting functionality is helpful to reveal structural information.

All in all, we think that integrating novel visual representation and analytical methods into the Table Lens bears much potential in terms of increasing its usefulness and applicability. With these thoughts, we are in line with other researchers in the field. [START_REF] Telea | Combining Extended Table Lens and Treemap Techniques for Visualizing Tabular Data[END_REF]

TWO-TONE PSEUDO COLORED TABLE LENS

In what follows, we introduce a way to join the advantages of The intuitive overview+detail concept allows users to change the width of columns and to adjust the height of rows (geometric zoom), and it is possible to expand table rows instantly to see print outs of the data in textual form.

Nonetheless, integrating TTPC into the Table Lens is not without problems. First of all, we have to support users in learning the new way of representing data as TTPC-coded columns. Second, since TTPC involves color, we have to take care to use color scales that integrate smoothly into the data exploration process. In the next paragraphs, we will address these problems.

Learning the TTPC Table Lens

To support the process of learning TTPC Table Lens representations we build a bridge between the classic Table Lens and our approach This is accomplished by so to say interpolating between both techniques. Interestingly, this is possible, because the original Table Lens is a special case of a TTPC Table Lens where n, the number of color segments, equals one. We allow users to set n to any value between one and eight (see Fig. 2). Precision is increased for high values of n, but then users need some experience to interpret the TTPC Table Lens. By setting n to low values, novel users get easy access to a TTPC Table Lens and can get accustomed to the new presentation technique.

Effective Use of Color

The choice of colors used as well as the segmentation of the color scale have major impact on the effectiveness of the TTPC Table Lens. Consider a color scale with four segments like black : [0.3-0.99], red : [0.99-1.95], yellow : [1.95 -2.91], green : [2.91 -3.87]. Users that are not familiar with wavelengths of colors will hardly understand the order in that color scale. Furthermore, switching colors at odd values impedes reading intermediate values from a TTPC representation. Which colors are better suited and how segmentations like the previous one can be automatically transformed into easier to understand ones will be shown next.

Color scales To visualize quantitative variables with color, it makes sense to map data values to the saturation component. [START_REF] Brewer | Color Use Guidelines for Mapping and Visualization[END_REF] When visualizing multiple quantitative variables it is furthermore very helpful to assign separate color scales to each variable to communicate that each variable has its own value range. Variables are then easier to recognize, which is particularly useful in tabular representations that allow for column reordering. Because we want to use saturation to visualize data values, our only choice to generate separate color scales for the variables is to vary the hue component. We carefully chose eight color scales provided by the ColorBrewer 9 and applied those in a way such that adjacent columns of the Table Lens can be distinguished effortlessly. However, if the TTPC Table Lens is used in specific application scenarios, user preferences or application background (e.g., context-specific color associations or per se correlated variables) may raise the need to use color scales tailored to the requirements at hand.

Color scale segmentation As already indicated, carelessly set color segment borders are a problem. To generate more intuitive color scales, we propose a novel color scale segmentation heuristic. The basic idea is to spread the value range mapped to a color scale such that the subdivision into color segments results in transition points that are easier to interpret. How the heuristic method works can be best explained by our introducing example, where the value range was [0.3-3.87]. Assumed the number of color segments is n = 4, then a simplistic subdivision of the color scale would result in the following awkward transition points (i.e., segment borders): p 0 = 0.3, p 1 = 0.99, p 2 = 1.95, p 3 = 2.91, p 4 = 3.87. In contrast to that, our heuristic method does the following: In the first step, the original value range is scaled to a fixed range of 10 to 100. This step is necessary, because our heuristic has to cope with arbitrary magnitudes. In the second step, we set the maximum to a more intuitive value by increasing the scaled maximum to the next multiple of 10. The third step computes the distance between transition points. In step four, we consider the fact that people commonly tend to cut things into halves when they have to make estimations. To support this habit, we increase the computed distance to the next multiple of 2 and determine a new minimum value. Finally, the fifth step is to bring the values back to their original magnitude. Applied to the previously mentioned example, our heuristic spreads the value range such that the transition points read: p 0 = 0, p 1 = 1, p 2 = 2, p 3 = 3, p 4 = 4, which is much easier to interpret compared to straightforward color segmentation.

We tested our method with various quantitative variables and recognized that it usually helps to increase the readability of the TTPC Table Lens. Even though we believe that our approach can also be applied to other color-based visualization techniques, a heuristic will never be able to foresee all possible requirements and constraints that may be implied by certain application scenarios. Therefore, we allow users to adjust the color scale segmentation interactively. Fig. 3 shows a comparison of a simplistic segmentation, the result of our heuristic method, and a segmentation adjusted by a user. It is quite obvious, that the heuristic has advantages over the simplistic method. Whether the heuristic can compete with user-defined segmentations depends on the application context and tasks at hand.

HYBRID DATA CLUSTERING

So far, we have seen that improved color scales are very useful in terms of readability. A further aspect must be taken into account to make our approach applicable for arbitrary datasets: Stronger oscillation in the data hinder users in spotting relations between variables, and thus degrades efficiency. This is due to the larger number of color switches present in the columns of the tabular representation (see Fig. 4(a)). A color switch is considered if vertically adjacent table cells do not share common colors. In such cases, the cognitive effort is higher, because each cell must be interpreted separately. On the other hand, if adjacent cells show the same colors, the overview (step one in interpreting TTPC) is the same for these cells, only their precise values vary (detail). Users can derive data values more easily. We call this cell-to-cell coherence.

To utilize cell-to-cell coherence, that is to reduce the number of color switches, a mechanism is needed that assigns similar data items (in terms of Euclidean distance) to rows that are close in the visual representation. A first possibility, which is also provided in the original Table Lens, is to sort data values. Sorting data successively with respect to user-chosen columns opens up the chance to arrange similar table rows close together. However, chances decrease with a growing number of variables. Moreover, this way of sorting is time consuming in terms of required interaction steps.

Clustering methods for multivariate data are a promising alternative. The result of clustering is a set of clusters, which contain similar data items. Clustering yields two advantages: First, knowing clusters and their associated data items makes it easy for us to rearrange table rows to avoid strong data oscillation, and second, clustering methods crystallize structure within the data.

Because our interests regard larger datasets, we have to choose a clustering method that generates results within reasonable time. Self-organizing maps (SOM) are well-suited for finding clusters in multivariate data. [START_REF] Kreuseler | A Flexible Approach for Visual Data Mining[END_REF] Yet, their computational complexity does not fully fit our needs. Therefore, we opted for implementing a hybrid clustering that integrates SOM clustering and hierarchical clustering. Our hybrid clustering performs the following two steps:

1. Coarse SOM clustering to generate a reasonable number of SOM clusters, 2. Refine each of the SOM clusters with hierarchical clustering.

We will take a closer look at our hybrid clustering method and provide some test results in the following paragraphs. Step one -SOM clustering Self-organizing maps 11 are artificial neuronal networks whose neurons are typically arranged on a rectangular lattice. For our purposes we use a 1D variant of SOM, for which the lattice is reduced to a chain. This enables us to use the clustering result directly to rearrange rows in the Table Lens. In a SOM, each neuron is associated with an m-dimensional neuron vector, where m is the number of variables of the dataset. After a network training with randomly chosen training vectors, each neuron vector can be considered the representative of a cluster of similar data vectors. Moreover, adjacent cluster representatives (in terms of their position in the chain of neurons) show similarities. Hence, we have achieved our goal to bring similar data vectors close to each other.

Although SOM are straightforward to implement, the problem is to estimate the number of neurons (i.e., clusters) required to generate results that are balanced in terms of clustering quality and computation time. The proper number of neurons is not known in the first place and depends heavily on the given data. Therefore, we take a closer look at the influence of the number of neurons used. Notice that once the number of neurons is known, the number of training iterations can be determined by the following rule of thumb: [START_REF] Kohonen | Self-organizing Maps[END_REF] training iterations ≈ 500 • number of neurons.

Besides extracting structural information, for us the primary goal of clustering is to reduce data oscillation to improve the visual representation. So, roughly speaking, clustering quality can be related to the number of color switches. The occurrence of a single color switch within a certain area is not really a problem, yet, the sequence of many is. Therefore we introduce the term crucial color switch. It describes a color switch whose distance to another switch located above in the spreadsheet is less than a certain number of rows. The maximum distance for a crucial color switch is approximately five rows. Color switches with larger distances do not seem to interfere with value reading. Yet, perceptual tests may yield a slightly different threshold. To evaluate clustering results, we created a metric that counts all crucial color switches and relates them to the total number of color switches. For a dataset with m variables and k data vectors, we compute our metric M (c) for each column c of the TTPC Table Lens as follows:

M (c) = k i=1 switch(i, c) • crucial(i, c) n -1 (1) switch(i, c) = 1 | S(v i-1,c ) = S(v i,c ) 0 | else (2) crucial(i, c) = 1 | ∃j : i < j ≤ i + d ∧ S(v j,c ) = S(v i,c ) 0 | else (3) d ≈ 5 (4)
where v i,c denotes the data value in row i and column c, and S(v i,c ) is its associated color segment. To get a single value for clustering quality, we compute the average of M (c) of all m columns. We tested pure SOM clustering quality for different numbers of neurons on datasets of different sizes (see Tab. 1). Since SOM is a nondeterministic method, we conducted the evaluation five times and averaged the final result. Tab. 2 shows that reducing the occurrences of crucial color switches means to significantly increase the number of neurons. That is a problem, because the training of large SOM networks (e.g., 1000 neurons) takes an unreasonable amount of time. This is where the second step of our hybrid clustering comes into the picture. Step two -Hierarchical clustering Hierarchical methods introduce a step by step procedure to cluster data vectors. [START_REF] Kaufmann | Finding Groups in Data -An Introduction to Cluster Analysis[END_REF] The result is a hierarchical structure with several cluster levels, where the root of the hierarchy is a representative for the whole dataset and the leaves are the original data vectors. The hierarchical structure is very useful from a visualization perspective, because it communicates the process of generating clusters very well.

climate
Hierarchical clustering methods usually require a dissimilarity matrix, which describes the distances between any two clusters. This matrix can become huge and its computation may be too expensive if we consider larger datasets for visualization. Therefore, we do not apply hierarchical clustering to the entire dataset, but only to the generated SOM clusters. We tested this hybrid approach and derived the following interesting conclusion from Tab. 3. Since hierarchical clustering refines SOM clusters, the number of neurons used for the SOM step is rather independent from clustering quality. This independence allows us to reduce the number of neurons without compromising quality and, thus, to accelerate the hybrid clustering. We propose the following rationale: First, we can reduce the number of neurons required for the SOM step. Second, we avoid huge dissimilarity matrices, since the SOM clusters to be refined contain only a limited number of data vectors.

number of neurons = √ k ,
We can conclude that by integrating our hybrid clustering into the Table Lens strong oscillation in the data can be reduced significantly. Fig. 4(b) shows the visual effect: The overview is much easier to grasp and precise values can be read with less effort. Moreover, hidden multivariate relationships within the data can be revealed, which is an advantage over plain sorting functionality. It is worth mentioning that even for larger datasets, the hybrid approach is still capable of reflecting the underlying data characteristics very well, without putting too much burden neither in terms of computation nor in terms of memory usage.

FURTHER EXTENSIONS & EXAMPLES

Further Extensions

In order to widen the range of applicability of the TTPC Table Lens, we considered two further aspects:

• Communication of underlying data characteristics

• Reproducibility of visualization results

The hybrid clustering process induces a hierarchical structure that is very well suited to communicate underlying data characteristics. We make use of that fact by providing an additional Icicle plot to visualize the cluster hierarchy. An Icicle Plot requires less space than other designs [START_REF] Barlow | A Comparison of 2-D Visualizations of Hierarchies[END_REF] and can be integrated easily by inserting an additional column into the Table Lens (see Fig. 4(b)). Table Lens and Icicle Plot are linked, so that we achieve a very tight coupling of data representation and visualization of underlying structure. We also integrated homogeneity information derived from the hierarchical clustering step. This information is valuable when grading similarity between cluster elements. Lastly, the Icicle Plot drives a semantic zoom facility: Each element in the Icicle Plot and its header row can be clicked to expand single clusters or entire levels of the cluster hierarchy respectively. Thus, starting with a completely collapsed Icicle Plot, users can interactively drill-down to areas of interest; the Table Lens is updated accordingly (see Fig. 4

(c)).

A second extension regards reproducibility. We think that being able to reproduce visualization results is a key to success of any visual method. Since our approach involves a nondeterministic analysis step (SOM clustering), reproducing a particular result can be an awful task. Therefore, we incorporated a mechanism that enables users to save the current state of the visualization, including color scale, clustering results, current state of the Icicle Plot. That allows users to come back later, either to revisit the visualization result or to continue with their analysis task. The saving functionality is also useful when comparing visualization results that were generated using different parameterizations.

Examples

In the previous sections, several extensions to the Table Lens approach were introduced. The potential of the extended Table Lens evolves from its different features, including overview+detail concepts, functionality to increase readability (i.e., conscious use of color scales and color segmentation), integration of a hybrid clustering linked to an interactive drill-down facility, and mechanisms to support reproducibility. We implemented our approach in a tool [START_REF] John | Prototype: Two-Tone Pseudo Colored TableLens[END_REF] to get an impression of how well the features can be applied to certain problems in practice. In what follows, three different examples will be explained. Let us sort the data with respect to petal length (4th column). At first view a discontinuity in the middle of that column can be recognized (see Fig. 5(a)). The area describes a value step which is silhouetted against a smooth value development. Supported by TTPC overview and an appropriate color scale this area can be found immediately. Using Table Lens' mechanism to expand rows of interest, we can easily confirm a numeric gap (see Fig. 5(b)), meaning that the dataset holds no Iris exemplar with a petal length between 1.9 cm and 3.0 cm.

Countries dataset "Countries" is a dataset describing the census of 110 countries in the early 90s. As can be seen in Fig. 6(a) the representation of the original dataset looks rather abstruse. Therefore, we applied hybrid clustering using ten neurons as suggested by our tool. The resulting visualization looks smoother and is easier to interpret. By taking a closer look, one can find an interesting cluster that contains exclusively former east block countries (see Fig. 6(b)). The Table Lens view clearly shows, that only these countries share an exclusive property, that is the combination of high literacy rate with low GDP, low birth rate and rather high death rate (see Fig. Climate dataset We now consider a larger dataset related to climate research. [START_REF] Nocke | Methods for the Visualization of Clustered Climate Data[END_REF] It contains 40419 data vectors about daily weather conditions, including for instance temperature, air pressure, precipitation, sunshine duration etc., at the Potsdam observatory. By looking at Fig. 7(a) one can see a smooth value development in most of the columns, which is due to the fact that weather usually changes gradually. From looking at just the first 365 data vectors, we may guess that a correlation exists between temperature and vapor content (4th and 5th column). However, since only a fraction of the data is visible, we cannot be confident. What we could do in such a case is to scroll all the way through the data to confirm our guess. Yet, that would be rather inconvenient in terms of navigation, since we are dealing with a rather large dataset. The alternative is to apply clustering to the data first, and then to collapse all the elements of the Icicle Plot by a single click. This results in an overview that shows average values for each cluster (see Fig. 7(b)). This overview is much easier to scroll, and by doing so, user can see that there is indeed a correlation between temperature and vapor content throughout the whole dataset. Hypotheses about possible correlations can be confirmed in detail later on by applying the drill-down mechanism introduced in Section 5.1. In that case, the free space on the right hand side of Fig. 7(b) will again be filled with the structural view of expanded clusters.

Expert Experience Report

In this section, we report on an initial attempt to evaluate our approach by informal expert interviews. In particular, we gave our tool to the expert biologist of our graduate school, Ben. Ben is interested in our approach since he needs visual support for his data analysis tasks. In one of his experiments he tested the genetic reaction of immune cells on two different parameters of environmental influence. The resulting data spanned a table of 14 columns (regulation values) and over 15.000 rows (genes). After a brief introduction to our tool Ben needed a little time to get used to the Two-Tone representation. Therefore, he modified the number of color segments, starting with the TableLens representation (one interval) and then stepwise increasing it to five. Once feeling confident, Ben applied hybrid clustering to the data and entirely collapsed the Icicle Plot for a general overview. Ben was mostly interested in genes that are specifically upregulated for one parameter. As it turned out, the SOM clustering generated exactly two clusters with genes that fulfill this requirement. To gain more confidence Ben expanded the Icicle Plot and concentrated on the respective values. He pointed out two very high values as being incorrect. He adjusted the color scale segment borders such that there was one large segment for the incorrect outliers and several smaller ones for the other values, which thus could be read off more precisely.

In Ben's opinion the most helpful feature of the extended TableLens is that many rows but also many columns can be explored at the same time while still being able to read off specific values. He said, that the clustering facilitates the search for interesting data elements. Additionally, he suggested that functionality for extracting table parts into new tables should be introduced. Thus, he could explore and sort different table parts separately.

SUMMARY & FUTURE WORK

Our work shows that the Table Lens can be subject to visual and analytical extensions improving the usefulness of the original approach. We proposed two major extensions. The first is the integration of Two-Tone Pseudo Coloring (TTPC) as an alternative visual encoding. In that context, we presented a heuristic for color scale segmentation that enables users to read precise values more easily. This heuristic can also be applied to other color-based visualizations. The second extension is a hybrid clustering algorithm that helps to reduce strong oscillation in the representation and to communicate the underlying structure of the data. Strong oscillation is reduced by automatic rearrangement of data vectors with respect to clustering results. The exploration of underlying structure is facilitated with the help of an Icicle Plot that represents generated clusters. Additionally, a drill-down mechanism (collapsing and expanding Icicle Plot elements) provides overview and detail on demand.

In future work, we plan to address the following issues related to the visualization time-oriented data. Our current approach to resolve strong data oscillation sacrifices the natural order in time-series. A first idea to solve that problem is to omit rearrangement of data vectors, but instead extend the heuristic for color scale segmentation. The goal is to succeed in finding a segmentation such that both the segmentation is expressive and the number of crucial color switches is reduced. Secondly, since people are used to look at horizontal time axes, we plan to implement a flipping functionality such that the table can be rotated by 90 degrees counterclockwise. This would provide the option to switch between a table and a chart-like representation on demand.

It is planned to achieve these goals in cooperation with Ben and other expert biologists. A large set of biological wet lab data will be available, which provides ample opportunity to apply, evaluate, and refine our Table Lens extensions.

Figure 1 .

 1 Figure 1. Two-Tone Pseudo Coloring.

Figure 2 .

 2 Figure 2. TTPC Table Lens with different numbers of color segments (left: n = 1 to right: n = 4).

1 .

 1 Find smallest k such that 10 ≤ r * 10 k < 100, with r = max -min. Scale value range by setting max = max * 10 k and min = min * 10 k ⇒ (max = 38.7, min = 3). 2. Determine smallest i with (max + i) mod 10 = 0. Set max = max + i ⇒ (max = 40). 3. Calculate expanded value range r = max -min and determine distance between transition points d = r /n ⇒ (r = 37, d = 9.25). 4. Determine smallest j with (d + j) mod 2 = 0. Set min = max -n * (d + j) ⇒ (min = 0). 5. Scale back to original magnitude min= min * 10 -k and max= max * 10 -k ⇒ (max= 4, min= 0).

Figure 3 .

 3 Figure 3. Color scale segmentation -(a) Simplistic segmentations are difficult to interpret; (b) Our heuristic helps in generating more intuitive segmentations; (c) Users can always intervene and adjust the color scale to their particular needs.

Figure 4 .

 4 Figure 4. Hybrid clustering to reduce strong data oscillation -(a) Original arrangement of rows with many color switches; (b) Hybrid clustering reduced the number of color switches, the cluster hierarchy is represented in an Icicle Plot (fully expanded); (c) Starting from a fully collapsed Icicle Plot, users can perform information drill-down.

Figure 5 .

 5 Figure 5. Iris dataset -(a) Data sorted with respect to petal length; (b) Rows adjacent to discontinuity are expanded to get detail information.

  6(c) columns 4,6, 7, and 8).

Figure 6 .

 6 Figure 6. Countries dataset -(a) Unsorted view; (b) After clustering has been applied, a cluster of interest can be selected (denoted by blue bar); (c) Expanding that cluster reveals former east block countries.

Figure 7 .

 7 Figure 7. Climate dataset -(a) A correlation between temperature and vapor content (4th and 5th column) can be guessed just from looking at the first 365 (out of 40419) data vectors; (b) Applying clustering and collapsing the entire cluster hierarchy helps in confirming the correlation.

  Table Lens and TTPC. The basic idea is to replace the original Table Lens bars by TTPC value representations. The result is a Table Lens that shows TTPC views stacked as table columns. The integration of TTPC allows for more columns to be displayed, because TTPC supports precise reading of data values, even if the table cells are narrow. 3 We understand TTPC as a complementary technique for the Table Lens. All features of the original Table Lens are still supported:

Table 1 .

 1 Test datasets and their sizes

			health countries
	variables	10	10	10
	data items	40419	1096	110

Table 2

 2 where k is the number of data vectors. The tests we conducted indicate that the costs of SOM training and hierarchical clustering are balanced in that case. This has basically two advantages.

	neurons climate health countries
	0	27.41 23.57	51.20
	5	26.52 12.86	30.69
	10	25.74 12.71	25.87
	20	25.08 12.39	20.03
	50	23.50 11.42	-
	100	21.63 10.98	-
	500	17.19	-	-
	1000	15.54	-	-

. Pure SOM clustering quality -Lower values indicate a lower percentage of crucial color switches, and hence, higher quality; "-" denote test series that were not expected to lead to significant result variations.

Table 3 .

 3 Hybrid clustering quality -Lower values indicate a lower percentage of crucial color switches, and hence, higher quality; "-" denote test series that were not expected to lead to significant result variations.

ACKNOWLEDGMENTS

We would like to thank Prof. Takafumi Saito from Tokyo University of Agriculture and Technology for providing the original version of Figure 1.