Multivariate transfer function design

Visualization and exploration of volumetric datasets has been an active area of research for over two decades. During this period, volumetric datasets used by domain users have evolved from univariate to multivariate. The volume datasets are typically explored and classified via transfer function design and visualized using direct volume rendering. To improve classification results and to enable the exploration of multivariate volume datasets, multivariate transfer functions emerge. In this dissertation, we describe our research on multivariate transfer function design. To improve the classification of univariate volumes, various one-dimensional (1D) or two-dimensional (2D) transfer function spaces have been proposed; however, these methods work on only some datasets. We propose a novel transfer function method that provides better classifications by combining different transfer function spaces. Methods have been proposed for exploring multivariate simulations; however, these approaches are not suitable for complex real-world datasets and may be unintuitive for domain users. To this end, we propose a method based on user-selected samples in the spatial domain to make complex multivariate volume data visualization more accessible for domain users. However, this method still requires users to fine-tune transfer functions in parameter space transfer function widgets, which may not be familiar to them. We therefore propose GuideME, a novel slice-guided semiautomatic multivariate volume exploration approach. GuideME provides the user, an easy-to-use, slice-based user interface that suggests the feature boundaries and allows the user to select features via click and drag, and then an optimal transfer function is automatically generated by optimizing a response function. Throughout the exploration process, the user does not need to interact with the parameter views at all. Finally, real-world multivariate volume datasets are also usually of large size, which is larger than the GPU memory and even the main memory of standard work stations. We propose a ray-guided out-of-core, interactive volume rendering and efficient query method to support large and complex multivariate volumes on standard work stations.

