Field of invention

The present invention generally relates to a method and a system for carrying out visual queries with respect to volumetric data which are to be visualized together with a data of geometric objects, wherein data are stored in a database and

computer means comprising visualization means are provided for visualization of structures and data.

The invention is in particular applicable with respect to the visualization of neurobiological structures, or neuronal data, respectively, or comparable large amounts of (volumetric) data .

Background of invention

A major goal of neuroscientists ' work is to explain specific behavior of living beings, especially humans. However, human behavioral traits are complex and difficult to comprehend. For this purpose, researchers explore the anatomy and morphology of neuronal circuits of simpler species to identify their meaning and functionality. The fruit fly Drosophila melanogaster is a favorite organism in neurobiology research because it facilitates studies of complex systems on a relatively simple model. For this purpose, large databases of neuronal structures acquired by microscopy scans were built and adapted for computer- aided exploration and visualization, and visualization techniques tailored for exploration of biological structures are used .

However, orientation in large collections of spatial structures still poses a problem. The table-view database interfaces in use, although allowing filtering of items and accessing known subsets of data, do not support guide selection based on spatial relationships. Furthermore, it is not possible at present to calculate and provide specific objects of interest in a short time, bearing in mind that such calculations should be performed as runtime calculations; due to the large amount of data, these runtime calculations with a view to quick visualization are not possible at present.

Summary of invention

It is therefore an object of the invention to provide a method and a system for visual queries which facilitate visual exploration of a large collection of e.g. neuroanatomical structures, and allow rapid calculation of spatial data of structures of interest, for visual exploration and queries.

In particular, it is an object to give access to structures located in the proximity of a sketched path, and this in comparatively, short response times.

Furthermore, it is an object to select objects based on their mutual spatial distances.

A further object is to allow fast browsing using semantic relationships stored in the database.

Thus, it is intended to provide rapid runtime calculation to obtain quick visualizations, with response time e.g. in the ms- range .

To solve the given problem, the invention provides a method and a system as defined in the independent claims. Preferred embodiments are subject of the dependent claims.

According to the present invention, a specific data retrieval technique is provided to be able to quickly access and read the wished data from the database associated with a server, this in spite of the large amount of data stored in the database. Accordingly, a technique is proposed to query a - large - database in a visual manner, by direct interaction, with a 2D (slice) or 3D visualization. In connection therewith, in a preprocessing step, for a row of volume elements of a region of interest, a respective distance table is calculated; then, on the basis of said distance table, an associated look-up table is derived which can be handled in a memory with limited capacity, in a runtime memory, and which nevertheless allow rapid access to data by looking into this look-up tables and, thereafter, into the corresponding distance table. The look-up table contains a count of pairs of structure distance values and structure IDs as contained in the distance tables, as well as an offset of the row in the distance table; according to this, the look-up table contains less amount of data which, nevertheless, allow to read out respective pairs in the distance table, so that it is then possible to finally read out the structure data from the database, this all in a very short time.

Preferably, the calculation of the distance table is performed on the basis of a 3D Hilbert scan, to obtain neighbourhood-preserving storage patterns of rows. It should be mentioned here that Hilbert scans, whether in 2D or 3D, are known per se.

Preferably, the distances between the structures and the respective volume elements are limited to a maximum distance value, for instance to 40 voxels. That is, when reading out data from the database, only structures are taken in consideration which have a negative distance, i.e. are intersecting, or are within this maximum distance.

As already mentioned, data of neurobiological structures are very voluminous, and it is useful to store such voluminous data in a central database, connected to a server. Furthermore, it is preferred that the distance tables, too, are stored in an off- core memory.

Normally, after an amount of data, in particular neuronal data, have been retrieved, further data will be retrieved later on, and it is suitable to add this new data to the database already existing. Accordingly, it is advantageous to pairwise merge distance tables and looking-up tables referring to data retrieved at different times to a common distance table and to a common look-up table.

Preferably, also object queries are carried out on the basis of special distribution of objects around a selected object of interest. Here, preferably, each object query provides a list of objects close to, or intersecting the selected object as well as the minimal distances in voxels.

Furthermore, it is possible by semantic queries to retrieve semantic information associated with objects in the database and to display this information together with the object query result .

Brief description of the drawings :

The invention is described in the following by referring to preferred examples to which the invention is not limited, as well as with reference to the drawings. In these drawings,

Fig. 1 schematically shows a central neuronal system of a fruit fly;

Fig. 2 illustrates a template volume in conjunction with a Gal4-volume in the neuronal cluster of a fruit fly;

Fig. 3 schematically illustrates a 2D-Hilbert curve s(t), with a map of a 2D space, shown in part (a) , which is mapped on a ID curve shown in part (b) ;

Fig. 4 shows several examples of curves for mapping a 2D space onto a ID curve; namely so called "Peano curves", in line (1), in parts (a), (b) , (c) ; "Z-order curves" in parts (d) , (e) ,

(f) (line 2); and "Hilbert curves" in line (3), that is in parts

(g) , (h) and (i) of fig. 4; in particular, here the curves are shown, in the principle form, then in a first iteration and finally in a second iteration;

Fig. 5 schematically shows a system for carrying out visual queries according to the present invention;

Fig. 5A shows a screen shot on the visualization display of such a system, illustrating a list of structures in the selection panel and further semantic information in connection with the visualization of a fly brain;

Fig. 6 schematically illustrates, in parts (a) and (b) , how distance tables may be derived from data objects, here in a simplified 2D-version;

Fig. 7 schematically shows the sketching of an interface- path comprised of a row of points, for carrying out path queries in the neuronal data of a fruit fly;

Fig. 8 shows a scheme illustrating the principle concept of a point query;

Fig. 9 illustrates, in upper line (A), the calculation of a distance table and a look-up table, here again in simplified 2D manner, and in lower line (B) the use of the tables for a point query according a path query as outlined in Fig. 7;

Fig. 10 shows in a manner similar to Fig. 9 the calculation of a distance table and a look-up table using a 2D Hilbert scan (first line A) and the use of the tables for a point query (line B) ;

Fig. 11 illustrates a scheme for a preprocessing pipeline;

Fig. 12 schematically illustrates 3D Hilbert scans of different size (2x2x2 voxels in part (a) , 4x4x4 voxels in part (b) and 8x8x8 voxels in part (c) ) ;

Fig. 13 schematically shows the indirections in the case of 2D Hilbert scanning;

Fig. 14 shows, in parts (a) and (b) , schemes representing the merging of distance tables and look-up tables, in the form of flow charts, where it is also shown that an intermediate look-up table is stored since in the future it could be needed to merge a given distance table with other distance tables;

Fig. 15 and 16 show schematic illustrations of examples of merging using Hilbert indirection (fig. 15) and z-order indirections (Fig. 16), where element-wise addition is used to concatenate look-up tables and a union operation is used to

concatenate the distance tables;

Fig. 17 illustrates a scheme showing a concept of an object query; and

Fig. 18 schematically illustrates the investigation of voxel pairs for two objects i and j.

Detailed description of the drawings

Orientation in large collections of spatial structures still poses a problem. Traditional table-view database interfaces allow filtering items and accessing known subsets of data, but do not support selection based on spatial relationships.

In the following, a technique of visual queries, for rapidly accessing structures from the database based on the spatial information, will be described.

To give an introduction, it is first referred to Fig. 1 which schematically illustrates a part of a fruit fly 1, showing the brain 2 and the ventral nerve cord (VNC) 3 of the fruit fly 1. These two organs 2, 3 are composed by neuronal cells called neurons. Neurons transmit excitements by electrochemical signaling. The communication between individual neurons happens via connections called synapses. The short processes called dendrites are excitement-receptive and lead the signal into the cell body. The long process called axon transmits the excitement from the cell body to other cells. The neuronal tissue is supported and nourished by tissues called neuroglia. The network of dendrites, axon and glial branchings, which form the bulk matter called neuropil, embeds the neuronal cell bodies.

This rather simple organism of the fruit fly is a favorite model for research because it facilitates studies of intricate systems on a simple example; studies of the simple organism of the fruit fly can help to demystify relevant complex human traits such as sleep, alcohol dependency, and neurodegenerative dis- eases such as Parkinson and Alzheimer.

The Gal4/UAS system is used to label and manipulate specific neurons in the fly brain and ventral nerve cord. The brain and nerve cord are separately dissected. Specific neurons are stained. Additionally, separate neuropil staining is used to facilitate registration - it highlights regions of high synaptic density which provide a .stable morphological reference. After preparation and staining, the tissues are scanned using a laser scanning confocal microscope to generate 3D images (i.e. stacks of 2D images) of the respective object.

The image showing the neuropil staining is then used to perform non-rigid registration of the scans to a corresponding template for either brain or ventral nerve cord. The template was generated by averaging a representative set of scans registered against a reference scan. The registration process itself is automatic, but results are manually verified and additional image processing operations may be applied to reduce noise. Only scans considered to be registered with sufficient accuracy are used for the database.

The following table 1 shows typical types of neuronal data:

Table 1: Types of Neuronal Data

Each neuron is characterized by three types of features: a cell body, the axon tract, which is an elongated structure that spreads over large areas, and arborizations, which contain synapses where communication with other neurons occurs. Neurons are classified based on the morphology or shape of these features.

Overall, four different objects are segmented: neuropils, cell bodies, axon tracts and arborizations. The segmented surfaces of neuropils, cell bodies and arborizations are stored as triangular meshes. In contrast, axon tracts are represented as skeleton graphs.

All acquired volumes and objects, i.e. the corresponding data, are stored on a file server and their references are inserted into the relational database. The database also facilitates storage of annotations of individual objects and allows the user to associate different objects into neural clusters.

Fig. 2 depicts an example of a template 4 of the brain of a fruit fly 1 which provides anatomical context in conjunction with a Gal4-volume 5 and a neural cluster which consists of three arborizations 6, two axon tracts 7 and one cell body 8.

There are on-line atlases available which offer illustrated and annotated descriptions of organs. They help to localize neuronal structures and bring a system and coherency into the nomenclature of anatomical structures. Coherency, a logical system, and localization are crucial in computer-aided applications. In neuroscience, the construction of such tools is helpful for easy and quick localization, search and analysis of neurological structures. In the article aye et al., "Visualization, reconstruction, and integration of neuronal structures in digital brain atlases", International Journal of Neuroscience, 116 (4 ): 31-459, 2006; an overview of methods for neurological data reconstruction and visualization, creating atlases and registering the data is presented. A survey on multidimensional and metric data structures and query methods for spatial databases is given in Gaede et al., "Multidimensional access methods", ACM Computing Surveys, 30 (2 ): 170-231, 1998. Spatial data are defined by an explicit information about their extent and position in space. Gaede et al. discuss specific types of queries on these types of data - point and region queries. A point query returns all objects intersecting a selected point in space. A region query selects all objects intersecting a selected region in space. Users can interact with the query results, e.g., make a union or intersection.

There exists no linear mapping of a multidimensional space into ID space which fully preserves spatial locality. Voxels which are direct neighbors in a 3D volume might not be direct neighbors when serialized on a disk. However, there exist techniques which attempt to maximize locality, as would be important for fast search operations in spatial databases and to ensure database scalability. Among numerous well-known approaches to the organization of spatial data on a disk, Kd-trees, Octrees, BSP-trees and space-filling curves may be mentioned as prominent ones, cf. e.g. Donald Hearn and M. Pauline Baker, Computer

Graphics with OpenGL. Prentice Hall, 3rd edition, 2004, with respect to tree configurations. Space-filling curves s(t) are fractal curves which provide one-to-one mapping between points in a multidimensional and a ID space R ⁿ → R and well-preserve the spatial distance between points. A family of 2D space-filling curves was introduced by Giuseppe Peano in 1890.

Fig.3 illustrates a space-filling Hilbert curve s (t) which maps a 2D square into a ID curve R ² → R.

A discrete representation of these curves can be computed using an L-System, cf. Przemyslaw Prusinkiewicz and Aristid Lindenmayer, "The Algorithmic Beauty of Plants", Springer-Ver- lag, 1st edition, 1990. An L-System consists of an initial shape, alphabet and grammar building rules. Iterative application of grammar rules transforms the initial shape and converges to the solution.