MULTIVARIATE TRANSFER FUNCTION DESIGN by Liang Zhou A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Computing School of Computing The University of Utah August 2014 Copyright © Liang Zhou 2014 All Rights Reserved The Unive r si t y of Utah Graduat e School STATEMENT OF DISSERTATION APPROVAL The dissertation of Liang Zhou has been approved by the following supervisory committee members: Charles Hansen Chair 04/29/2014 Date Approved Christopher Johnson Member 04/29/2014 Date Approved Ross Whitaker Member 04/29/2014 Date Approved Valerio Pascucci Member 04/29/2014 Date Approved Claudio Silva Member 04/29/2014 Date Approved and by Ross Whitaker Chair of the Department of _____________________School of Computing and by David B. Kieda, Dean of The Graduate School. ABSTRACT Visualization and exploration of volumetric datasets has been an active area of research for over two decades. During this period, volumetric datasets used by domain users have evolved from univariate to multivariate. The volume datasets are typically explored and classified via transfer function design and visualized using direct volume rendering. To improve classification results and to enable the exploration of multivariate volume datasets, multivariate transfer functions emerge. In this dissertation, we describe our research on multivariate transfer function design. To improve the classification of univariate volumes, various one-dimensional (1D) or two-dimensional (2D) transfer function spaces have been proposed; however, these methods work on only some datasets. We propose a novel transfer function method that provides better classifications by combining different transfer function spaces. Methods have been proposed for exploring multivariate simulations; however, these approaches are not suitable for complex real-world datasets and may be unintuitive for domain users. To this end, we propose a method based on user-selected samples in the spatial domain to make complex multivariate volume data visualization more accessible for domain users. However, this method still requires users to fine-tune transfer functions in parameter space transfer function widgets, which may not be familiar to them. We therefore propose GuideME, a novel slice-guided semiautomatic multivari­ate volume exploration approach. GuideME provides the user, an easy-to-use, slice-based user interface that suggests the feature boundaries and allows the user to select features via click and drag, and then an optimal transfer function is automatically generated by optimizing a response function. Throughout the exploration process, the user does not need to interact with the parameter views at all. Finally, real-world multivariate volume datasets are also usually of large size, which is larger than the GPU memory and even the main memory of standard work stations. We propose a ray-guided out-of-core, interactive volume rendering and efficient query method to support large and complex multivariate volumes on standard work stations. iv For my wife, family and friends. CONTENTS ABSTRACT....................................................................................................................... iii LIST OF F IG U R E S ............................................................................................................ ix LIST OF TABLES .............................................................................................................. x NOTATION....................................................................................................................... xi ACKNOWLEDGEMENTS ......................................................................................... xii CHAPTERS 1...... INTRODUCTION..................................................................................................... 1 1.1 Motivation............................................................................................................ 1 1.2 Transfer Function Combinations................................................................ 2 1.3 Transfer Function Design Based on User-Selected Samples for Intuitive Multivariate Volume Exploration.............................................. 3 1.4 GuideME: Slice-guided Semiautomatic Multivariate Volume Exploration....................................................................................................... 5 1.5 Interactive Multivariate Volume Rendering and Efficient Query on Standard PCs....................................................................................................... 7 1.6 Dissertation Statement..................................................................................... 7 1.7 Dissertation Contributions............................................................................ 8 1.8 Outline................................................................................................................ 9 2. BACKGROUND..................................................................................................... 10 2.1 Transfer Function Spaces.............................................................................. 10 2.2 User Interfaces for Transfer Function Design......................................... 12 2.3 Multivariate Data Visualization................................................................... 14 2.4 Interactive Linked View Multivariate Volume Visualization System 16 2.5 Interactive Large Volume Rendering......................................................... 18 3. TRANSFER FUNCTION COMBINATIONS................................................ 20 3.1 Combining Transfer Functions ................................................................... 20 3.1.1 Formulating Transfer Function Combinations.............................. 20 3.1.2 Selecting Combinations....................................................................... 22 3.1.3 User Interface......................................................................................... 24 3.2 Specific Transfer Function Spaces U s e d .................................................. 28 3.2.1 Size Information Computation......................................................... 28 3.2.2 Occlusion Information Computation.............................................. 28 3.2.3 Statistical Properties ............................................................................ 30 3.3 Results and Discussion.......................................................................................31 3.3.1 Synthetic Dataset....................................................................................... 33 3.3.2 CT Scan of a Chest: "Artifix" ................................................................ 34 3.3.3 CT Scan of a Back P a c k ....................................................................... .... 36 3.3.4 MRI Scan of a Brain: "CerebrixCrop" ................................................37 3.3.5 Multivariate Dataset: Hurricane I s a b e l......................................... .... 39 4. TRANSFER FUNCTION DESIGN BASED ON USER-SELECTED SAMPLES FOR INTUITIVE MULTIVARIATE VOLUME EXPLORATION ..41 4.1 Method Overview.............................................................................................. ..41 4.2 Voxel Query and PCP Generation.............................................................. .... 43 4.2.1 GPU-based Voxel Query via Conditional Histogram Computation .... 43 4.2.2 Parallel Coordinate Plots Generation.................................................. 45 4.3 Transfer Function Generation from User-Selected S am p le s.................. 45 4.3.1 Sample Selection in the Multipanel View....................................... .... 46 4.3.2 Kernel Density Estimation-based Transfer Function Generation .... 46 4.3.3 Automated Gaussian Transfer Functions on Dimensionality Reduced Space .... 49 4.4 Feature Refinement in the Spatial Domain.................................................. 52 4.4.1 Screen Space Brush in the 3D V iew ................................................. .... 52 4.4.2 Screen Space Lasso in the 3D View. ................................................ .... 53 4.4.3 Refinement Brush in the Panel View................................................... 53 4.5 Rendering.............................................................................................................. 53 4.6 User Interface ....................................................................................................... 54 4.6.1 Multipanel Viewer................................................................................ .... 54 4.6.2 HDTF Editor................................................................................................ 54 4.6.3 Projection Viewer....................................................................................... 55 4.7 Use Cases .............................................................................................................. 56 4.7.1 Hurricane Isabel Simulation.............................................................. .... 56 4.7.2 3D Seismic Dataset................................................................................ .... 58 4.8 Implementation ................................................................................................ ..61 5. GUIDEME: SLICE-GUIDED MULTIVARIATE EXPLORATION OF VOLUMES................................................................................................................ .... 62 5.1 Method Overview ......................................................................................... .... 62 5.2 Guided Uncertainty-aware Lasso .............................................................. .... 64 5.2.1 Boundary Extraction................................................................................ 64 5.2.2 Boundary Confidence Image ................................................................ 66 5.2.3 Guided Uncertainty-aware Lasso......................................................... 66 5.3 Automated Feature Extraction....................................................................... ..67 5.3.1 Automated Transfer Function T u n in g ................................................67 5.3.2 3D Connected Component Extraction................................................71 5.4 Volume Rendering and Spatial Fine Tuning................................................ 72 5.5 Implementation .............................................................................................. .... 73 5.6 Examples ................................................................................................................ 73 5.6.1 Seismic Dataset ..................................................................................... .... 74 vii 5.6.2 Brain S c a n .............................................................................................. .... 75 5.6.3 Performance ..............................................................................................77 5.6.4 Discussion .............................................................................................. .... 78 6. INTERACTIVE RENDERING AND EFFICIENT QUERYING FOR LARGE MULTIVARIATE VOLUMES ON CONSUMER LEVEL PCS . 80 6.1 Multivariate Out-of-Core Volume Rendering......................................... .... 80 6.1.1 Virtual Memory Structure for Multivariate Volumes................ .... 80 6.1.2 Multivariate Transfer Functions......................................................... ..81 6.2 Efficient Multivariate Q u e ry ....................................................................... .... 82 6.2.1 Per-block Gaussian Mixture Model Computation....................... .... 82 6.2.2 Runtime Ellipse-Polygon Intersection Test .................................. .... 82 6.3 Result ....................................................................................................................... 83 7. CONCLUSIONS AND FUTURE W O RK ....................................................... .... 85 APPENDIX: RELATED PUBLICATIONS .............................................................. .... 88 REFERENCES...................................................................................................................... 89 viii LIST OF FIGURES 3.1 The separable 3D transfer function editor.................................................... .... 26 3.2 Comparisons of transfer function combinations for a synthetic dataset. 32 3.3 Transfer function combinations for a CT chest scan...................................... 35 3.4 Transfer function combinations for a back pack CT scan........................ .... 36 3.5 Transfer function combinations for an MRI brain scan data....................... 38 3.6 Classifying a multivariate dataset using the combined transfer function space....................................................................................................... .... 40 4.1 The user interface and the work flow of the system implementing our proposed method............................................................................................. 42 4.2 Generating a PCP from a joint histogram......................................................... 45 4.3 Sample selection in the multipanel view. 47 4.4 Transfer function generated from user-selected sample points............. .... 48 4.5 Feature refinement tools.................................................................................... .... 52 4.6 The high-dimensional transfer function editor........................................... .... 56 4.7 Results of a hurricane simulation dataset 57 4.8 Extracted geological features in a seismic dataset......................................... 59 4.9 User selected samples and classified features of a seismic dataset . . . 60 5.1 The workflow of GuideME............................................................................... .... 63 5.2 The inspection window and the boundary confidence image............... .... 65 5.3 Illustrations of images involved in the automated transfer function tuning process........................................................................................................... 69 5.4 Illustrations of the transfer function modification process and effects of the process............................................................................................................. 72 5.5 Results of the New Zealand seismic dataset............................................... .... 74 5.6 Visualization results of multimodal brain scans using GuideME 77 6.1 Transfer functions for the channel system of a large seismic dataset. . . 81 6.2 Our proposed method allows interactive visualization and efficient query of large multivariate seismic datasets on consumer level PCs. 84 LIST OF TABLES 3.1 The correlation coefficients and the entropies of the properties computed from Artifix CT chest scan dataset............................................. 25 5.1 Quantitative comparisons and timing results for the features extracted in the example datasets................................................................... 76 NOTATION One-dimensional (1D) Two-dimensional (2D) Three-dimensional (3D) Transfer Function (TF) Graphics Processing Unit (GPU) Open Graphics Library (OpenGL) OpenGL Shading Language (GLSL) Parallel Coordinate Plot (PCP) Personal Computer (PC) Computed Tomography (CT) Magnetic Resonance Imaging (MRI) ACKNOWLEDGEMENTS It was the interest of computer graphics that brought me to University of Utah in August 2008. I took the interactive computer graphics class taught by Professor Charles (Chuck) Hansen, my advisor-to-be. The class gave me an opportunity to become a research assistant in Chuck's research group, which tackles challenging scientific visualization problems. It actually took me a while to fully understand and appreciate the importance and beauty of visualization and to shift my interest from graphics to visualization. I would like to thank Chuck for being so nice as an advisor and for ensuring the funding for my study. Moreover, I wish to thank him for his style of advising: encouraging students to do research projects independently from the beginning to the end. He never tells you what to do next, but expects you to do it by yourself first and then helps to make it better with a full heart. I have to admit that the learning curve is steep but the gains from this process are tremendous and life-long. A collaboration with ExxonMobil on multivariate volume visualization and analysis was built in the early days in my Ph.D. study and lasted throughout the process and is still active. It was my great pleasure to meet and know Mark Dobin and Guo-shi Li who enabled my research to help real users in the petroleum industry. They also made it possible for me to do an internship in the Upstream Research Company from which several of my research topics stemmed. I would like to thank Mark for his passion, creativity, generosity and helpfulness. Mark, as a geophysicist, always thinks out of the box and gives suggestions for my work. I would like to thank Guo-shi for being such a great mentor. It has been my privilege to work in the SCI Institute for the past five years. SCI has allowed me to meet and interact with peer students in different research areas, and to learn cutting edge techniques from the talks given by distinguished researchers. SCI has given me the unique opportunity to work in such a great environment where I can enjoy the beauty of the slopes of the city and aromatic coffee. I wish to thank Chris Johnson for making it all possible and being a member of my committee. I would also like to thank Ross Whitaker, Valerio Pascucci and Claudio Silva for kindly being my committee members. It was a joint effort from all of you that introduced me to the world of visualization and helped me to prepare for the research challenges. I also wish to thank the entire administration team of SCI. I wish to thank my colleagues for their helps. It would have been impossible to finish this dissertation without the help from members of Chuck's research group. I thank Mathias Schott who helped me write my first manuscript and answered all problems I have on coding and graphics hardware. I would also like to thank Pascal Grosset, Yong Wan, Carson Brownlee, Mark Kim, Guoning Chen and Aaron Knoll for their help in academic and life throughout my study. I remember the inspiring and fruitful discussions with Bo Wang and Wei Liu on many topics related to image processing and machine learning. I would like to thank them and Xianzong Xie, Zhan Wang and Dafang Wang for being such good friends who have helped me out more times than I can count. My family has always been my greatest support and the biggest reason why I have been working hard. It is my parents who have made it possible for me to pursue my goals in life without any hesitation, and it is also their support that helped me to overcome all the difficulties in this place an ocean away from home. I wish to thank my wife Peng Lin who is such a wonderful companion of mine and helps me not only keep everything in good shape in everyday life but also to find the next goal in my career. Also, I am fortunate to have such caring and supportive parents-in-law and all other family members. Finally, my research was made possible by the following grants: Award No. KUS-C1-016-04, made by King Abdullah University of Science and Technology (KAUST), DOE SciDAC:VACET, NSF OCI-0906379, NIH-1R01GM098151-01. DOE NNSA Award DE-NA0000740, DOE SciDAC Institute of Scalable Data Manage­ment Analysis and Visualization DOE DE-SC0007446, NSF ACI-1339881, NSF IIS-1162013. xiii CHAPTER 1 INTRODUCTION 1.1 Motivation Volumetric data are important to many application domains, and are usually visualized using direct volume rendering. In recent years, the volume datasets utilized by domain users have evolved from univariate to multivariate. Multi­variate volume datasets have become increasingly important and popular. Me­teorologists simulate atmosphere movements, e.g., hurricanes, which involves the interaction of dozen of physical parameters to study their causes and try to forecast them. Combustion simulations involving multiple chemical elements and physical measurements can help engineers to improve the efficiency of an engine. Geophysicists compute several attributes derived from the amplitude of three-dimensional (3D) seismic survey data to facilitate the exploration of oil and gas [12]. Physicians determine the location of lesions in the brain using MRI scans that contain multiple channels. Transfer function design is the major means for exploring volume data and extracting features in direct volume rendering. Transfer functions of 1D or 2D are commonly used for the exploration of uni­variate volume datasets, and multivariate transfer functions are used to explore multivariate volume datasets. However, extracting clear features of interest for univariate volumes using 1D or 2D transfer functions is usually difficult due to the limited classification ability of these low-dimensional transfer functions. On the other hand, although methods and systems have been proposed for multivariate volume exploration and visualization, domain users urgently need an intuitive and flexible multivariate visualization tool that is able to extract meaningful features in complicated real-world datasets. Moreover, due to the boost of accuracies of acquisition devices, storage and computational capacities of computers, multivariate volume datasets are generated with ever increasing size. 2 Giga-scale or even tera-scale multivariate volume datasets have become common, and fast volume rendering and query techniques are vital for the visualization of these datasets. In this dissertation, we present our proposed approaches to the multivariate transfer function design problem. This chapter first introduces transfer function combinations, in which we propose to combine the best features of existing transfer function spaces to create a transfer function space that provides better classification. It then introduces a transfer function design approach based on user-selected samples in the spatial domain to make multivariate volumetric data visualization more accessible for domain users. Next, it introduces GuideME: an automated technique for designing optimal multivariate transfer functions with a simple and easy-to-use slice-based user interface for highly complex volumes. Finally, it introduces our work on scalable out-of-core methods for interactive rendering and efficient querying for large multivariate seismic volumes on consumer level PCs. 1.2 Transfer Function Combinations Direct volume rendering is an active area of research. Mapping of data values to optical properties, known as classification, remains a challenging problem. Transfer functions are most commonly used for classification in volume rendering, yet finding good transfer functions remains a difficult problem. For material boundaries, it has been shown that 2D transfer functions provide greater speci-ficity1 than 1D transfer functions [54,60]. In many datasets, separate features may share the same scalar value and gradient magnitudes and as such scalar value, gradient magnitude tuples are not sufficient for separating such features. Recently, many new 2D transfer function spaces have been proposed to improve the classification from different metrics. The size-based transfer function is a transfer function space [14] built upon blob detection techniques using scale space theories to classify objects based on their sizes. The occlusion spectrum is another 2D transfer function space [15] that takes into consideration ambient occlusion 1 We use the disambiguation definition of specificity rather than the statistical definition, which means the proportion of negatives in a binary classification test that are correctly identified. 3 within the volume for discriminating between features of similar scalar values. It is also possible to compute statistical measurements such as mean value and standard deviation in a local region around a voxel [38] to form a 2D transfer function space. All these methods are effective on some datasets. Other datasets, however, may contain materials that have similar statistical properties but occlude each other, or have materials that share similar statistical properties and occlusion measurements but differ in size. We propose, therefore, to combine the best features of these transfer functions to create a transfer function space that provides better specificity. Our contributions in this work are threefold: 1) combining 2D transfer function space with 1D transfer function spaces with a basic approach for selecting combinations, 2) a user interface supports transfer function design in the combined transfer function space, 3) experiments and detailed discussions of different transfer function combinations and original 2D transfer functions on various datasets. We experimented with combinations of these transfer function spaces and discuss a basic approach for selecting combinations that improve classification and show that this combined transfer function space provides better classification than 2D gradient magnitude transfer functions, 2D statistical transfer functions, 2D occlusion-based transfer functions or 2D size-based transfer functions. 1.3 Transfer Function Design Based on User-Selected Samples for Intuitive Multivariate Volume Exploration Multivariate dataset visualization has been an active research area for the past decade and remains a challenging topic. A linked-view visualization system that enables the users to explore the datasets in both the transfer function domain and the spatial domain may boost their understanding of the data. In recent years, visualization researchers have been studying this topic and some solutions have been proposed [21, 1, 5, 31]. These linked-view systems provide users the ability to explore the dataset with closely linked scientific visualization views, e.g., volume rendering or isosurface rendering, and information visualization views, e.g., scatter plots, parallel coordinate plots (PCP) or dimensionality reduction 4 views. Typically, the user explores and extracts features of interest by interactively designing transfer functions (TFs) in the value domain over the information visualization contexts and examining classified results in the spatial domain from the scientific visualization view. Successful examples using these systems are clearly shown for simulation datasets. However, extracting meaningful features in real-world measurement datasets, e.g., multivariate 3D seismic survey, via these systems is not trivial. Features inside the seismic dataset have to be recognized in the spatial domain by a geology expert, and the features have complicated combinations of attribute values and subtle differences from their surroundings. Therefore, it is too laborious to extract features by iterating between TF design in the value domain and getting feedback from the results rendered in the spatial do­main, especially when the dimensionality is high. Our geophysicists collaborators have found extracting features with only the value domain TF widgets, e.g., on a PCP to be cumbersome, and specifically asked for more automated methods. We propose a TF design approach based on user-selected samples from the spatial domain represented as slices for more intuitive exploration of multivariate volume datasets. Specifically, the user starts the visualization by probing fea­tures of interest in a panel view, which simultaneously displays associated data attributes in slices. Then, the data values of these features can be instantly and conveniently queried by drawing lassos around the features or, more easily, by applying "magic wand" strokes. High-dimensional transfer functions (HDTFs) can then be automatically and robustly generated from the queried data samples via the kernel density estimation (KDE) [87] method. The TFs are represented by parallel coordinate plots (PCPs) and can be interactively modified in an HDTF editor. Automatically generated Gaussian TFs in dimensionality reduced 2D view can also be utilized to extract features. The extracted features are rendered in the volume rendering view using directional occlusion shading to overcome artifacts from Phong shading in the multivariate case. To further refine features, which share similar data value ranges, direct volume selection tools on the volume rendering view or the panel view can be applied. The contributions of this work are the following: First, we propose a transfer 5 function design method for multivariate volume visualization based on user-selected samples, specifically an HDTF generation method based on KDE and a Gaussian mixture model based 2D Gaussian TF generation method. Second, we also propose an interactive multivariate volume visualization system based on the proposed method that has been implemented to allow domain users to extract refined features in very complicated multivariate volume datasets more intuitively. 1.4 GuideME: Slice-guided Semiautomatic Multivariate Volume Exploration The state-of-the-art method for exploring multivariate volumes is user inter­action with multiple linked view systems. This method requires the user to explore the volume using parameter views, e.g., parallel coordinate plots (PCP) or histograms, in a trial-and-error manner [1, 5, 31]. Although these systems have been successful in simulation datasets where the user understands the "recipe" of the parameter space, i.e., knows what combinations of value ranges of attributes may result in interesting features, it is difficult for the user to explore complex measured datasets, e.g., seismic datasets. To this end, research efforts have addressed the exploration of complex datasets such as seismic data [103, 43, 44]. In [103], the approach allows domain users to apply their expertise to the finding of features by directly selecting a region of interest in a multipanel slice view. However, these methods either work only on univariate seismic data for a certain type of features [43, 44] or still require transfer function tuning with a PCP-based or a histogram-based editor [103], which can be unintuitive and time consuming for domain users. Moreover, switching between multiple views may be somewhat distracting. Seismic datasets are an important tool for the petroleum industry which is the driving application of our method. Geophysicists interpret features that indicate potential oil and gas reservoirs, including channels and salt domes, on 3D seismic data slices. To interpret the seismic data, they first identify and locate geological fea­tures on slices from the 3D seismic data through examination and selection. With the advancement of multivariate 3D seismic [12] interpretation, attributes derived 6 from the seismic amplitude are used to aid in the extraction of relevant features. Interpretation is done mainly by free-hand drawing on slices and refinement of the features through multiattribute transfer functions [103, 90]. Users typically perform the following tasks during interpretation: selecting features by drawing on 2D slices, refining those features through transfer function manipulation, and examining results in 3D renderings. Other domains use similar tasks and we demonstrate the generality of our approach with multimodal MR brain scans. We therefore propose GuideME, a novel method for multivariate volume exploration that strives to provide the user with a very simple and intuitive exploration process for highly complex datasets. Instead of multiple linked views, our method has only one slice view coupled with a focus overlay and a volume rendering view, and the tedious trial-and-error interactions are largely replaced by a guided uncertainty-aware lasso and automated feature extraction. The user starts the exploration by browsing through the slices and detects a feature of interest using his/her domain knowledge. A focus window that allows the user to inspect other attributes that can be placed over the feature of interest. Through the inspection, the user determines attributes that best describe the boundaries of the feature. A boundary confidence image is then blended with the slices, and the user can easily select the region with a guided uncertainty-aware lasso that automatically snaps to the detected feature boundaries. The selected region is then used as input for the automated feature extraction. Eventually, the feature is volume rendered and may be optionally edited directly in the 3D view. Using a highly complex real-world seismic dataset and multimodal MR brain scans, we show our approach is efficient, and is able to create results comparable to those given by previous method and ground-truth segmentations. In this work, we make the following contributions: • A novel slice-guided semiautomatic multivariate volume exploration work­flow. The user is freed from unfamiliar parameter space views and tedious trial-and-error transfer function tunings. • A guided uncertainty-aware lasso for region selection, based on edge detec­tion and Dijkstra shortest path algorithms. 7 • A technique to automatically fine-tune a multivariate transfer function given the lasso, based on the optimization of a response function. 1.5 Interactive Multivariate Volume Rendering and Efficient Query on Standard PCs Due to the advance in 3D seismic imaging techniques, resolution of 3D seis­mic datasets used in the petroleum industry is usually of giga-bytes. Multiple attributes derived from the original seismic amplitude dataset have been used for aiding the interpretation of the seismic survey. With these additional attributes, however, the size of the entire dataset may become far larger than the capacity of the GPU memory and even the main memory of a typical workstation. Recently, GPU-based multiresolution out-of-core volume rendering systems have been pro­posed. The Gigavoxel approach by Crassin et al. [16] and CERA-TVR by Engel [23] divide the dataset into multiresolution bricks and utilize an octree structure and ray guided paging system to efficiently render large sparse volume datasets. Had-wiger et al. [35] propose a rendering system that uses a virtual memory system and 2D mipmapping to support dense and noisy petascale microscopy scans. However, none of these methods are able to handle large multivariate volume datasets. In this work, we extend the approach by Hadwiger et al. to support interactive rendering of large multivariate seismic datasets on a consumer-level PC. On the other hand, data value querying raises another challenging issue for multivariate datasets especially when the data are very dense as hierarchical acceleration techniques, e.g., octrees may not be beneficial. We propose an efficient data querying technique based on precomputed per-block Gaussian mixture models and run-time ellipse-polygon intersection detection. An interactive exploration system has been built to allow the user to visualize the multivariate volumes as well as to edit multivariate transfer functions with the query feedback on parallel coordinate plots. 1.6 Dissertation Statement Multivariate transfer function design is vital for the exploration and visual­ization of univariate and multivariate volumetric datasets, which are important 8 to many real-world applications. The design of multivariate transfer functions requires both existing and novel techniques in volume rendering, user interaction, image processing, data analysis and optimization. 1.7 Dissertation Contributions The main contributions of this dissertation are: • More intuitive multivariate volume exploration and visualization work­flows designed for domain users. Our proposed methods strive to maintain the power and flexibility of existing workflows, while overcoming the short­comings of existing approaches that are designed from the perspective of visualization experts. We provide domain users with new work flows that are designed as a result of close collaborations. The resulting workflows give domain users a central role and try to reduce unfamiliar widgets and views that may hassle users. For example, we allow the user to select samples directly on data slices to generate the initial transfer function, and then fine-tune the transfer function using parameter space widgets. A further developed workflow can be seen in GuideME, in which the parameter space is completely hidden from the user and the manual transfer function tuning process is replaced by an optimization approach. Through real-world examples, domain users find the new workflows more intuitive and efficient. • Improved classification ability of transfer functions demonstrated on com­plex real-world datasets. We propose a transfer function combination method that generates multidimensional transfer function spaces by ranking and selecting the most helpful 1D or 2D transfer function spaces computed from a univariate volume. Moreover, transfer functions do not contain spatial information. Easy-to-use spatial fine tuning tools are therefore provided. Diffusion-based region growing tools and lasso tools are provided for effi­cient volume editing on transfer function classified results on either 2D or 3D image spaces. Complex real-world datasets, including CT lung scan, MR brain scans and multivariate seismic datasets, have been successfully used to demonstrate the improved classification ability of transfer functions. 9 • Improved user interactions in the visualization process. The goal for easier visualization experience for domain users is considered in all our proposed work. An easy-to-understand user interface that combines 2D and 1D transfer function widgets is built for the transfer function combination work. The user interface of the multivariate volume visualization system that is based on user-selected samples provides synchronized panel view, click-and-drag transfer function editors and hand drawing-based region selection and spatial fine tuning tools on 2D and 3D views. In GuideME, the user is given suggestions of feature boundaries and uncertainty information, and the selection of region is aided by lassos that automatically snap to feature boundaries. 1.8 Outline The background and important related work are explained in Chapter 2. Chap­ter 3 details the work of transfer function combinations. The work of multivariate transfer function design based on user-selected samples in the spatial domain is presented in Chapter 4. Chapter 5 presents the method of GuideME, slice-guided multivariate exploration of volumes. In Chapter 6, we present the work of interactive rendering and efficient query of large multivariate seismic volumes on consumer-level PCs. Finally, conclusions and future work are given in Chapter 7. CHAPTER 2 BACKGROUND 2.1 Transfer Function Spaces Volume datasets can be explored using transfer functions. The most frequently used transfer function for volume rendering is a 1D transfer function that uses scalar values for classification. Realizing the poor classification ability of that transfer function space, Levoy [60] and Kindlmann et al. [54] use the gradient magnitude of the volume as another property for better classification. Kniss et al. [56] advocate and implement multidimensional transfer functions widgets, making the 2D transfer function a standard method in modern volume renderers. By far, the 1D and 2D transfer functions are the most popular and practical techniques for classification in volume rendering; however, great efforts have been made to define new transfer function spaces to improve the classification ability. Due to noise and partial volume effects, selecting a boundary in the arches of a gradient magnitude-based transfer function is not easy and sometimes even impossible. To resolve this problem, Lundstrom et al. [64] employ the local histograms to better discern tissues in medical datasets and propose a 2D transfer function space that uses competitive classification certainty measure in addition to scalar values. Sereda et al. [86] use a 2D LH histogram-based transfer function for easier boundary identification and selection and further use this boundary information for a region growing segmentation schema. Higher degree derivatives of the original scalar volume can also be used as transfer function spaces. Kindlmann et al. [55] use curvature as a second dimen­sion of their transfer function domain. The curvature-based transfer functions allow nonphotorealistic renderings that highlight the contours of the volume. 11 The theory of scale-space, developed originally by the computer vision and image processing communities, can be used to classify objects based on their sizes. A commonly used scale-space representation is the linear Gaussian scale-space, which is essentially a convolution of a volume with differently sized Gaussian filters. Lum et al. [63] combine it with an image pyramid representation of different scales to improve classification. Correa and Ma [14] create a continuous scale-space for the volume and use anisotropic diffusion to detect "blobs" in the volume. The size of these defines an additional metric of the volume, which is then used to create size-based transfer functions. Shape is another important aspect to classify an object. Sato et al. [83] use eigen­value analysis on 3D local intensity structures to classify tissues in medical datasets with 2D transfer function spaces created using shape measurements: sheet, line or blob along with the scalar value. Prassni et al. [78] propose shape-based transfer functions by computing shape descriptors over presegmented volume to provide a manageable set of shape classified volumetric features with an intuitive optical properties assignment interface. In many cases, different features occlude with each other but share similar scalar values, e.g., the skin and the gray matter in an MR brain scan. Correa and Ma [15] use the occlusion of a voxel as an additional dimension of the transfer function domain to classify features of similar scalar value, but different local neighborhoods. Volumes can also be classified based on their statistical metrics, such as mean value or standard deviation of voxels in a certain neighborhood. Caban et al. [10] compute local statistical metrics and use their linear combinations to classify fine structures. Patel et al. [74] use a dynamically changing neighborhood to compute mean value and variance for voxels, thus defining a transfer function domain. A user interface then allows a selection of features based on the mean value, variance and radius of the neighborhood. Haidacher et al. [38] further extend this approach by selecting the radius semiautomatically via an adaptive sample selection technique. Transforming the volume data into frequency domain is another idea for 12 generating transfer function spaces. Vucini et al. [94] utilize GPU-based fast Fourier transformation to support interactive frequency-based transfer function design that enhances conventional volume visualization. 2.2 User Interfaces for Transfer Function Design A 1D transfer function that uses scalar values of the volume or a 2D transfer function that has the gradient magnitude of the volume as a second property for better classification [54] are most frequently used. As in most volume visualization systems available nowadays, e.g., Voreen [93], VisIt [62] and ImageVis3D [48], the transfer functions can be interactively defined by 1D transfer function widgets or 2D transfer function widgets proposed by Kniss et al. [56]. However, to design a good transfer function, the user has to manipulate the transfer function widgets in the transfer function space and check the result in the volume rendered image, which is laborious and time consuming. To address this issue, researchers have proposed to automate the transfer function design process. Many researchers focus on the automation of user interactions on the transfer function space. Maciejewski et al. [65] utilize KDE to structure the data value space to generate initial transfer functions. Wang et al. [96] initialize transfer functions by modeling the data value space using the Gaussian mixture model and render the extracted volume with preintegrated volume rendering. Most recently, Ip et al. [50] propose a hierarchical visual segmentation method using normalized-cut on the intensity-gradient magnitude 2D transfer function space to assist the volume exploration process. Although these automated methods applied on the transfer function space significantly reduce the time a user spends in the volume exploration process [50] compared to the most commonly used transfer function widgets, interacting with the transfer function space is not intuitive. Therefore, the potential of transfer function design on more intuitive spaces has also been studied. One strategy is to provide the user with a gallery of predefined transfer functions, and then the user can easily design customized transfer functions by picking features of interest from the gallery and refining them. Marks et al. [67] propose Design Gallery as a general user interface for 13 computer graphics applications. As for volume rendering, the design gallery automatically defines a set of random 1D transfer functions and generates the resulting thumbnails by the dispersion heuristic. The resulting thumbnails are arranged using multidimensional scaling algorithm. The user can then select a preferred thumbnail and fine tune its associated transfer function. Similarly, a spreadsheet-like interface is proposed by Jankun-Kelly and Ma [51] where the user can explore a range of parameter combinations at the same time and compare the results side-by-side. Tory et al. [89] propose a parallel coordinates style interface that provides an overview of rendering options and transfer function settings. The user can easily explore different parameters and backtrack the visualization history using the proposed interface. Tzeng and Ma [92] propose to use ISODATA (Iterative Self-Organizing Data Analysis) clustering on a small subset of all data voxels to classify the volume, and then these classified features are displayed in a cluster-space user interface. With this user interface, the user is able to design transfer functions without knowing the transfer function space, simply by picking features from the gallery of preclassified features and refining them using the clustering operation buttons and rendering property adjustment widgets. Interacting with the spatial domain views, e.g., sketching on the volume rendering view, is another intuitive option for the user. Several methods have been proposed to enable the user to design transfer functions by sketching on the spatial domain. Tzeng et al. [90, 91] propose systems that allow users to sketch the volume slices to assign color and transparency, and then high-dimensional transfer functions are generated using artificial neural network. Based on user sketches on the rendered images, Wu and Qu [99] fuse multiple features in distinct rendering results into a comprehensive visualization. Ropinski et al. [80] propose a sketch-based user interface for 1D transfer function design where the user draws strokes to define foreground and background to extract layers in the volume, and then the transfer functions are refined by adjusting the color and opacity of each layer. A more convenient sketch based volume exploration system is proposed by Guo et al. [30] where a full set of tools have been developed to enable direct manipulation of color, transparency, contrast and other optical properties on the 14 volume rendering view by means of user drawn strokes. In its essence, the system intelligently defines 1D transfer functions based on user sketches. Other researchers have proposed methods to design transfer functions by semantics. Salama et al. [82] propose a framework that allows visualization experts to design high-level transfer function user interface with semantic information. Given specific features of interest by the domain user, e.g., bones, skin and blood vessels, the visualization expert creates transfer function models described by sets of control points from many data instances of the same type. Each transfer function model can be written as a data point in a high-dimensional space, and the semantic parameters can be generated using principal component analysis (PCA) on these high-dimensional data points. A simple semantics editor user interface can then be created from the semantic parameters that modify the control points of the low level 2D transfer function widgets. The user simply needs to adjust a set of sliders for each semantic parameter. 2.3 Multivariate Data Visualization Visualizing and understanding multidimensional datasets has been an active research topic in information visualization. The scatter plot matrix is a straight­forward yet scalable way that utilizes pairwise scatterplots as matrix elements to represent multidimensional datasets. The scatter plot matrix trades the resolution of each scatter plot to display more plots. Due to the large amount of plots shown in a scatter plot matrix, exploration of the multivariate space becomes cumbersome and time consuming. Elmqvist et al. [22] therefore propose a method for navigating through the plots. Other researchers have studied how to choose the display order of the scatter plots [85]. Parallel coordinate plots (PCPs) [49] is a popular multivariate visualization technique that overcomes the two-variable limit of the scatter plot. PCP arranges individual variable axes parallel to each other and represents individual samples as a polyline passing through all axes. With increasing samples, the PCP will become more cluttered and the rendering cost will be prohibitively expensive. A large amount of research efforts have been done to address this over-plotting 15 issue. Fua et al. [26] propose to cluster the data values and render only the representing polylines of each cluster. Novotny and Hauser [70] separate the number of polylines to render from the number of data samples by generating the PCP from 2D histograms of adjacent variable pairs, and outliers are also identified through histogram analysis. Zhou et al. [102] perform visual clustering on the PCPs by drawing curved edges instead of polylines and optimizing the arrangement of these curved edges. Heinrich and Weiskopf [41] show how PCPs can be made continuous, which gives a smooth and uncluttered representation. The over-plotting issue can also be resolved by edge bundling, and researchers have proposed a variety of techniques: e.g., geometry-based edge bundling [17], hierarchical edge bundling [45] and force-directed edge bundling [46]. McDonnell and Mueller [68] address the over-plotting issue using illustrative rendering which applies opacity and shading effects, silhouettes emphasizing, shadows and halos to edge bundled PCPs. Dimensional reduction and projection are other techniques for multidimen­sional data visualization. These techniques provide a similarity-based overview for multidimensional data. Numerous research efforts have been focused on this topic, and popular methods include principal component analysis (PCA) [53], mul­tidimensional scaling (MDS) and Isomap [88]. To reduce computation complexity, techniques that apply classical dimensional reduction and projection methods to only a small subset of representative samples and project the remaining samples via interpolation have been proposed. These techniques include Landmarks MDS [19] and Pivot MDS [8]. Alternatively, Faloutsos and Lin propose the Fastmap algorithm [24], which has exactly O(n) complexity. Fastmap utilizes dissimilarities between each sample and two pivot elements per coordinate axis to make its distance computation O(n). More recently, Paulovich et al. [75] propose the part-linear multidimensional projection (PLMP) method, which requires only distance information between pairs of representative samples, and therefore is faster than previous methods for large datasets. Moreover, with a representative sample positioning strategy, PLMP is able to conduct dimensional projection for out-of-core datasets. 16 Parallel coordinate plots capture individual dimensional information well but suffer from the clutter problem and such plots require expertise to interpret, because data points are transformed into polylines and the polylines occlude with each other. It is hard and sometimes even impossible to check the data correlation between a single pair of attributes, let alone for multiple attributes on the PCPs. On the other hand, pairwise scatter plots provide a clear correlation between a pair of attributes. A 2D scatter plot of the dimensional reduced or projected space of the high-dimensional data provides an easy-to-understand overview of the high-dimensional space at the cost of losing individual dimensional information. It is nontrivial if not impossible to use only one of these techniques to provide the user with proficient insights to a multidimensional dataset. Unfortunately, providing several linked views: one for PCPs, one for pairwise scatter plots and yet another for dimensional reduced/projected scatter plot would cause a context jump for the user. Researchers therefore have proposed to take the advantages of these techniques and combine them in a unified plot. Yuan et al. [100] propose SPPS (scattering points in parallel coordinates), which draws pairwise scatter plots between each pair of PCP axes or adopts a DIMDS (dimensional incremental multidimensional scaling) scatter plot between a selected pair of PCP axes. They convert parallel coordinates segments into point plots and draw the PCPs as curves that pass through their associated points. To provide a seemless integration, a uniform brushing tool that allows linked brushing on either the PCPs or the scatter plots is also proposed. 2.4 Interactive Linked View Multivariate Volume Visualization System Multivariate volume datasets can be explored using linked view systems that have shown to be useful for multivariate simulation data exploration. Early studies utilize multiple linked scatter plots as the data value view and the user brushes regions of interest on these plots to design transfer functions. The SimVis system [21, 77] allows the user to interact with several 2D scatter plot views using linked brushes to select features of interest in particle simulations rendered as polygons and particles. 17 As parallel coordinate plots have become a widely accepted method for mul­tidimensional data visualization, researchers propose to build data value view together with transfer function widgets based on parallel coordinate plots. Akiba and Ma [1] propose a tri-space exploration technique involving parallel coordi­nate plots together with time histograms to help the design of high-dimensional transfer functions for time-varying multivariate volume datasets. Blaas et al. [5] extend parallel coordinate plots for interactive exploration of large multitime point datasets rendered as isosurfaces. Rubel et al. [81] build a cluster-based multivariate visualization system centered on the histogram-based parallel coordinates for very large multivariate time-varying particle simulation. Given a user-selected multivariate value range on any attribute axis on the parallel coordinates, a fast multivariate query using bitmap indexing is conducted. The query result is then represented as parallel coordinate plots using the method proposed by Novotny and Hauser [70]. Finally, the particles that satisfy the query are rendered in the spatial view using particle systems. A parallel coordinate plot provides a good context for the definition of each attribute value in high-dimensional transfer functions. However, due to data point to polyline transformation and the occlusion issue by its nature, it is hard to observe high-dimensional features and check correlations between attributes in a parallel coordinate plot as stated in Section 2.3. Researchers therefore resort to the dimensional reduction and projection techniques in conjunction with parallel coordinate plots to provide the user with more insight. Zhao and Kaufman [101] combine multidimensional reduction and transfer function design using parallel coordinates but their system is able to handle only very small datasets. Guo et al. [31] propose an interactive HDTF design framework using both continuous PCPs and MDS technique accelerated by employing an octree structure. Guo et al. [32] develop parallel algorithms for multivariate volume rendering, continuous PCPs computation and MDS computation to make their work [31] scalable. However, we have observed two limitations in the above systems: 1) the user has to explore the data via interactions on the transfer function view, which may be unintuitive for domain users and moreover makes exploration for real- 18 world datasets difficult, and 2) the visualization is merely produced with transfer functions and it is difficult to achieve a more refined result. 2.5 Interactive Large Volume Rendering Ray casting and slice-based volume rendering are the two methods used for direct volume rendering. Thanks to the computational power of GPUs, volume rendering has become interactive. However, with the ever increasing size of volume datasets, interactive direct volume rendering of large volume datasets that cannot fit into the GPU memory, i.e., out-of-core volume rendering, is still a challenging topic. Building hierarchical structure for a volumetric data is a common way to enable and accelerate out-of-core volume rendering. In the early work of hardware assisted volume rendering, LaMar et al. [59] and Weiler et al. [97] propose to use hierarchical bricking schemes. Boada et al. [6] build a mipmap-like structure based on an octree, and then they choose a cut through the tree and use the mipmap data of the leaves during rendering. Guethe and Strasser [33] use hierarchical wavelet representation and screen-space error estimation for level of detail selection. The ImageVis3D system [48] uses a kD tree to subdivide data and each brick in the tree is rendered in one rendering pass. Gobbetti et al. [28] determine the visibility of octree nodes using the corresponding partial octree on the CPU, which is then downloaded to the GPU. All work above requires a CPU-based traversal of an explicit hierarchical structure, which can be very expensive. More recently, due to the improvement of the GPU, GPU-based ray-guided volume rendering frameworks enable efficient rendering of gigascale and even petascale volume data on a single consumer level GPU. Crassin et al. [16] propose a GPU-based ray-guided octree volume rendering framework, called "gigavoxel," which uses ray casting information to directly guide the data streaming. The "gigavoxel" framework is efficient as the ray casting information naturally determines voxel visibility and view frustum. However, "gigavoxel" is intended for entertainment applications, which usually result in sparse octrees, and moreover, the kD restart octree traversal scheme requires a full path traversal from the root of the tree for each voxel, which can 19 be costly. Hadwiger et al. [35] present an interactive volume renderer that scales to petascale which uses a visualization-driven virtual memory approach. Similar to "gigavoxel," ray casting is also utilized to detect visible data, however, [35] avoids the potentially costly kD restart octree traversal. Unlike "gigavoxel,: [35] it is designed and optimized for dense anisotropic microscopy data. In contrast to all previous work, the virtual memory approach requires no precomputation of a multiresolution hierarchy. Instead, it constructs volume data from the 2D raw images based on a visualization-driven data construction scheme via on-the-fly stitching and re-sampling. CHAPTER 3 TRANSFER FUNCTION COMBINATIONS Direct volume rendering has been an active area of research for over two decades. Transfer function design remains a difficult task since current methods, such as traditional 1D and 2D transfer functions, are not always effective for all datasets. Various 1D or 2D transfer function spaces have been proposed to improve classification exploiting different aspects, such as using the gradient magnitude for boundary location and statistical, occlusion or size metrics. In this chapter, we present a novel transfer function method that can provide more specificity for data classification by combining different transfer function spaces. In this work, a 2D transfer function can be combined with 1D transfer functions and improves the classification. Specifically, we use the traditional 2D scalar/gradient magnitude, 2D statistical, and 2D occlusion spectrum transfer functions and combine these with occlusion and/or size-based transfer functions to provide better specificity. We demonstrate the usefulness of the new method by comparing it to the following previous techniques: 2D gradient magnitude, 2D occlusion spectrum, 2D statistical transfer functions and 2D-size based transfer functions. 3.1 Combining Transfer Functions 3.1.1 Formulating Transfer Function Combinations In practice, using just one or two metrics during volume classification makes it difficult to robustly classify and separate features in complex volumes. Using more properties in the transfer function space often can better describe features in the volume; however, user interaction becomes more difficult or even impossible when the number of properties, and thus the dimensionality, of the transfer function space increases. Gaussian transfer functions have been proposed by Kniss et al. [57] to provide analytical multidimensional transfer functions of arbitrary 21 dimensionality. Also a procedural high-dimensional transfer function model is proposed in [36]. However, in both works, how to provide an effective user interface remains unclear. The proposed transfer function combination sequentially applies two transfer functions, a 2D and a 1D one, to all voxels vx, where X is the 3D position, of the input dataset V that has l properties Y1rY2,... Yf. {C1 , C2, C3, ... Q } = TFw( V)where Cj := {vx|TFwj(Yp(X), Yq(X)) > 0,p,q e [1,l]} (3.1) {Wi c Cj} = TFr(Yr(X))where X e Cj(X), j e [1,k],r e { 1 ,...,l}/{p, q} (3.2) In Equation 3.1, a number of 2D transfer function widgets, k, are first applied to the volume, resulting in sets of classified voxels C1 , C2, C3 ... Ck, respectively. Then one from a set of r, which is typically 1 or 2, 1D transfer functions is applied to the classified region Cj, yielding the final classified volume region Wi. Each 2D transfer function widget has one associated 1D transfer function. Kniss et al. [57] clearly show a 2D example that separating high dimensional transfer functions into lower-dimensional transfer functions using multiplication can lead to misclassification, which gets worse when the dimensionality is ex­tended into 3D. Our proposed method, however, does not suffer from such issues as each 2D transfer function widget has a 1D transfer function that helps further separate features within the voxels selected by the 2D transfer function. This dimension reduction method, however, can cause classification inconsistencies compared to a true 3D transfer function. We believe that this is a reasonable compromise, considering that the losses in classification precision compared to using an equivalent higher dimensional transfer function are typically minor. Rezk-Salama [79] proposed a similar idea called local transfer functions to set transfer functions for segmented volumes, i.e., a transfer function is associated with a tag in the tagged volume; voxels are essentially preclassified and their tags are stored in a volume. Our method is more flexible as the user essentially inter­actively labels voxels using the 2D transfer functions and then further classifies 22 the features using the associated 1D transfer function. Bruckner and Groller [9] similarly use a 1D transfer function to index into a table of style transfer functions for flexible illustrative volume renderings. Their work conceptually differs from ours as our transfer function combination method is utilized to improve the specificity of transfer functions rather than producing illustrative visualizations. 3.1.2 Selecting Combinations We propose to separate the transfer function space into a 2D transfer function space with a set of 1D transfer function spaces as a trade-off between dimension­ality and usability. A problem naturally arises when more than three properties/attributes are provided, namely which properties contain salient features, which attributes are most effectively used as the 2D transfer function domain, and which are best classified by the associated 1D transfer functions. Thus, we provide a few simple rules to aid the user in selecting appropriate combinations. For a given set l properties of a dataset, the correlation coefficient matrix R of size l x l is computed, as well as the entropy vector E of size l, which contains all properties' entropy. The primary property Yp, is chosen that represents the original information of the dataset (e.g., original intensity dataset or the mean dataset computed from the statistical properties extraction process as shown in Section 3.2.3). A property that is intrinsically associated with Yp (e.g., gradient magnitude vs. original intensity dataset or standard deviation vs. mean value) is used as the secondary property Yq. The primary and secondary properties define the 2D transfer function space. For all remaining properties Y{, i e [1, l] and i ± p, q a score is computed as a linear interpolation between the correlation coefficient E(Y) Rpi and the normalized entropy mxE, as shown in Equation 3.3: E(Y) Si = -aRpil + (1 - a)------ e / 0 - a < 1 (3.3) max E The correlation coefficient depicts the similarity between properties: a lower correlation coefficient value indicates a higher independence of properties. By intuition, more independent properties correspond to more interesting features, which we hope can be extracted by combining them together. Therefore, we favor 23 Rpi = (3.4) properties that are less correlated with the already chosen properties and as such a negative relationship between the absolute value of correlation coefficient |Rpi| and the score Si is shown in Equation 3.3. Specifically, the coefficient matrix R of property Yp and Yi is computed by Equation 3.4. Cov(p, i) ■\JCov(pr p)Cov(i, i) where Cov(p, i) is the covariance matrix of property Yp and Yi. However, using the correlation coefficient alone could lead to situations where properties that do not increase classification ability can beat more meaningful properties in the scoring, and to remedy this, the entropy of a property is also considered in Equation 3.3. The entropy value of a property reflects the amount of information contained in that property, shown as a normalized form mE(xYE) in Equation 3.3. The entropy is defined as n E(Yi) = Yu p(yb)log2(p(yb)) (3.5) b=1 where n is the number of bins in the histogram of property Yi, b is the current bin and p(yb) is the probability of data value yb at current bin. E(Y;) describes the homogeneity of property Yi and is negatively proportional to the homogeneity, i.e., higher entropy represents less homogeneity. Properties that are less homogenous usually contain more features of interest compared to more homogenous ones. Therefore, low homogeneity can be used to rule out less contributing properties that have higher score from the correlation coefficient. As such, high entropy is desired in our scheme, i.e., properties that are less homogenous are favored over more homogenous ones. However, low entropy may also be of interest on some occasions, e.g., a property contains large homogenous regions but highlights a small feature that no other properties can. The classification ability of those properties, however, is hard to describe by mathematical quantities but can be rather easily determined subjectively. The parameter a is dataset dependent and allows the user to choose a balance between the correlation of two properties and the amount of information contained in an individual property. 24 The remaining properties are then ranked based on their scores Si and used as the tertiary attributes for the associated 1D transfer functions. We found that using one or two tertiary attributes provides a good compromise between complexity and effectiveness of the classification. One of the available tertiary attributes is selected as the active one for each widget in the 2D transfer function space. As an example, the process of combination selection for CT chest scan Artifix discussed in Section 3.3.2 is shown. Using the rules, we compute the correlation coefficients and the entropies of the five properties of the dataset as shown in Table 3.1. Choosing the scalar value x as the primary attribute suggests using the gradient magnitude |Vx| as the secondary attribute. Then scores s^,a,p for mean, standard deviation and occlusion properties are computed for the remaining attributes by setting a to 0.4, which yields s^,a,p = [-0.0347,0.1938,0.3086]. The occlusion property has the highest score, meaning it is the best property regarding both the correlation between it and the primary attribute and the information it contains. The occlusion property is used as the tertiary attribute to define a combined 3D gradient magnitude/occlusion transfer function space. Alternatively, choosing the mean value p and the standard deviation a as the primary and secondary attributes, the scores sx,|vx|,p are computed for the other attributes, yielding sx,|vx|,p = [0.0094,0.2061,0.3045] for scalar, gradient magnitude and occlusion properties, respectively. The occlusion property has the highest score and is thus used as the tertiary attribute to define a combined 3D statisti­cal/ occlusion transfer function space. 3.1.3 User Interface In general, true 3D transfer function widgets are relatively difficult to interact with, since robust and effective interaction with a 3D space is still an open research problem [7]. The proposed combined transfer function space, however, is separable into a 2D transfer function space and a set of 1D transfer function spaces. Haidacher et al. [37] propose a similar separation method for multimodal visualization. In contrast to their simple triangle shaped windowing function, our method provides more insights and flexible controls for the 1D transfer function 25 Table 3.1. The correlation coefficients and the entropies of the properties computed from Artifix CT chest scan dataset. The properties are the intensity value X, the gradient magnitude |Vx|, the mean value p, the standard deviation a and the occlusion metric p. The bottom row shows the entropy E of each attribute. X |Vx| a P X 1.0000 0.1654 0.9973 0.2435 0.7286 |Vx| 0.1654 1.0000 0.1690 0.9569 -0.0583 0.9973 0.1690 1.0000 0.2464 0.7388 a 0.2435 0.9569 0.2464 1.0000 0.0067 P 0.7286 -0.0583 0.7388 0.0067 1.0000 E 5.3828 3.6077 4.8012 3.8392 7.9090 spaces. This separation, as stated before, can cause decreased classification precision when the 1D transfer function spaces are not independent from the 2D transfer function space compared to a true 3D transfer function space. However, our combination selection rules proposed in Section 3.1.2 help to rule out highly dependent 1D transfer function spaces. Therefore, we believe this separation is a good trade-off between interactivity and classification precision. Each 1D transfer function is attached to every selected region in the 2D transfer function domain based on the usual transfer function widgets or selectors. Thus, features in the volume can be classified by selecting their voxels in the 2D domain defined by the primary and secondary attributes. In cases where those voxels represent multiple separate features, the additional 1D transfer function can be used to further separate such features within the voxels selected in the 2D domain using one of the tertiary attributes. While adding complexity to the manipulation of transfer functions, this technique provides familiar interaction with each of the 2D and 1D transfer functions (TF). We believe this additional interaction (com­bining familiar 2D TF manipulation with familiar 1D TF manipulation) provides a reasonable method for interacting with the higher dimensionality of transfer function combinations. However, it does require users to be familiar with such interaction techniques. Figure 3.1 illustrates the proposed 3D transfer function editor for a 2D gradient magnitude transfer function space with associated 1D occlusion transfer functions. The top part shows the 2D gradient magnitude transfer function domain x x 26 Figure 3.1. The separable 3D transfer function editor. The transfer function editor with the 2D gradient magnitude transfer function space x x ||Vx|| shown on top, and the 1D occlusion space pc attached to the currently selected widget c in the 2D space, shown below. In this example, the blue widget is active, and the 1D histogram represents the occlusion information of all the voxels with statistical properties selected by the widget in the 2D statistical domain. ||Vx||, where the user can place and interact with traditional 2D transfer function widgets [56] TFw2d and a more generic lasso tool. The occlusion volume space pc or the size volume space tc of the region c selected by the currently active 2D transfer function widget TFwc is represented by a 1D transfer function editor, shown at the bottom, along with a 1D histogram of the occlusion information of all voxels selected by c. That is, the 1D transfer function editor operates strictly on voxels selected by a 2D transfer function widget (the blue widget in Figure 3.1). 27 The 2D transfer function widgets, such as ellipse, rectangle or triangle widgets as proposed by Kniss et al. [56], typically include some default shapes with few degrees of freedom. Users are able to set colors, opacities and different fall-offs for each of these widgets. These tools provide facilities to the user for a general exploration of transfer function spaces using easy to manipulate high-level widgets. However, it is difficult for the user to precisely select arbitrary regions. This often prevents a user from exploring the subtle structures in the transfer function domain, which may make a significant difference in the final visualization. Thus, similarly to commonly used image processing applications, we also include a lasso tool to allow the user to intuitively and easily select arbitrary regions by drawing the region boundaries directly into the transfer function space. In Figure 3.1, the red curve illustrates the hand drawn boundary path with a spherical fall-off for the color and opacity. A box on the left hand side of the 1D transfer function editor allows the user to select which tertiary attribute is used as the 1D transfer function space for each 2D widget. The proposed user interface allows the user to interact with the 3D transfer function space intuitively. Whenever the user creates a transfer function widget on the 2D transfer function space, the histogram of the voxels selected by that widget is computed and immediately shown in the 1D transfer function editor. Initially, the 1D transfer function maps, as visible, all voxels that are selected by the 2D transfer function widget. With the help of the 1D histogram, one can then design the 1D transfer function intuitively. As such, users are provided with a familiar interface thus providing intuitive interaction. This user interface adds minimal complexity to the standard 1D and 2D transfer function editors in existing volume visualization systems, e.g., Voreen [93] and ImageVis3D [48]. With a 3D transfer function space we are able to leverage the usability of the user interface; however, we are also interested in extending it for higher dimensional transfer function spaces in the future. 28 3.2 Specific Transfer Function Spaces Used In addition to the well-known 2D gradient magnitude and scalar value transfer function, we include several recently proposed transfer functions to be used in combinations. Creation of these transfer functions is generally based on the methods described in the respective papers, but with slight modifications, which are discussed in the following subsections: size-based transfer functions [14] in Section 3.2.1, occlusion-based transfer functions [15] in Section 3.2.2 and statistical transfer functions [38] in Section 3.2.3. 3.2.1 Size Information Computation Correa and Ma [14] proposed a three-step method to create a size volume S from an input volume. The three steps are scale-space computation, scale detection and back projection. Correa and Ma use anisotropic diffusion to create the scale space with better localization. The classical normalized Laplacian kernel is used to detect the blobs as local maxima both in spatial and scale domains. A back projection step utilizing Shepard's interpolation with Wendland polynomials is then conducted for the detected blob tuple (x, y,z, t). A single voxel can be part of features with multiple sizes; however, only the largest size value is kept at each voxel. Thus smaller features get masked out by larger ones, which happens in the brain MRI example shown in Section 3.3.4. To avoid this situation, we allow the user to specify an intensity range to compute a scale space specifically for that range. 3.2.2 Occlusion Information Computation Correa and Ma [15] suggest using an extended ambient occlusion metric to measure the occlusion of the volume. One can view the occlusion information p as a weighted sample mean value for a spherical neighborhood with certain radius R centered at each voxel, which results in an isotropic blurring effect that does not preserve the boundaries of the structures. Sometimes, overly smoothed volumes that lose all their boundary information are not desired, thus we derive a gradient based equation for computing the occlusion information, inspired by work done by Perona and Malik [76]. 29 For a sphere of radius R, we compute the occlusion information of the N voxels xi surrounding the current voxel x as shown in Equation 3.6: N mx = (3.6) i=1 2 n gx = llv/xl I2 + n2 <3,7) In Equation 3.7, gx is a term based on the gradient magnitude of the current voxel x. The dataset dependent parameter n G R+ handles gradients of zero magnitude e.g., for n G [0.001, 0.01], essentially helping to preserve boundaries of different structures. If n ^ 1, the filter behaves similar to a box filter. Computing mx is equivalent to convolving the volume with a spherical filter B r of radius R, and then modulating it with gx: mx = gx • <BR * Ix) <3.8) The complexity of this operation is O(mn), where m = |nR3 +1, and thus very costly, since the radius should be large enough to maximize the variance of the result [15]. This computational scheme is infeasible in practice, due to its computational complexity. However, since each sample inside the sphere is treated equally, a box filter of width 2R can be used to approximate the sphere. Exploiting the separability of convolving with a box filter and the performance of modern GPUs allows the computation of mx within seconds. The 3D convolution is then separated into three consecutive convolutions with a 1D box filter b2R+1 of width 2R +1, as Equation 3.9 shows: mx = gx • {b2R+1 * [b2R+1 * <b2R+1 * Ix)]} <3.9) This separation considerably reduces the computation time. Such an occlusion metric is view-independent and thus can be precomputed and stored, and there­fore does not affect the speed of visualization. 30 3.2.3 Statistical Properties We construct the statistical feature space with a procedure similar to that presented by M. Haidacher et al. [38]. They propose to grow a sphere over the neighborhood of each voxel and to compute the following statistical metrics: mean value ^, standard deviation o, skewness as well as kurtosis. It is a multistage process: first, extract statistical metrics and second, conduct the normality test. If the test is passed, continue with the similarity test. After the similarity test, if the new samples are similar to the old ones, we combine the statistical metrics. If any of the above tests fail or a user-defined maximum radius rmax is reached, the procedure is terminated, otherwise we increase the neighborhood by one voxel. Haidacher et al. [38] use the Jarque-Bera test [52] for normality since it is easily implementable on a GPU. It, however, requires a relatively large set of samples in order yield results of sufficient quality. Therefore, various other normality tests have been proposed in the literature; we chose D'Agostino's K-squared test [18] as a state-of-the art method. Its robustness with respect to identical values in the dataset makes it a good fit for CT and MRI datasets, which can contain large homogeneous regions. Utilizing the transformations Z1(^ b ) and Z2(b2) of the sample skewness Vb! and the sample kurtosis b2, the K-squared test (Equation 3.10) is then defined as: K2 = Z l (^ b ) 1 + Z2(b2)2 (3.10) K2 is approximately x 2-distributed with 2 degrees of freedom; we can test its null hypothesis by looking up the x 2-distribution table. The entry for test level 1 - a = 0.999 with 2 degrees of freedom in the x 2-distribution table is 13.82. Therefore, the normality test will be passed if K2 < 13.82 (3.11) If the samples in the spherical neighborhood pass the normality test, it is necessary to further test whether they have the same distribution as that of the samples computed in the previous iteration. As done by Haidacher et al. [38], Welch's T-test [98] is used to compare the similarity of the sample distributions. 31 3.3 Results and Discussion The statistical properties, the occlusion information and the size information are all precomputed on the GPU, and those volumes are then used in the interactive visualization stage to define the transfer function space. Users interact with an ex­tended slice-based volume renderer implemented in OpenGL and Qt that supports combined 3D transfer functions to explore and generate final visualizations. The scoring process is not part of our volume renderer and is conducted in MATLAB only once for a dataset. The input is a matrix where each of its columns is a property volume that is flattened into an 1D array. The correlation coefficient matrix is computed by the MATLAB function corrcoe f, which uses Pearson's correlation. The MATLAB function entropy, which implements Equation 3.5, is applied to compute the entropy of each property by taking the histograms of the properties in a column of the input matrix. The number of bins of the histograms is determined by the number of bits of the data, e.g., an unsigned 8 bit volume has 256 bins. Finally, Equation 3.3 is evaluated for the corresponding row of the major property in matrix R and the normalized entropy vector mY E . The whole process takes about 10 seconds for each of the examples shown. The following discussion compares 2D gradient magnitude, 2D statistical, 2D occlusion, 2D size with 3D combined statistical/occlusion, statistical/size, occlu­sion spectrum/size or statistical/(occlusion, size) transfer functions applied to a synthetic dataset and real-world datasets. The combined 3D transfer functions for each dataset were typically designed within 15 to 20 minutes, similar to the time required to design the traditional 2D transfer functions. The synthetic dataset models a filled shell encompassing varying sized spheres; the "Artifix" dataset has been retrieved from the OsiriX DICOM archive [72]. The back pack and the "Artifix" datasets are CT scans of a back pack and chest, respectively, "CerebrixCrop" is the T1 channel of an MRI scan of a brain. The parameters used to create the transfer function spaces are chosen by trial-and- error on each dataset. The synthetic dataset and CT datasets are computed with confidence level 0.1, whereas the MRI dataset with a confidence level 0.001 when generating the statistical transfer function space. The radius is set to 40 for all 32 datasets when creating the occlusion volumes. The synthetic dataset is processed with a boundary preserving parameter n = 1.0 in order to overcome the noisiness, whereas all other datasets use n = 0.005 to preserve the boundary details. The size property computed for the MRI dataset is limited to the intensity range [250,500] in order to classify the tumor. The transfer function combinations shown below are chosen by applying the algorithm described in Section 3.1.2 with varying parameter a. An exception to this are the results shown in Figure 3.2, where the extremely noisy nature makes Equation 3.3 ineffective. © ©i (a) (b) ■-* i k . (c) (d) (e) (f) (g) (h) Figure 3.2. Comparisons of transfer function combinations for a synthetic dataset. The synthetic dataset was created as a mixture of overlapping Gaussian distri­butions with varying parameters to model a filled shell encompassing varying sized spheres as shown in a). The dataset has been classified, from left to right, using transfer functions (shown right below the rendered images) based on b) 2D gradient magnitude, c) 2D occlusion spectrum and d) 2D size-based transfer function, e) 2D statistical, f) combined statistical/occlusion transfer function, g) combined statistical/size transfer function and h) combined statistical/(occlusion, size) transfer function. 33 3.3.1 Synthetic Dataset A synthetic dataset was created, as illustrated in Figure 3.2(a), in order to mimic a common scenario in real life medical datasets, such as chest CT scans or head MRI scans, where different structures overlap both spatially and in the scalar values. Often, the outer structures occlude with the inner ones, but they also can have different sizes. The synthetic dataset contains six different materials: the environment with p0 = 0.20,a0 = 0.14, the middle hull with ^1 = 0.40,a1 = 0.16, the outer hull and the upper small inner sphere with p2 = 0.60,a2 = 0.11, and both the remaining larger and smaller inner spheres have ^3 = 0.80,a3 = 0.13. In addition, low amplitude noise following a Gaussian distribution has been added across the whole domain to simulate noise introduced by acquiring a volumetric image with a scanner. Various transfer functions have been applied to the synthetic dataset, as shown in Figure 3.2. Traditional 2D gradient magnitude-based transfer functions, as Figure 3.2(b) illustrates, suffer severely from the overlapping scalar values in the transfer function domain. There, features are indistinguishable due to noise, which makes it hard to separate features based on their gradient magnitude, as seen in the joint histogram in Figure 3.2(b). Occlusion spectrum 2D transfer functions, shown in Figure 3.2(c), are able to separate the inner and outer structures based on their occlusion property as in the transfer function shown in Figure 3.2(c). The three inner spheres, however, cannot be separated clearly due to the similarity in their occlusion information as well as their scalar values. Also, the center of the inner yellow region overlaps with all spheres in the occlusion spectrum, thus causing misclassification. The size-based 2D transfer function applied to the dataset (Figure 3.2(d)) separates the inner spheres from each other and the outer rings; however, there are classification artifacts at the top and right part of the green outer ring. The small sphere at the bottom right cannot be properly separated from the purple sphere, since they both overlap in their scalar values. Statistical 2D transfer functions, as demonstrated in Figure 3.2(e), are able to separate the overlap in the (p,a) transfer function domain. It is thus possible to 34 classify them using different properties. However, both spheres at the lower center have the same statistical properties, and similarly, the outermost shell shares the statistical properties with the upper central sphere, yet they represent different structures. Supplementing the statistical information with occlusion information, as shown in Figure 3.2(f), makes it possible to separate the inner purple sphere, compared to Figure 3.2(e). The transfer function in Figure 3.2(f) shows that the 1D occlusion histogram for the highlighted 2D widget can be used to separate the purple sphere with its low amount of highly occluded voxels from the green outer shell, which has a higher amount of less occluded voxels. However, the two red spheres at the bottom are not separated from each other. On the other hand, supplementing the statistical information with size infor­mation, as shown in Figure 3.2(g), makes it possible to separate the two spheres at the bottom into the cyan small one and the larger red one, when compared to Figure 3.2(e). Noticeable are the purple artifacts in the green outer shell at the right side, since that region has a similar feature size compared to the purple sphere. Figure 3.2(h) shows that occlusion and size information together are able to classify all the features of the dataset without ambiguity. The 1D transfer function associated with each widget in the 2D statistical transfer function space uses either size information or occlusion information to further classify the voxels selected in the statistical 2D transfer function domain, thus allowing the user to exploit the benefits of either method, while being able to interact with 1D and 2D transfer functions, instead of 3D or 4D transfer functions. 3.3.2 CT Scan of a Chest: "Artifix" In the chest CT scan "Artifix" (Figure 3.3), both traditional 2D and combined 3D transfer functions were used to classify the lung (blue), bones (shades of gray), blood vessels (red), aorta (dark orange), kidney (brown) and the skin (transparent gray). The gradient magnitude transfer function fails to correctly separate the blood vessels and the kidneys from the bones. Also noticeable is the relatively high amount of noise distributed across the volume. 35 (a) (b) (c) (d) (e) Figure 3.3. Transfer function combinations for a CT chest scan. The chest CT scan "Artifix" classified using transfer functions based on a) 2D gradient magnitude, b) 2D occlusion spectrum, c) 2D statistical, d) combined 3D gradient magnitude/occlusion transfer function and e) combined 3D statistical/occlusion transfer function. The occlusion spectrum can be used to separate the kidney from the surround­ing tissue. However, the aorta is similarly classified, since the aorta and the kidney are overlapping in the occlusion spectrum. Also, details of the lung are lost, since its tissue has similar occlusion values compared to the surrounding tissue, due to the intricacy and delicacy of the alveoli and bronchioles. A statistical transfer function (Figure 3.3(c)) removes a noticable amount of that noise, but still leaves some areas, such as the front part of the ribs, and the kidney misclassified, since they are close with respect to their statistical properties. Experimentation with the size-based transfer function as the associated transfer function space did not measurably improve the classification since the relative similarity of the scalar values in this CT scan mapped them to similar size values. However, combining occlusion information with either a 2D gradient magni­tude transfer function (Figure 3.3(d)) or a statistical transfer function (Figure 3.3(e)) 36 increases the ability to correctly separate the kidneys from the aorta. The fine structures of the lung's surface are identifiable, since they have different statistical properties compared to the surrounding tissues. There are only slight differences between the combined transfer functions since they are similar without consider­ing occlusion information. 3.3.3 CT Scan of a Back Pack Figure 3.4 shows the CT scan of a back pack filled with liquids <in red, green, blue), a battery <in purple) and a box <in cyan) classified with various transfer functions. The scoring with a = 0.6 conducted on the back pack dataset with scalar value chosen as the main property and gradient magnitude as the intrinsically associated secondary property results in: s^,a,p,S = [-0.5208,-0.3882,0.2165, -0.1595] <a) <b) <c) <d) <e) <f) <g) Figure 3.4. Transfer function combinations for a back pack CT scan. The back pack CT scan classified using transfer functions based on a) 2D gradient magnitude, b) 2D statistical, c) 2D size, d) 2D occlusion, e) combined 3D gradient magnitude /occlusion transfer function, f) combined 3D statistical/size transfer function and g) combined occlusion spectrum and size transfer function. 37 suggests that the occlusion volume p and size volume S should be considered for tertiary attributes. Changing the main property to mean volume with standard deviation volume as the secondary attribute gives the scoring for the rest of the properties: s*,|v*|/P/s = [-0.4847, -0.2140,0.2153, -0.1679] also hints to us that the occlusion volume p and size volume S should be used as tertiary attributes. The 2D transfer functions separate the different liquids to varying degrees, but they fail to identify the battery properly. Other features, such as the wires or the small circular shapes, are mapped to the same color yellow (Figures 3.4(a), 3.4(b)), or the same feature is mapped to different colors (Figures 3.4(c), 3.4(d)). Notable is the 2D occlusion transfer function, which allows the extraction of the cyan box but classifies the liquids with less specificity. Adding occlusion as the third axis did not yield meaningful results, since the dataset itself has many features that are similarly occluded by the clothing articles (showing in transparent gray) inside the back pack, thus reducing the separability in the occlusion channel. Utilizing a size transfer function as the third axis allows the clear separation of the battery (purple color). The 3D occlusion spectrum/size transfer function (Figure 3.4(g)) is additionally able to visualize the cyan box, which is difficult to do using gradient magnitude (Figure 3.4(e)) and statistical information (Figure 3.4(f)) as the 2D transfer function domain. However, all the 3D transfer functions have problems in classifying the wires as features both connected spatially and with respect to their colors, suggesting further investigations of alternative volumetric attributes as the third axis. 3.3.4 MRI Scan of a Brain: "CerebrixCrop" MRI datasets, occurring in clinical and research studies where separating the brain from the surrounding tissue is of particular interest, are typically challenging to classify, since they often contain ubiquitous noise [27]. Figure 3.5 shows such a dataset containing a tumor in the center of the brain. Transfer functions are 38 Figure 3.5. Transfer function combinations for an MRI brain scan data. The "CerebrixCrop" MRI dataset shown with focus on the brain tissue shown as yellow (top row) and a tumor shown in red (bottom row). The following transfer functions were applied: a,f) 2D gradient magnitude, b,g) 2D occlusion, c,h) 2D statistical, d,i) 3D statistical/occlusion, e,j) 3D statistical/(occlusion,size). applied to classify the brain tissue (in yellow) and the fluid inside the tumor (in red). Note that although both features can be shown simultaneously by setting transparency of the brain, we set the brain to be completely transparent in the second row of images for clear visualizations of the tumor. We apply the scoring process with a = 0.6 to the MRI dataset: set scalar value as primary and gradient magnitude as secondary yields: sF,o,p,s = [-0.3188,-0.1247,0.1674, -0.0224] meaning that the occlusion volume p and size volume S once again should be considered for tertiary attributes. Substituting the main attribute with the mean value with the standard deviation as the secondary attribute gives sX/|VX|,P/S = [-0.2841,-0.0591,0.1505, -0.0298] and leads us to the same decision. Gradient magnitude-based 2D transfer functions (Figure 3.5(a)) fail to properly separate the brain from the skin, since they both share similar ranges of scalar values and gradient magnitudes. 39 Figure 3.5(f) demonstrates the inability for the gradient magnitude based 2D transfer functions to clearly pull out the tumor, since similar scalar values and gradient magnitudes appear universally across the dataset. The occlusion spectrum (Figures 3.5(b) and 3.5(g) helps to better separate the brain from its surrounding tissues as well as remove noise with scalar values similar to the tumor. However, the surface of the brain tissue is still incorrectly classified and a large amount of noise still appears around the tumor due to similar occlusion values in these regions. Statistical transfer functions (Figures 3.5(c) and 3.5(h)) significantly smooth the dataset, making the creases and recesses of the brain tissue clearly show up, however, noise that heavily affects the visual quality is still seen across the dataset, especially in Figure 3.5(h). Combining the occlusion information with statistical information, as shown in Figures 3.5(d) and 3.5(i), classifies the brain tissue properly, but fails to clearly extract the tumor. However, a transfer function combination with two tertiary attributes, as shown (Figures 3.5(e) and 3.5(j)), clearly separates both the brain tissue and the tumor. The statistical attributes are used as the primary and secondary attributes, and the occlusion and size information are used as the tertiary attributes. The widget that classifies the yellow brain tissue uses the occlusion attribute to further classify it with the associated 1D transfer function; however the widget classifying the red tumor uses the size attribute instead to further remove the noise via its associated 1D transfer function. 3.3.5 Multivariate Dataset: Hurricane Isabel One time step (time step 30) of the VisContest 2004 Hurricane Isabel [47] multivariate dataset is used to demonstrate the generality of our method. The dataset is a simulation of a hurricane from the National Center for Atmospheric Research in the United States. The original dataset contains 100 time steps and each with 12 attributes. Many of these attributes, however, are redundant or contain little amount of information. 40 The three most salient attributes are selected, namely pressure, temperature and a water vapor mixing ratio measurement, QVAPOR, by evaluating the entropy of each attribute. The 2D transfer function domain is pressure and temperature. We use QVAPOR as the associated 1D transfer function. These attributes are then used to classify significant features in meteorology, including hurricane eye and spiral arms. As shown in Figure 3.6, each colored widget in the 2D domain uses a different QVAPOR 1D transfer function for classification and the classified volume are visualized using volume rendering. The eye of the hurricane (shown in red) has a lower pressure but higher temperature than the blue outer bands and a lower temperature compared to the yellow and green spiraling bands. The QVAPOR attribute allows us to see the spiraling bands with fine details in the dataset. (a) Hurricane Isabel Visualization (b) Hurricane Isabel UI Figure 3.6. Classifying a multivariate dataset using the combined transfer function space. Visualization of the multivariate Hurricane Isabel dataset using pressure and temperature in the 2D transfer function with different 1D transfer functions using QVAPOR for each 2D transfer function widget shown in different colors. CHAPTER 4 TRANSFER FUNCTION DESIGN BASED ON USER-SELECTED SAMPLES FOR INTUITIVE MULTIVARIATE VOLUME EXPLORATION Multivariate volumetric datasets are important to both science and medicine. We propose a transfer function (TF) design approach based on user-selected samples in the spatial domain to make multivariate volumetric data visualization more accessible for domain users. Specifically, the user starts the visualization by probing features of interest on slices and the data values are instantly queried by user selection. The queried sample values are then used to automatically and robustly generate high-dimensional transfer functions (HDTFs) via kernel density estimation (KDE). Alternatively, 2D Gaussian TFs can be automatically generated in the dimensionality reduced space using these samples. With the extracted features rendered in the volume rendering view, the user can further refine these features using segmentation brushes. Interactivity is achieved in our system and different views are tightly linked. Use cases show that our system has been successfully applied for simulation and complicated seismic datasets. 4.1 Method Overview The workflow of our proposed method as shown in Figure 4.1 is comprised of three major stages: (A) data probing, (B) qualitative analysis and (C) optional feature refinement. Data probing is the process where the user discovers regions of interest by examining multivariate data slices. The regions of interest can be conveniently selected using a lasso tool or a "magic wand" tool. Once the regions of interest are selected, a simple, yet efficient, voxel query operation that inquires the multivariate data values is performed. The user then performs a qualitative analysis, i.e., Input T Multivariate Volume Dataset 1 Voxel Q u e r y (A) Data Probing -^ ----- r V 2D Gaussian TFs Creation Output KDE based Fine m HDTF Creation Tuning*-------- - 1 J ) Figure 4.1. The user interface and the work flow of the system implementing our proposed method. Four closely linked views are shown and labeled, namely: (1) multipanel view, (2) volume rendering view, (3) projection view and (4) high-dimensional transfer function view. Three stages: (A) data probing, (B) qualitative analysis and (C) optional feature refinement comprise our work flow. With the proposed method and user interface, domain users are able to explore and extract meaningful features in highly complex multivariate dataset, e.g., the 3D seismic survey shown above. K) S e gm e n ta tio n L a sso o n 3D 1 S eg Bru sh o n 3D S eg Bru sh on S lice s (C) Feature Refinement (1) Multi-Panel View (2) Volume Rendering View (3) Projection View (4) High-Dimensional Transfer Function View 43 extracting and rendering volumetric features by means of designing HDTFs or 2D TFs on dimensionality reduced spaces. KDE is utilized to automatically generate the HDTFs and to robustly discard outliers from the queried samples. Similar to previous methods [1,5, 31], the transfer functions can also be directly modified in the high dimensional transfer function editor. The high dimensional data space as well as the transfer functions are represented by parallel coordinate plots and pairwise scatter plots. In addition, automated 2D Gaussian TFs on the projection view offer a simpler alternative for more distinct features. The HDTFs can then be fine-tuned directly in a PCP-based HDTF editor while the 2D Gaussian TFs can be manipulated by 2D Gaussian TF widgets. On many occasions, however, different features share similar data values and thus an optional feature refinement stage is introduced to refine the features classified by the TFs. Features are refined by the user via segmentation brushes or lassos that are applied directly on the volume rendering view or the multipanel view. After several iterations of the three stages, the user can choose to output the classified result as labeled volume for further processing. We have implemented an interactive multivariate volume visualization system based on the proposed method that has been implemented to allow domain users to extract refined features in very complicated multivariate volume datasets more intuitively. 4.2 Voxel Query and PCP Generation Our proposed method is based on user-selected multivariate voxel samples through interactive selection, which requires efficient voxel query. The multivari­ate values of the queried samples should be immediately presented to the user by means of PCPs, and so a fast PCP generation method is needed. 4.2.1 GPU-based Voxel Query via Conditional Histogram Computation Voxel query can be accelerated by spatial hierarchy structures that group similar neighboring voxels into nodes, e.g., an octree structure adopted by Guo et al. [31]. However, Knoll et al. [58] report that, "Conversely, volumes with uniformly high variance yield little consolidation; due to the overhead of the octree hierarchy they 44 could potentially occupy greater space than the original 3D array." Our initial experiment on the seismic data with the code from [58] agrees with this statement. We therefore propose to efficiently conduct the voxel query by computing sets of joint conditional histograms via a simple GPU-based volume traversal. A joint conditional histogram jch<a, b)f of two attributes a and b is a 2D histogram showing the joint distribution of attribute values Ya and Yb of voxels V whose evaluated result from a certain boolean function f<Y(V)) <Y(V) being the attribute values of V) is true. If f is always true, the joint conditional histogram degenerates to an unconditional joint histogram. Note that the values of user selected samples are queried via an unconditional joint histogram computation over the user-selected region on the given slice. For a multivariate volume of N attributes, given an N-dimensional TF as the condition, a set of N - 1 joint conditional histograms can be computed to record the query results. The values of the joint conditional histograms are accumulated by first evaluating the N-dimensional TF for all voxels in the volume, and then transforming the voxels that have positive opacities from the TF into bins in the conditional histogram space, and finally incrementing the joint conditional histogram count at those bins. Specifically, given a voxel vX of N attributes Y1, Y2,...,Yn <to be concise, we use yi to denote the attribute value Yi<vX)) located at 3D position X in the spatial domain, and an N-dimensional TF TF. vx ^ {< y1, y2),<y2, y3), - , <yN-1, yN)} where TF<y1,y2, . . . , yN).a > 0 <4.1) <y1,y2) ,<y2, y3),...,<yN-1, yN) being the bins of joint conditional histograms jch<Y1 , Y2), jch<Y2, Y3) ,..., jch< Yn-1, Yn ), respectively. Equation 4.1 and the accumulation of the conditional histograms, which are stored aggregately as a 2D texture array of N - 1 slices, can be easily implemented on the GPU via geometry shader and ADD blending or read-write textures with atomic operations that are supported on recent GPUs. 45 4.2.2 Parallel Coordinate Plots Generation As proposed in [70], Figure 4.2 shows that each nonzero pixel P(i, j) in the joint histogram of attribute x and y yields a quad starting at the position of i on PC axis x and ending at the position of j on PC axis y. The highly parallel process can be implemented on the GPU using geometry shader and transform feedback buffers. The algorithm loops through all pairs of conditional histograms after setting up the transform feedback buffer for recording the resulting geometry. In each iteration, a regular grid of the same size of a slice of the input conditional histogram texture texcond is drawn and a geometry shader generates a colored quad for each vertex whose texcond value is not 0. The dynamic range of the data values is usually high and thus the ratio of the natural logarithm of the data value versus natural logarithm of the total voxel number is computed and then modulated with the input color C0(i, j) at grid position (i, j) to give the final color C(i, j). j)) C fa j) j) b g (£ v) (4.2) Finally, all quads are stored in the transform feedback buffer, and they can be rendered directly from the transform feedback buffer without being read back to the CPU. 4.3 Transfer Function Generation from User-Selected Samples In this section, the actual TF generation method will be explained. Section 4.3.1 introduces the method for interactive voxel sample selection, Section 4.3.2 dis- Figure 4.2. Generating a PCP from a joint histogram. 46 cusses the KDE-based HDTF generation method and Section 4.3.3 presents details on the automated 2D Gaussian TF on the dimensionality reduced space. 4.3.1 Sample Selection in the Multipanel View The user can interactively select an arbitrary region of interest in any attribute by either drawing a lasso or using the magic wand tool. The lasso tool is a simple free hand drawing tool that allows the user to select regions by manually drawing over the boundary of a feature. Although very flexible, the user has to be very careful when drawing on the boundary using the lasso tool. To alleviate the difficulty of perfectly drawing over the boundary of a feature, a more intuitive and easier to use magic wand tool is introduced. The magic wand tool is essentially a 2D segmentation tool based on Perona and Malik's anisotropic diffusion [76]. Equation 4.3 describes the diffusion equation where S(t,X, y) is the number of seeds at position (x, y) at time t, with V(t,x, y) being the intensity of the chosen attribute at the same point, |VV(t,x, y)| is its gradient magnitude, and g(s) a conductivity term. dS(t, X, y) dt = div(g(|VV(t, x, y)|)VS(t, x, y)) (4.3) where g(s) = v • exp K2 Parameter K governs how fast g(s) goes to zero for high gradients, regular term v is chosen as 1 and normalization term h is set to n+y for numerical stability, n being the number of neighbors of a pixel, which is 8 in our case. Equation 4.3 can be solved numerically using the finite difference method with a given iteration number T. The iteration number T, parameter K and seeding brush size are user controllable. Figure 4.3 shows the panel view of a six-attribute seismic volume dataset where attributes are co-rendered with the seismic amplitude volume. Note that a user drawn magic wand in dark blue highlights a potential salt dome structure. 4.3.2 Kernel Density Estimation-based Transfer Function Generation We would like to generate HDTFs from the samples selected using method described in Section 4.3.1. To reduce the computational complexity, we separate -s 47 Figure 4.3. Sample selection in the multipanel view. The user draws on a salt dome (stroke shown in light blue) over the fifth attribute in the panel view, resulting in the dark blue region of selection. the N-dimensional value space into N -1 2D value spaces, i.e., a 2D + 2D +------+ 2D (N - 1 of 2D) space. A naive approach is to generate a TF by taking the convex hull of these 2D sample points. Although useful when the user intends to select exact sample points, it is conceivable that the outliers in the samples can greatly bias the generated TF and result in unwanted regions selected in the value space. Figure 4.4(a) clearly demonstrates such a situation where a red 2D TF widget is generated as the convex hull of the sample points with the red boundary. Also notable is that the color gradient of the TF widget is arbitrarily defined by the user that may not follow the underlying distribution of data. Kernel density estimation (KDE) [87] seen in Equation 4.4 is a nonparametric method for estimating the density function fh(x) at location x of an arbitrary dimensional domain Q with given samples {xi},i e {1, 2, 3 , . . n}. 1 n 1 n _ f h(x) = n Y j Kh(X - Xi) = n k ^ K( )r X' Xi e Q (44) i=1 i=1 where K(x) is the kernel function and h is the bandwidth. Thanks to the separation of the value space, instead of computing the KDE for Q of N dimension, we 48 (a) (b) (c) Figure 4.4. Transfer function generated from user-selected sample points. User-se­lected sample points (shown in green) over a joint histogram. TF widget generated from the samples as (a) convex hull and (b) KDE. In (c): a point cloud (left) and its KDE result color coded with a "jet" color map. compute N - 1 KDE for Q in 2D spaces. In our case, each Q is set to the same size of the 2D joint histogram, which is typically 256 X 256. An empirical optimal bandwidth estimator is suggested in [87], which can be extended to 2D: h = 1.06 VdetE • n- 5 (4.5) where det E is the determinant of the 2D covariance matrix E of current attribute pairs. The kernel function K(x) we used is the 2D Gaussian kernel: 1 ||x ||2 K(x) = ,__e 2 (4.6) V2n With the Gaussian kernel, each sample xi contributes to the estimate in accordance with its distance from X. Therefore, in the region near the intended samples more short distanced samples are contributing to fh(x) compared to the region near the outliers. As a result, the density value fh(x) around the outliers is lower than that of the intended samples. Figure 4.4(c) shows the density function generated by the KDE method of the given samples with the above settings. This result verifies our expectation that the outliers have a lower density than the intended sample regions. We can discard the outliers by setting a threshold for the density value fh(x). Figure 4.4(b) shows the yellow TF widget generated by KDE with a density threshold of 0.15. Noticeable is that the outliers are excluded from the TF widget and the smooth color gradient that actually follows the underlying density. The 49 resulting TF can be represented by a set of 2D TFs or a PCP created using the method described in Section 4.2.2. In the presence of multiple HDTFs, ambiguity could arise: different HDTFs can cover the same regions of certain 2D attribute pairs. To differentiate the HDTFs, a unique ID is specified to each HDTF and an ID map of the same size of the N - 1 2D TF space is created by conducting bitwise OR for all HDTFs on each 2D attribute pair. The ID map is later decoded in the volume rendering shader to correctly select voxels. 4.3.3 Automated Gaussian Transfer Functions on Dimensionality Reduced Space Dimensional reduction is another popular method for visualizing high-dimensional data due to its ability to intrinsically generate visual representations that are easy to understand and interact with. Instances in an m-dimensional Cartesian space are projected into a lower p-dimensional visual space with preservation of the distances between instances as much as possible. In other words, voxels with similar m-dimensional attribute values are projected to be near each other in the p-dimensional space. With a projected visual space of p = 2, the user is able to better identify features by doing visual classification using a 2D TF widget, and moreover, automated clustering methods can be applied for classification. In our proposed method, the high-dimensional value space is projected into a 2D space using Fastmap [24] and then Gaussian TFs are generated via expectation maximization optimization with Gaussian mixture model. The user can choose to use either the 2D Gaussian TF or the HDTF for each feature by switching a button on the user interface. The 2D Gaussian TFs are preferred for more convenient extraction of several distinct features at the same time, whereas the HDTFs are better for features that have subtle differences in the HD value domain. We employ Fastmap [24] as the dimensional reduction technique since it is fast, stable and easy to implement. Fastmap is a recursive algorithm for multi­dimensional projection with an O<N) time complexity. Given target dimension k, a distance function D<) and object array O contains N objects of m dimension, the algorithm FastMap computes the fc-dimensional projected image X from the N 50 objects. The algorithm is summarized in Algorithm 4.1 Assuming that all attributes we are handling are continuous measurements, the dimensionality reduced 2D value space can be modeled by a Gaussian mixture model <GMM). GMM models point clouds by assigning each cluster a Gaussian distribution. For a point x in the 2D value space, a Gaussian distribution is shown in Equation 4.7 with mean value ^ being a 2D vector and covariance matrix E as a 2 x 2 matrix. N ^ ' E) = c- 2(x-p)^E-' (x-p) <4.7) Therefore, for a GMM with fc components, the distribution of the 2D value space can be written as fc p(x\0) = ^ a jN(x\j, Ej ) (4.8) j=1 where 0 is the parameter set of the fc-component GMM {aj ,^ j , Ej}j=v and a j is the prior probability of the jth Gaussian distribution. The optimal 0 can be found as d that maximizes the likelihood of p(X|0) n 0 = argmax p(X|0) = argmax f l p W ) (4.9) i=1 where n is the number of input points. Equation 4.9 can be solved by the expectation maximization (EM) algorithm [4]. Given an initial setup of 0, the EM Algorithm 4.1 FastMap(fc, D<), O) if fc < 0 then return else col = col +1 (col is initialized to 0) end if Choose and record the pair of pivot objects Oa, Ob. Project objects on line (Oa, Ob) using the cosine law: X[i,col] = Xi = D(OaA)2+DDOOOO)b2)-D<°bA)2, i 6 {0,1,2,...,N - 1} Call FastMap(fc - 1 ,D'<),O). Where D'(O',O')2 = D(Oi,Oj)2 - (xi - xj)2,i, j 6 {0,1,2,...,N - 1} 51 algorithm iterates between two steps: expectation step (E step) and maximization step (M step) until the log likelihood n n k ln p(X|0) = log(^| p(xlld)) = ^ { ^ a jN (xV , Ej )} i=1 i=1 j=1 converges. We initialize the EM algorithm using the K-means algorithm [40], which quickly gives a reasonable estimation of 0. With an initialization of k mean values {uj}k=1, K-means algorithm iteratively refines {uj }k=1 until convergence through assignment and update steps. The assignment step assigns each sample to the cluster with the closest mean, and the update step calculates the new means to be the centroid of each cluster. In our case, the initial means are k random samples in the input dimensional reduced 2D point cloud. Once the K-means algorithm terminates, {Ej }k=1 can be easily computed with the result means, and prior probabilities {a j}k=1 is given by the proportion of total samples inside each cluster. We use a modified TF generation scheme as in [96] but ours differs in that 1) the value space we use is the 2D dimensionality reduced space of high-dimensional attribute compared to the 2D intensity versus gradient magnitude space as in [96], and 2) we use the user-selected samples as the input point clouds, whereas they use all voxels in a volume. Given some user-provided sample data points and a class number k (which is set to 3 by default from our experiments), the EM algorithm computes the Gaussian distribution parameters 0. Each Gaussian distribution is managed by a Gaussian TF widget with a user-defined color C and opacity function a of location x: a = amax^- 2(x-U)T E-1(x-U) (4.10) The Gaussian TF widget is centered at the mean value u of the Gaussian distribu­tion and its boundary is generated by transforming a unit circle with the square root matrix E1/2 of covariance matrix E. E1/2 is calculated via eigen decomposition of E: E = VDV-1 (4.11) E1/2 = VD1/2 V-1 (4.12) 52 where D is a diagonal matrix holding the eigenvalues and V contains the eigen­vectors as columns. V is an orthogonal matrix, i.e., V-1 = VT, since E is symmetric. The eigenvalues 0 1 ,0 2 are the radii of the principal axes of the ellipse, whereas the eigenvectors a, b are the unit vectors of the principal axes. Transformations of the Gaussian widgets, i.e., translation, rotation and scaling, can be achieved using the eigenvalues and eigenvectors. The translation is done by shifting the ^ with an offset A^ given by user dragging. The rotation of the widget is achieved by rotating the eigenvectors in V with an angle j6. Finally, multiplying the eigenvalues 0 1 ,0 2 with a scaling factor (sa,Sb) results in the scaling of the widget. 4.4 Feature Refinement in the Spatial Domain The feature refinement stage is introduced to allow the user to directly ma­nipulate the features in the spatial domain. Various refinement tools have been implemented to handle different situations. All tools support three refinement modes: new, add and remove. 4.4.1 Screen Space Brush in the 3D View. The tool as seen in Figure 4.5(a) allows the user to draw strokes on the 3D view screen to set seeds in the visualization results, and then a GPU-based region growing is conducted to set the connected voxels to a given tag number. The seeding location is determined by casting rays from brush strokes on the image (a) (b) (c) Figure 4.5. Feature refinement tools. Three feature refinement tools are included in our method: (a) 3D brush, (b) 3D lasso and (c) 2D brush. 53 plane to the volume extracted by current TFs. A voxel along the ray is seeded when its opacity is greater than a user-defined threshold. 4.4.2 Screen Space Lasso in the 3D View. Alternatively, the user can directly indicate features of interest on the 3D view using a lasso as shown in Figure 4.5(b). A lasso is a simple tool that selects all voxels from the TF extracted volume that are inside the back projected volume of the screen space lasso covered area. 4.4.3 Refinement Brush in the Panel View. The refinement can also be done by seeding on the panel view via drawing strokes (Figure 4.5(c)), and this is useful when the features of interest are occluded in the 3D view or readily visible in a slice. A morphological closing, i.e., dilate the volume by one voxel and then erode the volume by one voxel, is performed after refinement in order to fill small holes and bridge tiny gaps. Note that all refined feature groups are managed in the group manager in the HDTF editor introduced in Section 4.6.2, and thus similar to TF groups, their colors can be changed, they can be deleted and their visibility can be toggled. 4.5 Rendering We employ the directional occlusion shading (DOS) [84], which is an efficient approximation to ambient occlusion as the rendering technique because the DOS is gradient-free and provides the user more insights into the dataset than local shading models as shown on seismic datasets [73]. A user study conducted by [61] shows that DOS outperforms other state-of-the-art shading techniques in relative depth and size perception correctness. Hardware supported trilinear interpolation cannot be used for tag volume rendering because false tag values will be generated. Instead, nearest neighbor sampling has to be used to correctly render the tag volume. However, a simple use of nearest neighbor sampling yields blocky looking results because of the voxel level filtering. Instead, using a manual trilinear 0-1 interpolation gives pixel level filtering. From our observations, the cases where multiple tags appear in a single 8 voxel neighborhood rarely occur 54 and so a simplified method of [34] is utilized. The largest tag value in the eight neighboring voxels around current pixel is mapped to 1 and all others to 0 and then a trilinear interpolation is conducted on these 0/1 values. The interpolated result is then compared against 0.5. If greater, the final tag value of the pixel is set to the pixel's nearest neighboring voxel's tag value, otherwise the tag value is set to 0. 4.6 User Interface The user interface of our system is seen in Figure 4.1 where a multipanel slice view for data probing is shown to the left (1), an interactive 3D view that shows volume rendering results and allows post feature manipulation is seen in the middle (2), a projection view shown to the upper right (3) and a high-dimensional transfer function view appears to its bottom (4). These four views are tightly linked and any updates in one view will be reflected in others. 4.6.1 Multipanel Viewer We have developed a multipanel view that shows all attributes of a slice by placing attributes into individual panels as seen in the left part of Figure 4.1 as well as in Figure 4.3. The multipanel viewer synchronizes user interactions across all attribute views, including: mouse positioning, panning, zooming, scrolling and aspect changing. To enhance the perception of attributes, each attribute can have a specifically designed color map that highlights features of interest. In order to better use the dynamic range of the color maps, the contrast of the attributes can be conveniently changed using the mouse wheel. Furthermore, a background volume can be co-rendered with the current attribute volume using transparency. This rendering mode is especially helpful for seismic volumes as our collaborating geologists suggest that it provides more insight into the attributes when the seismic amplitude volume is co-rendered as a context. 4.6.2 HDTF Editor The user can interact with the HDTF editor to manually modify the HDTFs. Figure 4.6 shows the HDTF editor where the PCP axes reorder button and attribute- 55 wise control panel can be seen on the top, the PCP TF editor is seen in the upper left, a group manager is shown to its right and the pairwise TF editor is shown in the lower part. The attribute-wise control buttons allow the user to specify a color map, toggle sampling between linear and nearest neighbor, and toggle lock/unlock for each attribute. A locked attribute is essentially an attribute with its entire value range used in TFs. In other words, it can be visualized in the panel view but is not contributing to classification. This setting is useful since not all attributes provide positive assistance in the extraction of specific features and this knowledge is usually not known beforehand. Als