In Fig. 4, building rules are shown for different types of space-filling curves - a Peano curve - first line (1); a Z-order curve - second line (2); and a Hilbert curve - third line (3).. After an infinite number of iterations on a 2D curve each point of a square belongs to the curve. Space-fillings curves are used to order discrete multidimensional data on a disk. Fig. 4 shows a serialization, i.e., a linear order, of a 2 ^{b χ} 2 ^b matrix of points. In the following, it is refered to linear ordering of discrete multidimensional data as to a "scan".

For linear ordering of multidimensional data on a disk, a Hilbert curve has proven to be a good choice because in the family of space-filling curves, it preserves the proximity the best. In W. J. Gilbert et al., "A cube-filling hilbert curve", The Mathematical Intelligencer, 6(3) :78, 1984, the construction of a cube-filling Hilbert curve is described. There exist many algorithms for 3D Hilbert-scanning-algorithms .

Hilbert' s curve provides a scanning order for multidimensional data which well preserves locality. Thus, it is widely used in digital image processing. However, sampling of spacefilling curves s(t) poses a difficulty. Given a cuboid region and a point s(t) 6 R ⁿ , the deduction of t is rather cumbersome. At present, there exists no solution for direct computation of the curve s(t) or inversely s ^-1 (s(t)). The whole curve s(t) must be computed for the whole region using one of the proposed algorithms .

The exploration of neuronal databases has been investigated by several researchers. So, direct queries on raster data such as visual queries on images have been provided. For instance, A. Sherbondy et al., "Exploring connectivity of the brain's white matter with dynamic queries", IEEE Transactions on Visualization and Computer Graphics, 11 ( 4 ): 419-430 , 2005, describe a set of interaction techniques that makes the exploration and interpretation of pathways in brain easy. The queries are interactive - the user marks a box or ellipse-shaped region and all pathways crossing this region are shown. The framework supports boolean operations on the regions of interest.

Visualization techniques for microscopy data are being explored since computer-aided microscopy is used in bioscience. Early work on microscopy data focused on image processing techniques, notably noise reduction, because the early microscopy scans were of worse quality than today. The microscopy scanning devices recently progressed, and the microscopy images of today are of much higher quality. The evolution of computer and microscopy hardware solved many of the problems addressed in the earliest works.

More recent approaches employ modern and innovative user interfaces and state-of-the art visualization techniques. For example, 0' Conor et al . , "3D visualisation of confocal

flourescence microscopy data", Proceedings of EUROGRAPHICS, pages 49-54, 2004, and Wang et al., "Gpu-based visualization techniques for 3D microscopic imaging data", Proceedings SPIE, page 64981H, 2007, used programmable graphics hardware to achieve higher performance.

In numerous cases, researchers want to analyze and process a collection of datasets simultaneously, e.g., to test their hypo ^¬ thesis. Processing each dataset individually is time-consuming. On the other side, in fully-automated batch processing of the whole collection of datasets, only final results are displayed, and the user cannot intervene into the process. Leeuw et al., "Visualization and analysis of large data collections: A case study applied to confocal microscopy data", IEEE Transactions on Visualization and Computer Graphics, 12 ( 5) : 1251-1258 , 2006, describe techniques for visualization and analysis of large data collections applied on confocal microscopy data. The presented system uses a hybrid approach combining the two processing approaches. Furthermore, a plethora of visualization tools such as histograms, scatter plots and parallel coordinates is incorporated.

In the past, direct volume rendering techniques were combined with shaded isosurfaces using techniques for mixing volumetric images and geometry, and even rendering of multichannel . volumes. The system described herein uses known visualization techniques in conjunction with a novel approach of quick data retrieval in a database of neurobiology structures, using also spatial and semantic information.

For the present method and system, it is preferred to use the open-source volume rendering framework so called

"VolumeShop", cf. Bruckner et al., "An interactive system for direct volume illustration", Proceedings of IEEE Visualization, pages 671-678, 2005, where specialized plugins as described hereinafter may be smoothly integrated. "VolumeShop" has a flexible architecture which allows to rapidly integrate new functionality. The architecture of the system consists of a

relational database connected to a file server, spatial indices of the data, and the visualization tools. The entities communicate via interfaces. The latter are responsible for the construction of the visual queries.

Fig. 5 illustrates the connections of entities and concept of the entire system 9 including - preferably local - computer means 10 (CPU) as well as - preferably remote - server means 11.

In this system 9 all data provided are stored on a file server 11. The application also supports offline access where the files are stored on a local medium and a local relational database 12. However, the database 12 and the server 11 are intended for access over a network, and the datasets are transferred to the local machines (e.g. 10) on demand. In this manner, researchers/users are able to provide and maintain their data on the server 11 and their clients view and explore them at local machines.

All data provided may be present in a native general-purpose file format. As shown in Table 1 above, templates, confocal images, and average-images are stored as volumes (i.e. 3D images) . Arborizations, neuropils, and cell bodies are represented as geometry and projections as skeleton graphs.

The relational database 11 may store additional information about the present neuronal structures. This additional informa- tion may comprise properties such as expert's comments, clusters, whether the structure was retrieved from a brain 2 or from a VNC 3, and the gender of the organism it was retrieved from. A user interface 13 contains a selection panel displayed on a visualization means (display) 14 which shows a list of st ctures in the database 12. A user 15 selects individual structures from the file list to be added into the visualization, which then also appear in the view panel. The file list in the selection panel on display 14 can be filtered using additional information stored in the database 12. In the search panel, the user 15 can make a set-up to display only the structures which were retrieved from a female brain. The structures which were retrieved from a VNC or from a male organism do not appear in the list then. The search panel may also contain a field for keyword-based search.

Fig. 5A depicts a typical screenshot of the user interface showing a visualization window, a view panel, a search panel, and a selection panel. In the search panel, a restricted list with the keywords "projection" and "only" is present. The restricted list in the selection panel contains only files containing both keywords in their annotations. This feature easies the orientation in the long file listings stored in the database .

In Fig. 5, then, data interpretation and data handling is shown at block 16, with a runtime memory 16' shown schematically, and a visual query module is schematically shown at block 17, that is within the (local) computer means 10. An arrow 18 illustrates direct queries in the database 12; finally, at 19, preprocessing of data is schematically illustrated, to obtain spatial indices, as shown at 20, and as described more in detail below .

The database 11 provides context-based search. The user 15 picks structures by their names or by using filtering options. However, up to now, search and exploration based on the spatial information is not supported; users cannot automatically select structures crossing particular points in the volume, the closest objects to an object of interest, and directly retrieve information about objects by clicking on them in the visualization. However, these tasks are of great interest of researchers. To this end, specific queries, here called "path-based queries", are proposed. A path-based query (or in short: path query) provides efficient retrieval of structures crossing a tight region around a path sketched over the visualization. The concept of a path-query is that, first, the user interactively sketches a path over the visualization as in Fig. 5A, cf. also Fig, 7. Consequently, a list of available structure crossing this region appears in the menu. Second, the user loads, e.g., all axon tracts from the list path.

It is preferred to provide the user with two additional types of visual queries: object queries and semantic queries.

Spatial indices, as shown in Fig. 5 at 20, represent packed spatial information about individual objects stored on the file server. They are advantageous for efficient path and object queries. At runtime, it is wished to rapidly calculate distances to n closest structures from any point in the volume. Furthermore, it is wished to quickly provide a list of n closest objects to any object of interest. Considering that the file server 11 stores several 10.000 or even 100.000 of structures, a runtime calculation normally would be impossible. Thus, two tables containing point and object distances are computed in the preprocessing step 14, using a separate tool.

Fig. 6 schematically illustrates fields of the mentioned distance tables. In part (a), each row of a distance table 21 refers to a point 22 in the space. The distances are measured using the chamfer metric. In Fig. 6, part (b) , each row 21' refers to an object Al, CI, C2. When two objects (e.g. CI, C2) intersect, the corresponding value in the table 21' is negative and its absolute value indicates the intersection volume. The application accesses these precomputed tables 21, 21' for distance look-up at runtime which makes the distance calculation fast .

In the following, the concept of visual queries - path, object and semantic queries, as implemented in the present system - is described more in detail, together with how a user may define each type of query. As mentioned, path (and object) queries require precalculation of spatial indices 20.

The present path queries are a powerful searching method; they enable searches for structures crossing a defined region. These regions are user-defined using an intuitive sketching interface, compare Fig. 7. Defining a region of interest while sketching an arbitrary path 25 is more flexible than selecting regions with some rectangular or elliptical markers. The query provides all structures crossing an area tight to the user- defined path 25. This area is suitably restricted to a maximum radius, e.g. a radius of 40 voxels around the path 25. In practice, each path query 25 falls apart into a set of point queries at points 26, as illustrated in Fig. 7.

Using the described data structures, it is possible to efficiently determine objects in the vicinity of a voxel. For performing path queries, it is therefore necessary to identify a corresponding 3D object-space position for each 2D point along the path. Whenever a new point is added to the path, a depth buffer is read at the corresponding 2D location, and the inverse viewing transformation is used to transform it into object space .

For geometry rendering, this results in the position on the surface of the object closest to the viewer.

For volume rendering, however, several samples along a ray may contribute to a pixel. For each viewing ray, it is therefore preferred to write the depth of the sample which contributes most to the final pixel color. Particularly in conjunction with the volume rendering technique this approach has proven useful - as stained tissues are presented visually more prominent, the respective depth will give access to objects in the vicinity of high-intensity volumetric structures. Using the depth buffer in this way also ensures that there is always a good correspondence between the selected query locations and the actual visualization when operations such as cropping have been applied or a slicing plane is displayed in the 3D visualization.

More in detail, individual point queries work as follows. For each point 26 (Fig. 7), the (x, y) -coordinates in the screen space are known, and the z-coordinate is retrieved from the depth buffer. The point 26 is transformed from screen space coordinates P _s into volume space coordinates p _v . As the volume space coordinates p _v are rounded, p _v refers also to a voxel in the volume. With p _v as input, the point/path query delivers a list of pairs < d, id > with a distance d to a respective structure and a structure identifier id, sorted by d. The distance d refers to the shortest distance from p _v to the surface of the structure id. When p _v lies inside of a structure, d is negative, as already indicated. The distance restriction of e.g. 40 voxels does not apply for negative distances because there is a great interest in every object if p _v lies inside. Fig. 8 shows the general concept of a point query 17'.

In the preprocessing step 19 (s. Fig. 5), lists for each point in the volume, i.e., distance tables 27, are prepared. The size of the distance table 27 file grows with the number of objects on the server 11. As it shall be possible to insert new files, a significant growth of the distance table file size is to be expected. For these reasons, this file is kept in an off- core memory and not in the runtime memory (RAM) . A precomputed look-up table (LUT) is used then to efficiently retrieve rows of the distance table 27 stored in the off-core memory. For each point in the volume p _v , the LUT provides an offset of the row corresponding to p _v in the distance table 27. Second, the LUT provides a count of the pairs < d, id > in each row. The size of the LUT is constant and independent from the number of structures, so it can be kept in the runtime memory.

In Fig. 9, the concept of the distance table 27 and LUT table and their use for a query is schematically demonstrated on a 2D example with three structures Al, CI, C2. With respect to a 3D image (volume), a 3D curve, e.g. a 3D Hilbert curve, is scanned; for each voxel reached, the neighbouring objects and the respective distances are calculated and stored in the distance table in the corresponding line as pair (distance; object id). The distance table is established, i.e. filled, in the order of the processing of the voxels,; therefore, the line index i of the matrix corresponds to the processing order of the voxel Vi of the 3 D image. In the case of a Hilbert curve, the positive neighbourhood properties are transferred thereby to the distance table; in this case, a temporary look-up table is generated, as will be described herein after with reference to Fig. 10.

In the example of Fig. 9, the distances d were restricted to maximum of d _max = 3 voxels. In Fig. 9, part (A), every point

(0,0), (1,0) ... (3,3) of a 4x4 region is investigated, the sorted list of pairs < d, id > are determined (distance table 27), and the LUT is created on the basis of the distance table 27.

In Fig. 9, part (B) , the user selects a point, e.g. (1, 2), in the region; on the basis thereof, the offset (3) and the count (24) are retrieved from the LUT. This is the key informa- tion to read the row corresponding to the list of pairs < d, id > in the distance table 27 which was created in the preprocessing step. In the example shown, the disk read head skips offset = 24 pairs of the file and reads subsequent count = 3 pairs.

The distance table 21 and the LUT are serialized, i.e., by their ascending linear indirections lin(p) as defined in the following Equation (1) with x x y x z = volumesize. Under "indirections", an ordering pattern of voxels is understood.

While drawing a path, point queries for points located close to each other are automatically called. The' point queries fetch the respective rows from the distance table 27.

To improve the performance, it is preferred to use neighborhood-preserving storage pattern of rows using a 3D Hilbert scan. Thus, the rows of the distance table 27' using Hilbert scan indirections are stored, as schematically demonstrated in Fig. 10. The < d, id >-pairs inside the rows of the distance table 27' are not permuted. From this distance table 27', a first look-up table LUT' is generated, as described above; here, line no. i (with i = 1, 2, 3,...) of the distance table 27' corresponds to line i in the LUT'. The two columns of LUT ' contain the following information: number of all objects which were neighbours of the voxels which have been processed before the actual voxel v _A (1 ^st column) ; and number of the objects which are neighbours of the actual voxel Vi (2 ^nd column) .

To create the final LUT, the lines of LUT' are permuted such that their order corresponds to the "natural" linear order of the corresponding voxel in the volume (3D image), cf. equ. (1), to enable quick - temporary - access to the information in the LUT. Accordingly, by reordering the (count and offset ) -pairs in accordance with the ordering (0,0), (1,0), (2,0) ... (3,3) of the 2D region as shown in Fig. 10, the final LUT is obtained. The final LUT must be accessible using linear indirections. It thus contains the offset and count information with relation to the permuted distance table. Comparing Figures 9 and 10, it can be seen that the fetch of a row from the distance table 27 or 27' has the same logic if the permutation is used when using 3D Hilbert scan indirections. Using the described storage pattern en- ables to keep the easy and' quick look-up in the LUT and to gain performance while reading from the distance table file. With other words, the information in the LUT makes it possible for each voxel of the image or volume to directly determine at which position of the memory containing the distance table the information with respect to the neighboured objects may be found, and to calculate how many objects are "neighboured", i.e. until where the memory must be read.

The results of all point queries along the sketched path are merged together and presented to the user in an in-window hypertext label. If there are more < d, id > pairs for one object id, only the pair with lowest distance d is kept.

For preprocessing, a preprocessing pipeline 30 is used which consists of six stages, as shown in Fig. 11.

The structures in the database 12 may be represented as geometry or skeleton graphs. According to Fig. 11, in a first step, distance field volumes 31 are gradually created for each structure. For skeleton graphs, it is possible to extrude the graph with some radius and to create the surface geometry representation. For geometry, it is possible to compute signed distance fields, to obtain distance tables 27/27', s. step 32 in Fig. 11, thereby choosing signed distance fields using chamfer metric. Each voxel of the distance field volume stores the distance to the surface in voxels. Signed distance fields enable to distinguish between the inside and the outside regions of each structure. The negative distances represent the inside area. It is possible to record the mentioned steps and automate the processing using batch files. This makes it possible to generate a distance field for each structure without individual intervention .

According to step 33 in Fig. 11, a 3D Hilbert scan for indirection indices is generated. In particular, a cuboid region with the size of the volume xxyxz is scanned. Fig. 12, parts (a) , (b) and (c) , schematically illustrate 3D Hilbert scans of differently sized regions.

In Table 2, the indirection indices corresponding to the scan shown in Fig. 12 are given.

Table 2: Indirections for 3D Hilbert scan of size 2x2x2.

In accordance with the above, a distance table with n = xxy x z empty rows may be prepared, and the distance field 3D Hilbert scan s(t) Volume coordinates p _v linear indirection lin(p _v )

volumes are processed voxel-by-voxel. For each voxel v with position v= ( JC, v _v , v _z ) ^T _r the distance d which is the voxel value is investigated. If d ≤ 40, the pair < d, id > is inserted into the r-th row of the distance table. In this case, simple linear indirection of the volume is used where r = lin{v) . The "id" is the identifier of the currently processed distance field volume

(structure) .

After the last distance field volume has been processed, the rows of the distance table 27/27' are sorted ascendingly with respect to the distance d of the pairs < d, id >. The entire rows of the distance table are written on the disk obeying 3D Hilbert scan indirections. For example, the rows of a distance table of a 2x2x2 volume would be stored in a permuted order (0, 2, 6, 4, 5, 7, 3, 1) as indicated in Table 2 above.

Then, the look-up table LUT' is generated (step 34 in Fig. 11) corresponding to the permuted distance table. Consequently, the LUT' is also permuted and cannot be accessed via simple linear indirections. However, the permuted LUT' can be needed for later merging of tables. Thus, this table LUT' is written as temporary data.

Thereafter, the LUT is generated, with the line permutation, as described above, cf. also steps 35, 36 in Fig. 11. The per ^¬ mutation (inverse indirections) map the 3D Hilbert scan back to the linear indirections.

In Fig. 13, the mapping from linear indirections to 2D Hilbert indirections and the mapping from 2D Hilbert indirections back to the linear indirections are illustrated. Both mappings are represented as permutations.

After the last stage of the preprocessing procedure, two tables, LUT and the distance table 27 or 27', have been generated which are ready to be used by the application. It is to be kept in mind that preferably, also the temporary permuted LUT' is kept which may be used for merging. Namely, expecting a significant growth of the database, the present system provides merging distance tables and LUTs. Furthermore, the preprocessing calculation requires a lot of RAM. However, the amount of available RAM is limited, so it cannot be guaranteed that all data can be successfully processed at once.

It is possible to overcome the RAM limitation by merging as follows: The data are split into n disjoint groups, n distance tables DT _a , DT _b ... and LUTs LUT _a , LUT _b ... are generated, and these are merged pairwise. Finally, a single distance table DT _ab and a single LUT table LUT _ab representing all structures are provided. The major requirement for the merging algorithm are handling of large data and requiring a low amount of RAM. Fig. 14 illustrates, in a flowchart, such merging of two distance tables DT _a and DT (part A) and two look-up tables LUT _a and LUT _b (part B) . Rows of the distance tables DT _a and DT _b can be concatenated using a union operation because they refer to two disjoint set of objects. The number of pairs in the new row equals the sum of the pairs in the corresponding rows in DT _a and DT _b . This property makes merging of look-up tables (LUT' and LUT) easy. Offsets and counts of two old look-up tables are simply added, and a new valid look-up table is obtained.

In Figures 15 and 16, there are given two simplified examples of pairwise merging of tables using Hilbert, or Z, indirection. It is shown that addition of the look-up tables LUT and union of distance tables 27 works also for different indirections .

The system described, with its specific data preprocessing and corresponding use of storage means, has several essential advantages. In the present approach, the precomputed distance tables 27 are used to determine closest objects for each point which has the following particular advantages. The table is precomputed only once in a preprocessing step and the tables are merged when new data are to be inserted. At runtime, however, only one look-up per point of the sketched path is performed. Moreover, the distance table stores the < d, id >-pairs in a sorted order, so that it is possible to perform no sorting and no distance calculation at runtime. The knowledge of the point distance d is essential as the key for sorting and for approxim ation of the real distance. The performance of the application also benefits from the locality-preserving storage pattern applied to the distance table using the 3D Hilbert scan indirections .

The present approach has been compared with the so-called Fibonacci coding of the spatial data which bases on the

Zeckendorf theorem. This theorem states that every positive integer n can be uniquely represented as a sum of nonconsecutive Fibonacci numbers F _k (k = 0 ... L) .

At the present system, every object in the database 11 is assigned an id = F _k . When objects i and j intersect, and idi = F then id _j ≠ F _k+i . The voxels store n which is the sum of the ids o all incident objects. This refrains from storing the whole list of ids for each voxel. However, the Fibonacci numbers increase rapidly - F ₈ exceeds the range of the 32-bit unsigned integer. As a number of objects to index L > 300 may be foreseen, this method becomes inconceivable regarding storage. Furthermore, this method requires decomposition of the indices into F _k s and calculation of the distances at runtime which is a disadvantage

In comparison to the present method, prior art approaches such as space-partitioning, e.g., kd-trees or octrees have the following disadvantages: Many neuronal structures in the database 12 have intersections. Space-partitioning methods would represent them as single tree-nodes. Consequently, for each node, saving a list of id-pairs of all incident objects is fore seen. In contrast, the distance table stores the lists of < d, id >-pairs for every point in the volume. The storage requirements of a kd-tree or of an octree are inferior to the storage requirements of the present distance table because the tree- nodes cover more than one voxel. If the tree-nodes store < d, i >-pairs, the storage requirements are equal to the present approach, so the space-partitioning method brings no advantage. Assuming that the tree-nodes do not implicitly store the distance information to the incident objects, the application must determine them at runtime while performing computational costly operations. Therefore, the performance of the space-partitionin methods is inferior to the performance of the present method. A one of the goals of this invention is to implement an interact- ive tool, the performance is the primary benchmark and storage the secondary. Furthermore, the implementation of the present method is intuitive and simpler as the implementation of a space-partitioning method. Also, merging of two space-partitioning trees is more complex and costly then merging of two distance tables as proposed in this section.

Hypertext labels may be used to present query results to the user on display 14. They pop up in-window which enables the user to keep focus on the visualization. A label appears in the corner of the window when user clicks with the mouse button and starts sketching the path. A color-coded slider may be attached to the label which enables control over the search radius around the path. The radius is by default set to 0 voxels. The label displays hyperlinks to whole subgroups, and to individual structures within each subgroup, and also the number of structures found for each subgroup. The hyperlinks use the same color-coding scheme as the slider which enables the user to quickly associate the structure some degree of proximity. This information is updated while the user keeps sketching.

When the user releases the mouse button, the path is finished, and the label moves to the center of the window. The user can adjust the search radius and inspect the results. Eventually, the user may click on hyperlinks to load structures individually or a whole subgroup of structures at once.

In addition to path queries related to spatial distribution of objects around points in the volume, the preferred set includes also object queries. Other than to path queries, they are related to spatial distribution of objects around a selected object of interest. In practice, this query provides a list of closest objects and their distances in voxels. If two objects intersect, the volume of their intersection is provided instead of the distance.

Object queries 17' (cf. Fig. 17) work as follows. It is preferable to provide an "id" of the selected objects as input and the query returns a list 37 of objects sorted by the distance. The negative distances again indicate intersection volumes. Similarly as for the path queries, a table, i.e., object table, is created in a preprocessing step.

In practice, it may be started with an η ^χ η-matrix . Μ ₁₃ stores the minimal distance between objects i and j, or eventu- ally their intersection volume. All matrix elements Mi _j are initialized to the highest distance which can occur, i.e., d _max = , with xyz = volumesize. Each pair of distance field volumes described above is inspected voxel-by-voxel. The matrix M _. is symmetric, so once the pair of distance fields ij has been inspected, it. is not necessary to inspect the pair ji.

Voxel pairs and vj , i.e., voxels of the object i and j at the position p in the volume, are gradually investigated, and the following cases (see also Fig. 18), are distinguished:

The value of v = 0: here, v belongs to the surface of the object i; the matrix M is modified as follows. If Vj > 0 and < Mi _j , the value of Mi is replaced by . For example, the sampling position A in Fig. 18 (a) fulfils this condition - it lies on the surface of the object i, but outside of the object j. If ≤ 0 then p belongs to the intersection volume. For example, sampling positions B and C in Fig. 18 (a) comply with this condition, because they both belong to the surface of the object i and lie either on the surface or inside of the object j . In this case, Mi _j is set to -1 if Mi ≥ 0, or M _Aj is decremented otherwise.

The value of ^v < 0, and < 0: both voxels belong to the intersection of the objects i and j, i.e., sampling position A in Fig. 18 (b) ; the matrix M is modified as follows. If M _±j ≥ 0 then Mi _j is set to -1, else M _1;j is decremented.

The value of > 0: here, v belongs neither to the surface of the object i, nor to its inside area, as it is the case of the sampling positions A and B in Fig. 18 (c) ; in this case, the matrix M is not modified.

Elements of M _u above the diagonal (i > j) are filled, and the symmetry is used to copy , to elements below the diagonal Mji . From each row of the matrix M , a list L consisting of < d, id >-pairs is generated where d = Mi _j and id = j . As there is no interest in elements M _Aj such that i = j, so we exclude them from L. In the final step, L is sorted by increasing distances d. The finalized object table stores minimal distances between objects, or the total volume of their intersection. In the case that objects have multiple disconnected intersections, the sum of the volumes of their partial intersections is stored.

An object table needs to be updated only when new objects are inserted into the database. Let A be the set of old objects and B the set of inserted objects. Then ADS = ø , |A| = n, and | B| = m. It may be proceeded as follows.

Initialize an (n + m) ^χ (n + m) -matrix M and use the old object table to fill the maximum of elements.

For all objects £ B, compute a m * m-matrix M' of distances and intersection volumes. M' i _j stores a distance or an intersection volume between the objects i and j, such that i, j e B. M is used tofill the additional elements of M.

For all objects e B, compute an m n-matrix M' ¹ . Matrix elements W ' i _j store distances or intersection volumes between objects i and j, such that i e B and j G A. M ¹ ' is used to fill the remaining elements of M; sorting is carried as described above .

Researchers sometimes annotate objects in the database and associate them to clusters. Semantic queries may retrieve this information and display it together with object query results in a hypertext label as mentioned above. Associated structures are listed as hyperlinks. Activating them triggers loading of the respective objects. This is a powerful method of browsing

through the data. Structures such as arborizations represent synapses which are often part of more clusters. Results of a semantic query allow the user to load them and explore the further connections. In this manner, neuroscientist can follow neural paths .

The visualization entity as indicated at 14 in Fig. 5 is an important node of the system. Visual queries are based on an interactive visualization. Numerous prior art visualization tools allow 3D data exploration which may include standards for viewpoint control such as rotation, zooming and panning, masking control such as clipping planes and cropping boxes and different renderers. The latter provide slice-rendering, direct volume rendering with interactive transfer function and geometry rendering. Other utilities supporting user's work comprise synchronized viewing windows, image and video-capture and saving of the current session. Such visualization techniques, need no detailed description as, in principle, they are state-of-the-art. The visualization entity 14, however, is essential to carry out a visual query.

One of the major requirements of researchers is the support for combined visualization of geometry and volumes and concurrent rendering of multiple volumes. A typical case is, e.g., locating more structures such as cell bodies and arborizations inside a template, which provides an anatomical context. At the same time, they may want to explore additional volumes. In this particular case, concurrent visualization of multiple volumes and geometry is necessary. Furthermore, exploration of numerous structures at the same time leads to a visual clutter. The aim is now to maximize the information the user can get from the visualization and to minimize the clutter. Thus, approaches to the concurrent visualization of volumes and to the geometry rendering are discussed in short.

Maximum Intensity Projection (MIP) is often used in biomedical visualization. This technique retains the highest values along the projection rays. Researchers also prefer MIP over direct volume rendering (DVR) . However, employing visualization using MIP destroys the depth perception. On the other hand, using DVR instead of MIP leads to occlusion of important maxima. This is a significant problem especially when a researcher explores several volumes and the template. The template should build an anatomical context but should not occlude the important information provided by volumes. Therefore, a hybrid approach called Maximum Intensity Difference Accumulation (MIDA) described in Bruckner et al., "Instant volume visualization using maximum intensity difference accumulation", Computer Graphics Forum,

28(3), 2009, may be used. This technique ensures that the strong features of the stained data remain unoccluded by the template volume .

The article of Bruckner et al. describes MIDA for mono-channel data, employing an additional modulation factor for the color and opacity composition equation used for DVR. Equations (2) present the color composition in DVR. The color and opacity are accumulated along the rays. £,· and ¾, ^■ denote the color and opacity of the i-th sample along the ray-traversal. C _L and A _L represent the accumulated color and opacity at the i-th step of the ray- traversal. The initial values C ₀ and A ₀ are both zero. A/ = A/_ i + ( l - A/_ i )Sl/, i = \ ..n

C _i = Ci- _l + { l - Ai- _l )% _i £i, i = l ..n (2)

Equations (3) denote a compositing method employed in MIDA. The additional parameter ¾/ represents the modulation factor at the i-th step of the ray traversal.

A/ = ¾ _i Ai- i + ( l - »/A/- i ).2l ₎ i = l ..n

d = ¾/Q_ i + ( 1 - ¾/A/_ i )¾,·£,·, = 1. i (3)

The modulation factors ¾, are calculated at each sampling step as in Equation (4) . The parameter 5 _L is introduced to make the later transition to multi-channel MIDA easier. Parameter denotes the density sampled at the i-th step and maxi (s _k ) denotes the maximal densit of samples taken at steps k=0..i - 1.

This composition method strengths maxima changing from a lower to higher value. This way, the maxima remain visible as in MIP. The other tissue along the ray also contributes, but with subtler transparency as in DVR.

It is wished to handle up to four different volumes. The volumetric data are scalar fields, thus it is possible to use four components (RGBA) to store four different volumes. In view of the given multi-channel data, the mono-channel approach of MIDA may be adjusted for multi-channel data as follows. For color and opacity mixing of each channel, the composition scheme described by Kniss et al., Gaussian transfer functions for multi-field volume visualization, Proceedings of IEEE Visualization, pages 497-504, 2003, may be used. Assume that there are n volumes Vj (x) , j = l..n, x e R ³ and n transfer functions assigning color and opacity to samples - Cj (Vj (x) ) and oij (V (x) ) , j = l..n, x E R ³ or simply c ₃ and οί ₃ . The opacity of each sample ¾/ may be calculated as the sum of opacities of each channel α-μ. The color of each sample C, equals the sum of all colors at each channel Cji weighted by opacities as denoted in Equation (5) : CjiCji

21 = ∑ α

7=1 (5)

Even though the compositing equation (Equation (3)) for multi-channel data is the same as for mono-channel data, the computation of parameters 23, is more complex. For each channel j and each sample along the ray i, δ-μ is calculated as in Equation (6) where Sji denotes the i-th sample of the channel j. The parameters 5 _j i are used for further calculations of the modulation factors 23,· in Equation (7) . The additional weighting of by opacity suppresses contribution of invisible samples. ji— maXj^Sjk), if sjj > max jj{s ) , k = 0../— 1

(6)

0, otherwise

23,· = l - max(Sji " ' ^' ) , j = 0..i - 1 (7)

Standard geometry rendering techniques may be used for all non-volumetric structures. Geometry rendering has many advantages. First, the standard per-pixel phong shading, is a realistic smooth surface representation. Second, this method is fast because it is hardware-accelerated. Third, surface representation with geometry allows easy rendering of outlines in slice views as mentioned above. However, e.g. the axon tracts are originally stored as skeleton graphs. To this end, a tessellation using cylinder-like extrusion of the graph may be provided.

Thus, the present visualization system is a rendering system based on a frame buffer/depth buffer. The 3D position p _v of a point in the resulting image can be determined by means of a transformation matrix as well as the depth buffer. Normally, for geometric objects, the depth as stored in the depth buffer cor ^¬ responds to that surface point which is the next to the viewer. However, this would not be apt for volumes, and therefore, the depth of that respective point along a viewing ray is written into the depth buffer which has the maximum contribution to the final pixel, as results from equations (3) above. The stored depth value corresponds to that point along the viewing ray for which (1— Q3i ,-_i) is a maximum.

As already mentioned, too, the respective query results are presented by means of interactive hypertext labels which are represented directly in the visualization window. With the aid of these labels, the query results may be inspected and loaded into, the display window. Furthermore, it is possible to interactively filter the query results (distance slider) .

To demonstrate the advantageous results achievable with the' present system, some use-case scenarios are presented now.

To illustrate a typical scenario of using a path-query, the user, first, explores an average Gal4-volume in conjunction with a template which provides anatomical context. The average Gal4- volume depicts a particular neuronal tissue of interest which includes, e.g., neurons and arborizations. An expert user is able to recognize structures within the average Gal4-volume.

Consequently, the user desires to further inspect these regions to exactly identify the structures. For example, the user

sketches a path over the region which was identified as a bundle of axon tracts. A hypertext label may show the axon tracts which cross this region, and the user loads them all by clicking on the axon tracts label. In the same fashion, the user marks a region and loads the closest cell body, and the closest arborization. Consequently, a proximity cloud indicating the extent of the search region appears.

As mentioned, object queries may provide information about distances and intersections between objects which enables fast and simple browsing. For example, the user is able to quickly determine which neurons cross a particular part of the brain. Regions of the brain which have been segmented from the template image of the brain can be loaded via the table-view interface to the database and selected with a mouse click on their visualization. The hypertext label displays the object query results. The user adjusts the slider so that only structures which have intersections are shown. Thus, only filtered structures appear in the filtered list. To load all filtered structures, the user may click on the hyperlink neuronal projections. Consequently, all axon tracts crossing the selected brain region appear in the visualization. Object queries also enable the user to quickly find a region of the brain to which a particular structure belongs to. To this end, the user may click on the depicted struc- ture in the visualization to prompt an object-query results label. The slider is adjusted to show only structures which intersect the selected structure. Subsequently, the user loads all tissue regions with intersections. In the view panel, the user may uncheck regions identified as not appropriate; then, for instance, only one tissue region may remain in the visualization for further inspection.

A typical case of applying a semantic-query is finding neural clusters which are related to a particular structure. In particular, arborizations are members of several neural clusters. Semantic-query enables the user to view and compare the entire clusters. Here the user may, at first, select a particular arborization. The semantic-query displays all cluster to which the selected arborization belongs to. Second, the user selects one of the proposed clusters, which is consequently loaded into the visualization .

In practice, all three types of queries may be used in conjunction with each other. For instance, the user desires to identify a particular region of the brain; then, the user uses the path-query to find the closest tissue region which includes the sketched path. Consequently, the tissue region displayed in the query results at the most left is loaded. By clicking on the loaded region, all information may be displayed which is stored in the database about this region and the results of the object query. When continuing browsing by loading all axon tracts crossing the selected tissue region, the highlighted axon tract is selected for further inspection. The user applies the semantic-query to load the associated neuronal cluster.

The performance of the present system has been evaluated as follows .

First, the latencies (in ms) while reading rows of the distance table from the disk were measured. To provide comparative values, the following criteria have been complied with. First, the latencies while reading rows of the distance table employed in practice (185MB) and while reading rows of a large distance table (2GB) were compared, s. Tables 3 and 4 below. For this purpose, an additional distance table with distance restriction of 256 voxels instead of 40 voxels was generated. Second, the latencies while reading rows of the distance tables which use the 3D Hilbert indirections as in Fig. 10 and the linear indir- ections as in Fig. 9 were compared.

Third, the latencies for sketching done on the 3D visualiz tion and on 2D slices were inspected. While the user sketches a path on a 2D slice, the path becomes a plane curve. However, while sketching on the 3D visualization, e.g., DVR, the path defines a space curve.

With respect to the results of the latency measurements as summarized in Tables 2 and 3, each comparative value presented in the tables is the average value over five independent measurements. During sketching, the system processes path queries and renders with frame rates up to 30Hz. All measurements were performed on a workstation equipped with an Intel ^® 6600 Core Duo CPU 2.4GHz, 2.0GB RAM, a hard disk VelociRaptor lOOOORPM with 16MB cache and 4.2ms seek time and a NVIDIA ^® 8800 GTX GPU.

Table 2: Latencies measured on the 185MB-large distance table .

Table 3: Latencies measured on the 2GB-large distance table.

3D Hilbert indirections Linear indirections

Slice XY 2.46ms 1.32ms

Slice YZ 3. I ms 4.95ms

Slice ZX 1.3ms 2.97ms

DVR 3.04ms 4. 15ms