3D base: A geometrical data base system for the analysis and visualisation of 3D-shapes obtained from parallel serial sections including three different geomerical representations

(1)

0895-61 I l/93 S6.00 + .I0

Copyright 0 1993 Pcrgamon Press Ltd.

3D BASE: A GEOMETRICAL DATA BASE SYSTEM FOR THE ANALYSIS AND VISUALISATION OF 3D-SHAPES OBTAINED FROM PARALLEL SERIAL SECTIONS INCLUDING THREE

DIFFERENT GEOMETRICAL REPRESENTATIONS

F J

. . Verbeek*+*,

M. M.

de Groat*,

D. P.

Huijsmans”, W. H. Lamers*, and I. T. Young+

* Department of Anatomy and Embryology, University of Amsterdam, The Netherlands + Pattern Recognition Group, Delft University of Technology, The Netherlands

5 computer Science Department, Leiden University, The Netherlands

(Received 21 October I992)

Abstract-In this paper we discuss a geometrical data base that includes three different geometrical representations of one and the same reconstructed 3D shape: the contour-pile, the voxel enumeration, and the triangulation of a surface. The data base is tailored for 3D shapes obtained from plan-parallel serial sections. It is explained how this geometrical data base is useful with the different processing approaches of a 3D shape, such as analysis and visualisation. Methods of conversion between the geometrical representations are discussed. Examples of the op- eration of the data base as it is embedded in a data base management system are given by illustrations of retrieval of geometrical information.

Key Words: Three-dimensional, 3D reconstruction, Serial sectioning, Geometrical data base, Visualisation, Triangulation-model, Voxel model, Contour model

INTRODUCIION

In biomedical research, three-dimensional (3D) data are collected in order to be able to analyse morphological parameters of a 3D phenomenon. A 3D phenomenon is understood to be anything that can be visual&d in the third dimension by means of an imaging technique. If clearly defined boundaries can be attributed to a 3D phenomenon, we define it as a 3D shape.

Analysis frequently aims at acquiring insight into the complexity of the 3D shape under study. This can be dealt with by means of a numerical analysis and/or a visual analysis of 3D data rendered on a computer screen.

Our research efforts focus on 3D data representing a 3D shape obtained from plan-parallel serial sections, especially those in the range of microscopical dimen- sions. These sections can either be virtual sections [i.e., made with scanning devices, e.g., the confocal laser scanning microscope (CLSM)] or physical sections (i.e., made with microtomes).

Collecting 3D data can thus be dealt with in the image-domain ( 1, 2) or by means of 3D-reconstruction

*Correspondence should be addressed to Fons J. Verbeek. Laboratory of Anatomy and Embryology, University of Amsterdam, Meibergdreef 15. 1105 AZ Amsterdam, The Netherlands.

techniques from polygonal lines (3,4) derived from an image of the 3D shape under study. If the 3D data consists of the plain data together with additional knowledge of the data (e.g., resolution in X, y, and z dimension or the structures the 3D-shape is composed of), we consider this as a 3D dataset.

In the input phase, the collected 3D dataset is represented by one particular geometrical representation. This might not be the ideal geometrical representation for a special task at hand (e.g., realistic display, image analysis). We state that the 3D dataset as collected during the input stage should not restrict the possibilities at the processing stage. Therefore, we have developed a data base that incorporates three different geometrical representations of one and the same 3D shape (5-8): contour pile, volume enumeration (voxel), and surface patch (triangulation). These three geometrical representations cover the geometrical representation domain needed for the different analytical approaches to the 3D shape, because each of the representations has its special advantages. The qualities of the geometrical representation models with respect to different processing approaches are listed in Table 1. A qualitative judgement of each of these approaches is given with respect to the geometrical representations. To come to such a judgement we have assumed an implementation on a standard graphical workstation. Additional hard-

(2)

152 Computerized Medical Imagmg and Graphics May-June/1993, Volume 17, Number 3

Table 1. Qualitative analysis of the different processing approaches of a 3D-shape with respect to the three geometrical representation models

Data input Boolean algebra Numerical analysis Fast rendering Realistic rendering Electronic knife

Contour model ++ f (high level) + ++

Volume model ++ ++ (low level) ++ _ ++ ++

Surface-patch model - - f + + ++ +

The explanations of the listed qualities are given in the introduction. The triangulation model is chosen as a surface-patch model. + + =kxcellent; + = good; + = moderate; - = poor, - - = bad.

ware (e.g., for visualisation) is not taken into account as we aim at a platform that is machine independent. For each of the data manipulations, we considered the goodness of behaviour with respect to a geometrical representation, taking into account the complexity of implementation, central processing unit (CPU)-load, and accuracy (with respect to numerical analysis). In addition, all considerations have been checked for the extent to which the specific functionality of that geometrical representation was appropriate for the processing approach. For the input of data, only two of the models are appropriate. Boolean operations can be applied to all representation models discussed. The selection of a defined entity in the 3D dataset, a structure or a slice (this is explained in more detail in Part I) is defined as a high level Boolean operation, whereas the intersection of an arbitrary Boolean shape with the shape as existing in the geometrical representation model is referred to as a low level Boolean operation. The electronic knife can be considered as a special case of low-level Boolean operations as the shape is inter- sected by one or more half-spaces. The voxel representation is the most favourable representation for these types of operations, as the whole space is defined on a cubic grid. For the same reasons, numerical analysis of a shape is favourably done on the voxel model. However, some morphological parameters (surface-related measures) can be calculated from the surface- patch model. Numerical analysis on the contour model is closely related to the concept of the spatial grid. Therefore, there is no unique use of the contour model with respect to numerical analysis in the third-dimension. Visualisation is divided in fast and realistic rendering. In doing this, we were able to qualify the merits for each of the geometrical representations best (see also Part III). The conclusions from Table I were in- cluded in an information analysis of geometrical representation models for computer-assisted reconstruction described in Part 1.

It should be noted that the data base described in this paper is meant for visualisation and analysis, and that it requires the data to be segmented and aligned. This means firstly that the components present in the

geometrical representations are separated from their background and must be uniquely labelled, and sec- ondly that 3D reconstruction of the serial sections with respect to the correct realignment and possibly other introduced deformations is already accomplished.

This geometric data base efficiently stores and re- trieves the administrative and spatial data of the different descriptions of the same object. In order to achieve this we have embedded the data base in a data base management system, called 3D Base.

This paper is organised into four parts. Part I ex- plains the design and implementation of the data base. As we have to account for the presence of the three different geometrical representations, the three types of conversion methods that had to be added to the data base are discussed in Part II. To be able to check re- trieved and/or converted geometrical data, visual inspection tools are added. The possibilities of visual inspection and the organisation of the management of the data base are described in Part III. In Part III we also pay attention to aspects of the programming environment. Finally, in Part IV, results are discussed together with the merits of the data base and the data base management. Furthermore, the direction of future developments is given.

PART I: DESIGN OF THE DATA BASE AND

IMPLEMENTATION REQUIREMENTS

This section addresses the design of a 3D data base containing three different geometrical representations as derived from an information analysis. From the results of the information analysis, the requirements for the implementation are established in terms of file for- mats and encoding of the data.

Data base design

Through a review of the literature (9) and inter- views with potential users we have provided ourselves with a basis for an information analysis of the needs of the 3D imaging community with respect to computer-assisted reconstruction, analysis and display of plan-parallel geometric data. The Nijssens Information

(3)

3D base ??F. J. VERBEEK et al. 153 Analysis Method (NIAM) ( 10) and normalisation ap-

proach (11) were used to come to an informal design of a geometric data base. After an informal logic design and feedback with the potential users we developed a coherent physical design and implementation of the geometrical data base. The physical design is the last phase of the information analysis after generalisation, differentiation, and aggregation and consists of making the decisions on the file- and record-structure to be derived from the formal descriptions. In this paragraph, we will firstly show the more informal aspects of the design of the 3D data base and gradually become more specific about its physical implementation.

The quintessence of the data base is the division of the data into geometrical representation-variant and -invariant data, together with a hierarchical description of the 3D shape. We considered the 3D shape under study as being one object that may be composed of one or more structures. As we describe the 3D shape as such, we can define a key entity for the description of a 3D shape in the data base. This key entity, the 3D object, is informally described as:

having a name;

being constructed from a series of slices [two-dimensional (2D) plan-parallel planes], defined as 2D

slice;

being possibly subdivided into structures, defined as 3D structure;

having numerical object data, such as the xy-sampling grid and the magnifications used at recording; and

having representation-specific (variant) data.

In the definition of the 3D object, two constitutive entities are introduced: the 2D slice and the 3D structure. These constitutive entities followed from the break- down of the 3D shape into elementary units. The 2D slice is informally described as:

??being uniquely represented by a plane-number;

??having a thickness;

??being present or absent;

??being selected or not selected; and

??having representation-specific (variant) data.

And the 31) structure is informally described as:

having a name;

being uniquely represented by a [red, green, blue (RGB) colour] label;

being constructed from structure-intersection data; occurring in a range of 2D slices;

being selected or not selected; and

having representation-specific (variant) data.

The key entity together with the two constitutive entities shall from now on be referred to as the basic entities. The invariant parts of the basic entities contain the generic administrative data.

The 3D structure is related to the 2D slices by the field that states the range of slices in which it occurs. If a slice is not available for reconstruction, it is said to be absent. As can be read from the informal descriptions, both constitutive entities, the 2D slice and the 3D structure, have a toggle status, which is used to make a selection from the 3D dataset. As we envisage the data base as sets in the mathematical sense, the union of the constitutive entities, either plane-wise or structure-wise, represents the 3D shape.

The variant part of the data base consists of administrative data of a specific representation as well as the data itself. The basic entities contain a field through which the representation-variant information is addressed. Depending on the representation chosen, this is an offset to either contours (contour pile), voxels (voxel), or triangle-strips (triangulation model). The variant information itself does not directly address the lowest addressable entities in the data base, but as we move down in the hierarchy the amount of available geometrical information grows rapidly until the basic elements are met. So starting at the level of the 3D object and proceeding via 2D slice information and 3D structure information to the voxels, the contours or the triangles, we accumulate more and more geometrical information about the 3D shape.

As we have chosen to work with three different geometrical representations, each of the three basic entities itself holds three entities that embody the representation-variant information. As such, we distinguish: 3D Object-Contour Model; 3D Object-Voxel Model; 3D Object-Triangulation Model, and: 3D Structure- Contour Model; 3D Structure-Voxel Model; 3D Struc- ture-Triangulation Model, and: 2D Slice-Contour Model; 2D Slice-Voxel Model: 2D Slice-Triangulation Model.

For the general understanding of the design of this geometrical data base we do not consider it necessary to list the complete description of the aforementioned entities. Therefore, we will select just those entities that are useful for the explanation of the principles of design and operation of this data base. To that end, we will elaborate on the description of the contour-representation. Each of the representations has a representation- specific entity, this is the 2D Contour for the contour model, the Voxel for the voxel model and the triangle for the triangulation model. Informally the 2D contour

is described as:

(4)

154 Computerized Medical Imaging and Graphics May-June/ 1993, Volume 17, Number 3

having a structure-membership reference; having a section-membership reference; having a start point;

being a list of connected contour-points; being either open or closed; and

being an inner-contour of a 3D structure, an outer- contour of a 3D structure or an inner-contour of one 3D structure and at the same time an outer contour of another 3D structure.

The last two fields of this description refer to the topology of the contour. In most cases a structure of interest (anatomical structure) consists of closed contours. The boundaries of the structures that lie on the border of the field of view often have significance as anatomical reference of the object under study. To enable the use of these boundaries the contour-lines of the partially visible structures have to be indicated as open contours. A structure-intersection may have, apart from its outer boundary, one or more internal boundaries. These internal boundaries result from cavities in and/or concavities of the structure under study. In a visualisation of such a structure they will appear as such. These external and internal boundaries are referred to as nested topologies. Nested topologies are stored in the data base as inner and outer contours. The bridge between the basic entities and the 2D contour are the representation-variant entities. These are described as follows:

3D object-contour model

??has a contour list;

??has numeric information in spatial distances in x, Y, z;

??has a coding method; and,

3D structure-contour model

??has a bounding-box; and,

2D slice-contour model

??has a first contour number and a last contour number; and

??has a bounding box.

From the given informal descriptions one can see how geometrical information is addressed while de- scending in the hierarchy of the structure under study. Imagine, for example, that we want to select all the contours of one particular 3D structure within the 3D

object, the following steps are to be made:

??on 3D object level the contour list is addressed;

??on 3D structure level the structure is selected;

??on 3D structure level the range of sections is selected,

??on 2D slice level within the range of sections the representation-variant information is addressed (per slice, all contours having the structure-membership of the selected structure are selected); and

??per selected 2D Contour a connected string of contour points is addressed.

The last example is strictly abstract. In the section on applications, the principles of operation of the data base are illustrated.

Data base implementation

The final implementation of the data base consists of the transformation of the informal descriptions via a formal design into a record structure and a file format. The basic entities and the representation-specific entities are stored in separate files. The files store the rec- ords derived from the aforementioned entities sequen- tially. The files containing the representation-variant information are referred to by their file names, which are to be found in the file containing only one record, being the 3D Object (the key-entity).

The lowest addressable geometrical information is stored differently for each of the representations. The representation-variant administration contains information on the coding technique that is used to meet the requirement of efficient storage of the data base. Coding techniques of representation-variant data were chosen in such a way that the use of memory and CPU- load were minimal at decoding.

The connected contour-strings should satisfy the requirement of being either 4-connected or 8-con- netted. Thus, the connected strings can be seen as either 4- or 8-connected Freeman code ( 12). Contour-strings that are 4-connected (Freeman crack code), are coded by 2-bits per coordinate, whereas 8-connected contour- strings are coded by 4-bits per coordinate. The 4-bits represent the permutations of the differential chain- code ( 13). Whenever the requirement of connectivity is not met, interpolation on a square grid is applied. To accomplish the interpolation we have implemented the Bresenham line algorithm (14), taking the type of connectivity as an argument and thus resulting in a 4- connected or an 8-connected string of coordinates.

Binary voxel planes of a single 3D structure can be coded with quadtrees ( 15) and these quadtrees can then be merged to octtrees (16). However, the voxels can also be stored as a labelled voxel image (8-bit values). Triangles are stored as vertices. As a vertex will be a part of at least three triangles, the storage is chosen such that redundant storage is minim&d. The type of coding chosen for a particular representation is stored on the representation-variant part of that specific representation, so that, upon retrieval, decoding is exe- cuted correspondingly.

(5)

3D base ??F. J. VERBEEK et al. 155

PART II: CONVERSION OF GEOMETRICAL INFORMATION

In this section the conversion methods of the geometrical representations are discussed. Apart from a brief description of the algorithms of a conversion, motivations are formulated.

Conversion strategies and methods

Of the three geometric representations mentioned, the contour model and the voxel model are the only representations that can be used at the phase at which the data are collected (l-4). For the data base, these representations serve as input channels. If and when either one of these two representations is available, the other two geometric representations are constructed from this initial shape representation by conversion routines. As we wanted to implement the minimum number of conversion routines necessary to derive all three representations, three conversion routines are needed:

1. from contour model to voxel model; 2. from voxel model to contour model; and 3. from contour model to triangulation model.

The triangulation cannot occur as an input model, and therefore no conversions from the triangulation representation to the other representations are called for. The complex topologies of biological shapes requires that the conversion routines are robust. No user interaction should be involved with the conversion routines, and therefore these routines should operate fully automatically. These constraints imply the use of heuristic methods in the representation conversions. Furthermore, we demand an efficient use of CPU-time and minimal use of memory. In the next part the methods used in the conversion routines are discussed.

Conversion from contour model to voxel model.

The conversion from contour model to voxel model could be dealt with in three ways. We distinguish:

a polygon-fill routine distributed with a graphics software library;

cellular logic image-processing routines; and a scanline-based algorithm, which is an improved version of an earlier algorithm developed by Huijs- mans (17).

The first method mentioned depends on the windowing environment supplied with the graphics library. Moreover, as the contours from the contour model are complex polygons and the structures may have complex topologies, it appeared that the implementation would probably be far from straightforward.

The second method is based upon image-processing routines and will, undoubtedly, work well. How- ever, as we have to address complex topologies, this method may lead to more memory requirements and processing time than would be strictly necessary.

The third method is the best suited for this conversion, because it was especially designed for the scan- conversion of closed contours of complex topologies to planes (17). The flavour of the algorithm can be summarised as follows:

define the bounding box of the slice, that is, the bounding box that encloses all contours in the section;

scan each line for the presence of contour-pixels, and accumulate those pixels. A contour point is considered a vector of three elements. An “x” and a “y” coordinate and a structure-membership code; group the pixels according to the structure-membership. Separate stacks for each of the structures are made. Everytime a new structure is met, a new stack is initiated;

code and sort the pixels;

remove redundant codes (only the end-points of connected pixels of a contour on a scanline are stored); and

fill a scanline in a voxel plane.

It should be stressed that the starting point of this algorithm is a connected contour (i.e., a string of pixels). The connectivity is either 4-connected or 8-connected. We prefer to work with 8-connected contours as these contour-strings consist of fewer points. It should finally be noted that this algorithm has to be regarded as a version of the “painters-algorithm” ( 14) for nested topologies and that it can therefore be used for hidden- line removal as well.

Conversion from voxel model to contour model.

Conversion from a voxel model to a contour model is realised through image-processing techniques. In con- trast to the conversion from contours to voxels this conversion method is applied in the image-domain, as voxels are defined in the image-domain. As the order of the consecutive voxel planes corresponds to the direction of scanning, the planes are processed in that order (i.e., plane by plane). Thus, a voxel plane is copied into a 2D image and thus the conversion is reduced to a 2D image processing problem.

The conversion of one voxel plane to contours in the contour model is described by the following algorithm:

??repeatedly select a labelled component (the labels are corresponding to the structures defined in the data base), then

(6)

156 Computerized Medical Imaging and Graphics May-June/1993, Volume 17, Number 3

label all the selected structure-intersections sepa- rately, then

for each label apply a threshold (at the label value), then

extract the contour in the resulting binary-image (the contour is here defined as the set of pixels that is connected to the background and contour extraction is realised through cellular logic image operations), then

apply a contour follower (here a contour in an image is converted to a connected-string of xy-coordinates), and

determine bounding-box of all contours (of the selected structure-intersection), and subsequently, determine the outer/inner contours (by comparison to the largest) as well as the direction of the contour from the bounding-boxes, and

store each contour accordingly.

The complexity of this conversion depends on the degree of nested topologies in the original 3D shape. However, from the given algorithm it can be concluded that nested topologies are easily found because the topology-state of the contour (being an inner or an outer contour) is determined by comparison to the largest bounding box of the structure-intersections present in the sub-image. The memory requirements of this algorithm are moderate, as it needs allocation of one image-buffer of the size of one voxel plane and a sub- image of variable magnitude. The maximum size of the variable sub-image equals the largest bounding box of the structure-intersection found in the voxel plane that is being processed.

Determination of the direction of the contour is necessary for the conversion algorithm used to convert to the triangle-representation. Thus, it is assured to have the right direction at storage. Contours can be stored either g-connected or 4-connected. At this stage we preferred to store the contours as g-connected strings.

Conversion from contour model to triangulation model. As surface-patch we have chosen the triangle

for its simple analytical qualities. The conversion to the triangulation representation is done from the contour description (contour model) of the object (18- 22) and therefore the only representation conversion to the triangulation representation emanates from the contour pile representation.

The conversion routine from contour representation to surface-patch representation (triangles) is the most difficult one because the reconstruction of the surface is based upon spatial information along the z axis, whereas the other conversions are only based upon spatial information in the xy plane. In the case of serial- sectioning, the sampling density in the xy plane does,

in most cases, not correspond with the sampling density along the direction of sectioning, the z axis. This means that connectivity between structure-intersections in successive planes is not per definition a fact. The ques- tion of connectivity arises in particular at locations of branching and/or concave structures. In order to be able to reconstruct the surface of the structure automatically, more information than present in consecutive planes of a contour model or a voxel model is needed. Conversion from the voxel model to the surface patch model is therefore real&d via the contour model, because in the contour model there are only contours in 3D space and connectivity can be calculated by rules. The algorithm for triangulation of a surface from a contour description developed by Boissonnat (23) provides such rules. The strong mathematical base’ of this method makes it possible to deal with complex topologies such as branching and concavities. The algorithm of Boissonnat is based on the maximisation of a volume between contours of the same structure in consecutive planes (23). A maximum volume is reconstructed from the contour-data by a 3D Delaunay triangulation, which can consequently be regarded as a 3D tetrahe- dralisation of a bounded space. From the volume reconstruction the object surface (triangles) is derived heuristically. Except for severely undersampled data, this method normally leads to an acceptable surface reconstruction. Results of this triangulation method are shown in Fig. 1C and discussed in the section on application (Part IV).

PART III: ANALYSIS OF THE DATA

In the third part of this paper we elaborate on the possibilities of the three different geometrical representations with respect to 3D image analysis and visualisation.

Considerations on visualisation

One of the aims of the geometrical data base is visualisation of the 3D shape under study. There are many ways in which a visualisation of a 3D shape can be realised. One of the objectives of the integration of three different geometrical representations into one data base is to be able to upgrade a simple visualisation of the 3D shape with a more advanced one. A qualitative evaluation of the possibilities of the visualisation (and analysis) with respect to the different geometrical representations is given in Table 1. Because all three geometrical representations are available for the display-generation, a correct viewpoint for the rendering

’ The method of Boissonnat is fully embedded in the concepts of Graph-Theory and Computational Geometry.

(7)

3D base ??F. J. VERBEEK ef al. 157

input layer:

check databane conrirtency I

conversion layer:

make all geometric representations available

inspection

Y

inspect database inspect geometry

shape analyeir

Fig. 1. Illustration of the processing environment of 3D base. The contour model and the voxel model can be used as input-representation model. These input data are translated as a data base to 3D base. It is advisable to check the consistency of the data at this stage. After having imported either the contour model or the voxel model the other geometric representations are created. After conversion, three main processes are distinguished: rendering, inspection and analysis (numerical). Rendering and inspection are sub-divided in processes that are typical to a special geometrical representation (surface rendering, volume rendering) or perform best with a particular geometrical representation (data base inspection). The processes listed (depicted by blocks) and the representation models (depicted by circles) are connected by lines. From (connected) process to model one should read: “is favourably

done by”; whereas from (connected) model to a process one should read: “can be used efficiently for.”

of the 3D shape under study can be generated using the contour model, as the contour model is the proper representation for quick rendering. Once a good viewpoint is found, the voxel model and/or the triangulation model can be used to creale a more realistic view of the 3D shape. Rendering of a voxel model on a computer display will clearly show the effect of noncubic voxels that are a result of the undersampling in the z axis during serial sectioning. Because the noncubic voxels directly reflect the thickness of the serial section from which the voxel planes originate, such a display is regarded as realistic. Rendering of a triangulation model on a computer display may result in a smooth surface. The smoothness depends on the smoothness of the contours at conversion to the triangulation model

as well as the type of smoothing (and illumination) that is applied subsequently to the surface (14,24,25). A smooth surface generally corresponds best with the concept of the user regarding a surface of an anatomical shape. This conceptual illustrative power of the display of a triangulation model is the reason why advanced rendering of the triangulation model is considered realistic.

Visualisation can also be used as an inspection tool of the contents of the data base. Because a 3D shape may consist of one or more structures, an inspection by means of display is used to make a selection of the structures necessary to illustrate the 3D shape under study. The same holds true for the selection of the sections. Again the contour model is the proper

(8)

158 Computerized Medical Imaging and Graphics May-June/1993, Volume 17, Number 3

representation for the assessment of the contents of a geometric model of the 3D shape.

As the contour model appears to be very useful for a preliminary visualisation of the 3D shape under study, we have explored its display possibilities. As such, we distinguish three display modes of the contour model:

1. direct rendering of the contours of the contour model;

2. hidden-lines calculation (see conversion from contour model to voxel model), and then rendering of the contours; and

3. hidden-lines calculation, as in mode 2, but now the surface, as enclosed by the structure-intersection and filled with the structure colour label, is rendered. This results in a pseudosolid appearance.

Without increasing processing time, the rendering quality of the contour model is improved by using the aforementioned display modes. In Figs. 2A, 2B, and 2D examples of visual inspection of the geometrical contents of the data base are given.

Data base management by the integration of the data base, geometric conversion and inspection visualisation

The implementation of the 3D data base together with the conversion routines and visualisation consti- tute a data base management system for manipulation of 3D geometric models from plan parallel slices. This management system consists of different layers. These layers are depicted in Fig. 1. The data base management system is meant to serve as a buffer that exports and/ or imports the data base.

The contour model and the voxel model serve as input models for the data base. Once the input model is available, the conversions to the other models can be effectuated. The selection fields of the 3D data base are used to make a selection of slices and structures to work with. The data base organisation is such that selection of a structure automatically means that the range of slices in which structure-intersections of this structure are found, is selected. One can further restrict the display of the structure by toggling a range of sections to the “nonselected” state. As these selections are representation-invariant, all three different geometric representations will use the same selection. In Table 1, selections of structures and slices are referred to as high- level Boolean algebra, because slices and/or structures are selected on the logic of their toggle status.

We concluded that the contour model is best qualified as a tool for the fast visual inspection. We have therefore added a visual inspection layer on top of the conversion layer using the different display

modes of the contour model. It should be noted that the three display modes of the contour model may be used simultaneously. A display mode is then attached to a 3D structure.

Besides visual inspection of the contents of the data base, an inspection of the consistency of the geometrical representation models, as obtained from conversion of geometrical information (see also Part II), is considered very useful. This is realised by a visual inspection of the results of the conversion routines. Such a visual inspection requires fast and informative displays of all geometrical representations. The display possibilities of the contour model satisfy this requirement. In order to attain a visual inspection possibility of the voxel model and the triangulation model, we have added simple display-routines for the voxel model and the triangulation model. Fast display of the voxel model is attained by a visualisation of the sides of the bounding-cube it is defined in. The voxel model can thus be inspected from six orientations. The triangulation model is visualised as a wire-frame of triangles and can be given any orientation. Because the contour model and the triangles of the triangulation model may be rendered on the 2D computer display in any orientation, all views of the 3D shape are generated by means of perspective transformation.

If preliminary inspection has proven the representation model to be correct and selections of orientation and regions of interest are made, an advanced display program may be used for the generation of a realistic display (Fig. 1) either by volume rendering (voxel model) or by surface rendering (triangulation model).

It should be noted that, apart from retrieval of geometrical information by visual inspection, the al- pha-numerical information in the data base can be re- trieved by the generation of lists.

Advanced visualisation and 30 image analysis It has not been our intention to build a software package including advanced graphics. Therefore, we interfaced the data base with programs that are specif- ically made for this purpose or that used 3D graphics libraries that support such advanced visualisation techniques.

With respect to 3D image analysis we have taken care that the data base can be interfaced to software packages which are better “equipped” for that special task. This is in accordance with the niche of 3D base as being a buffer in an analysis environment, which is depicted in Fig. 1. In most cases, the voxel model is used for the analysis of the 3D shape (Table 1). How- ever, as the triangulation model is in fact an estimation of the surface of a 3D shape (or a part thereof), this

(9)

3D base. F. J. VERBEEK et al. 159 Fig 3D All fix rep (a) (b) 64

2. Examples of retrieval of geometrical information from the data base. Different visualisation of parts of a “reconstruction of a rat of 1 I Embryonic Days. (a) The neural tube, the gut and the heart (pseudo-solid); (b) structures and sections displayed as pseudo-solid contours; (c) Triangulated surface of the neural tube; (d) Wire- me visualisation of the triangles from panel 2C in a different orientation; (e) Volume rendering of a voxel- lresentation of the same rat embryo; and (f) Volume rendering of the heart of the same embryo with lateral cut-

(10)

160 Computerized Medical Imaging and Graphics May-June/l993, Volume 17, Number 3

model might be useful in the estimation of surface- area. Especially in cases of undersampling along the z axis, other methods may be less accurate (26). The measures used for the description of a 3D shape are considered to be out of the scope of this paper. There- fore, the software environment that is used to generate/ calculate such measures will be discussed in the next section,

Programming environment and software platform of 30 base

The data base, the conversion routines, and the visualisation tools are programmed in the “C” language on a Sun Spare II running SunOs 4.1.1. (Sun Corp., Ichinomiya, Japan) Because 3D base is meant to serve a larger community of users, we have chosen X 11 (27) as a basis for the user interface, in order to guarantee portability amongst the majority of graphical worksta- tions (e.g., Sun, Apollo, HP, Silicon Graphics, IBM and DEC). The libraries supplied with the public domain version of Xl 1 also provide support for the 2D graphics needed for the visual inspection of the data base. The Programmars Hierarchical Interactive Graphics System (PHIGS) (28) subroutine libraries with support for 3D graphics as needed for realistic rendering were used to generate a smoothed surface display of the triangulation model. As from release 5 of the X 11 environment, PHIGS is an integral part of X 11 and often referred to as PEX (PHIGS Extensions to X) (28). Results that were generated by programs based on PHIGS are based on the SunView windowing environment. The volume rendering was accomplished with the libraries supplied with the SunVision software.

For the numerical analysis of 3D shapes we use the SCIL-Image environment (29). SCIGImage is an X 11 -based image processing package that incorporates a C-source interpreter. Furthermore, it provides the possibility to extend the program with “home-brewed” subroutine libraries. These features were used to their full extent as the integration of data base in the package (SCIL-Image) enabled us to use the C-interpreter for the prototyping of shape-analysis measures.

PART IV: THE APPLICATIONS AND DISCUSSION

In the last part of this paper we show results that are realised with the developments described in Part I (the data base), Part II (the conversion methods), and Part III (the visualisation possibilities).

Applications

At the2 Department of Anatomy and Embryology the changes in the spatial distribution of gene expres-

2 Applications are discussed in terms of research at the De- partment of Anatomy and Embryology, University of Amsterdam, The Netherlands.

sion in embryonic organs during development are studied. The gene expression patterns serve as molec- ular markers to elucidate the tissue remodelling that underlies normal and abnormal development (30,3 1). Gene-expression patterns are visualised in histological sections by immunological methods. Both input channels ( l-4) of the data base, the contour pile and labelled voxel model, are used with computer-assisted 3D reconstruction from the sections. The 3D reconstruction with contours is based on digitisation of contours by manual tracing (1, 2) and is used for quick 3D reconstruction from serial sections. A voxel model is generated by video digitisation and respective image reg- istration of the consecutive serial sections (3, 4).

In order to explain the various features of 3D Base mentioned in the preceding parts of this paper we have taken a straightforward example. A rat embryo of 11 embryonic days (ED) was dissected from the uterus and, after histological pretreatment, serially sectioned. Micrographs were made of a subset of 50 slices. Struc- tures of interest were digit&d by manually tracing the structure intersections in the successive slices. Realign- ment of the stack of contours resulted in a contour model. With this contour model as a basis, the results depicted in Fig. 2 were achieved. In terms of the data base the 3D shape under study is the 11 ED rat embryo, which in the data base is referred as the 3D object. As the embryo consists of a number of structures, such as the heart, gut, neural tube, veins, and arteries, it is of no use to display them all at once. Hence, selections of a few structures were made. Figure 2A illustrates such a selection, which also illustrates the use of a mixed display of contours. In Fig. 2B only the filled contours were used. Both displays clearly show contours oriented in 3D space.

In Figs. 2C and 2D, the results from the conversion of the contour model to the triangulation model are shown. In these two figures only the neural-tube structure is selected. Figure 2C shows the results of a display generated by a program for surface rendering based on PHIGS. No smoothing technique is applied in order to illustrate the surface as constructed from triangle- tiles. The conversion from contour model to triangulation model has accomplished a correct tiling of the surface as can be concluded from Figs. 2C and 2D. The branching of the optic-vesicle from the neural tube is handled correctly, as is the spiral shape of the neural tube. Figure 2D shows the wire frame of triangles that is used for the inspection of the conversion to the triangulation model.

The results regarding the visualisation of the voxel model are shown in Figs. 2E and 2F. Both displays were generated using volume rendering. The display of the voxel model, in correspondence with its nature,

(11)

3D base ??F. J. VERBEEK et al. 161

is solid. Cavities in structures are therefore not always visible. In order to visualise concealed structures, or parts thereof, the cut-away view is used. An electronic knife is used to leave out certain parts of the 3D data. The principle of the electronic knife is considered low- level Boolean algebra. As the voxel model is a continuous description of the 3D shape, this principle is very well applicable to this geometric representation (Table 1). In Fig. 2F an application of the electronic knife is shown. The embryonic heart is the only structure selected, and a subset is selected from the range of sections in which the heart is found. We look on the top of the range of selected slices and observe a lateral cut-away view of the ventricle of the heart. The principle of the electronic knife is also applicable to the triangulation model. However, as this representation is not contin-- uous, discontinuity of the surface has to be ac- counted for.

DISCUSSION

We have presented the development of a geometrical data base that integrates three different geometrical representation models. This data base enables us to use the advantageous aspects of each of the described geometrical representation models, which is achieved by a data base management system that includes conversion routines (described in Part II) and visual inspection tools (described in Part III). As we also realised inter- faces of the data base with special purpose software (advanced display, 3D Image analysis) a firm basis for different processing approaches necessary for our research-line was established.

Starting from either the contour model or the voxel model a maximum of two (out of three) conversion routines are necessary in order to have all three geometrical representations stored in the data base. The conversion of the voxel model to the triangulation model via the contour model can be considered a detour. This seems to be true if the visual inspection tools, based on the contour model, are not used. Conversions from the voxel model to a surface model use an extra processing step before the surfaces are derived (32,33). Most of these methods aim at rapid display generation of voxel images by reducing the number of voxels to be processed to just the surface voxels (surface tracking). Only the “marching cubes” algorithm (34) con- verts surface voxels to a tessellation of triangles. This algorithm was developed for rapid display of computer tomography (CT) scans and starts with the calculation of the surface voxels from its grey-valued voxel data. A tessellation of the surface voxels is realised but, for display-generation, the gradients derived from the grey- values are used in combination with the triangle-ver-

tices for the shading effect of the surface. In our case, edges are already defined as we start from a voxel model containing labelled data. Therefore, conversion of the (labelled) voxel model via the contour model is in fact an efficient detour, because in our case it omits the unnecessary climbing of a mountain.

The triangulation method of Boissonnat (23) was chosen for its ability to deal with complex situations. Other triangulation methods failed in the automatic reconstruction of the surface from contour lines in such complex cases ( 19-22) or required user-interaction (18). It should be noted that surface reconstruction is always dependent on the sampling in the z axis. At severe undersampling all automatic methods will fail. In order to produce a satisfactory result with the triangulation method of Boissonnat, resampling of the contours is needed, otherwise too many triangles will be generated. To that end, a resampling scheme that takes into account the local curvature must be used. One should be aware of the fact that resampling may result in a coarse polygonal approximation of a smooth contour, which in turn results in a coarse triangulation (surface patch representation). Starting form a pixelwise description of the contour, (string of connected coordinates) the method of resampling of the contour is a crucial step for the triangulation method.

As we have mentioned in the introduction, we start from voxel data that are already classified by seg- mentation and labelling. We prefer to refer to these voxel models as labelled binary voxel models. Such a model is thought of as constituted by the union of the binary voxel models of the subsequent labels. Appli- cation of volume rendering techniques to these labelled binary voxel models is referred to as binary voxel processing (35). The labels found in the voxel model should therefore be attached to a RGB colour for visualisation. As the representation invariant information of the 3D- structures contains a RGB colour-code of the structure and as the structures present in the data base correspond with the labels of the voxel model, the attach- ment of a RGB colour to a label is intrinsically established.

Fast display generation can be achieved with voxel models using special hardware (36). The inspection tools offer a software solution for the fast display generation. We have chosen this solution for reasons of portability of the software to other sites. As far as rapid display is concerned, we think much can be learned from the approaches taken with the rapid display of voxel data. The use of precalculated look-up tables will be studied for future addition of more rapid display facilities to the data base management system. At the same time, the upgraded functionality of Xl 1 with PEX (28) can be used for the fast generation of

(12)

162 Computerized Medical Imaging and Graphics May-June/l993, Volume 17, Number 3

smoothed triangulation models. Display possibilities using an opaque visualisation of surfaces will be in- corporated as well. If vendors supply X 11 -based drivers with their dedicated hardware, then in the future dedicated hardware might be used for interactive advanced display generation, while portability of the developed software can be maintained optimally.

for critical reading of the manuscript. This research is partially sup- ported by the Netherlands’ Project Team for Computer Science Re- search, SPIN (project Three-Dimensional Image Analysis).

REFERENCES

1.

Undersampling of the data along the z axis is al- most always present in serial sectioning, and current research is directed toward extensions of the method of Boissonnat with rules that enable the detection of a failure of the method as a result of undersampling. In combination with these rules we think that interpolation of consecutive planes (containing the contours or the voxels) to one or more intermediate planes will contribute to the solution of this pitfall. We think that our expertise with dynamic interpolation techniques (37, 38) regarding this problem will be valuable.

2.

3.

4. 5.

In conclusion, we think that the geometrical data base can be an essential instrument for the exchange of 3D shape information and algorithms for shape analysis. As we have chosen not to give a complete description of the geometrical data base in Part I, we offer, on request, the C-listings of the format of the data base. 6. 7. 8. 9. SUMMARY

The requirements of a system for the management of 3D data obtained from serial sectioning, have been discussed. A data base that integrates three different geometrical representations has been developed. The requirements for the management system, called 3D Base, are centred around the contents of the data base. The geometrical data base consists of representation- variant and -invariant data. Conversion routines are implemented to accomplish all three representations in the data base: the contour pile, the voxel enumeration, and the triangulated surface. The methods of conversion from one representation to another are discussed. For the visual inspection of the data base, sev- eral visualisation possibilities were evaluated and in- tegrated in the management system. Finally, results of the operation of the management system are shown in the scope of the developed data base, conversion, and inspection tools. Some hints are given about future developments of the system.

10.

Laan, A.C.; Lamers, W.H.; Huijsmans, D.P.; te Kortschot, A.; Smith, J.; Strackee, J.; Los, J.A. Deformation corrected computer- aided threedimensional reconstruction of immunohistochemi- tally stained organs: Application to the heart during early or- ganogenesis. Anat. Rec. 224443-457; 1989.

Verbeek, F.J.; Lamers, W.H.; Young, LT. Alignment and deformation correction of serial sections of microscopical objects. Cytometry (Suppl. 4):4; 1990.

Verbeek. F.J.: Schoutsen. N.J.C.: Huiismans. D.P. A PC-based package’for three-dimensional recon&uction of microscopical objects. Cytometry (Suppl. 4): 105; 1990.

Verbeek, F.J.; Schoutsen, N.J.C.; Huijsmans, D.P. A program for 3D-reconstruction from serial sections. Microsc. Res. Tech. (in press).

Verbeek, F.J.; de Groot, M.M.; Huijsmans, D.P. 3D-Base: Data base description for 3D shape exchange. Cytometry (Suppl. 5): 123; 1991.

Huijsmans, D.P.; Jense, G.J. Representations of 3D-objects reconstructed from series ofparallel 2D-slices. In: Earnshaw, R.A., ed. Theoretical foundation of computer graphics and CAD. Ber- lin: Springer; 1988: 1031-1038.

Badler, N.; Bajcsy, R. Three dimensional representations for computer graphics and computer vision. Comput. Graph.

12:153-160; 1978.

Barillot, C.; Gibaud, B.; Luo, L.M.; Scarabin, J.A. 3D representation of anatomic structures from CT examination. SPIE 602, Biostereometrics ‘85 307-3 14; 1985.

Huijsmans, D.P.; Lamers, W.H.; Los, J.A.; Strackee, J. Toward computerized morphometric facilities: A review of 58 software packages for computer-aided 3D reconstruction, quantification and picture generation from parallel serial sections. Anat. Rec. 216:449-470; 1986.

Verbeijen, G.M.A.; van Bekkum, J. NIAM: an information analysis method. In: Olle, Sol, Verrijn-Stuart, eds. Information systems design methodologies. Amsterdam: North Holland: 1982; 537-590. II. 12. 13. 14. 15. 16. 17. 18. 19.

Acknowledgments-We have greatly benefitted from the pilot studies on geometrical representation models of D. Dekker, V. van de Hoven, and C. Schoutsen, all former students of the Computer Science De- partment (Application Oriented Computer Science) of the Leiden University, The Netherlands. We further wish to thank Mr. C. Hers- bath for photographic assistance, Mr. R. W. A. M. Baeten and Mrs. C. de Gier-de Vries for technical assistance, and Dr. A. S. Thomson

20. 21. 22. 23.

Date, C.J. An introduction to data base systems. Reading, MA; Addison-Wesley; 197 1.

Freeman, H. On the encoding of arbritary geometric configu- rations. IRE Trans. Electronic Comput. EC 10:260-268; 196 1. Bowie, J.E.; Young, IT. An analysis technique for biological shape. Acta Cytologica 21(3):455-464; 1977.

Foley, J.D.; van Dam, R.A. Fundamentals of interactive computer graphics. Reading, MA: Addison-Wesley; 1982. Samet, H. The quadtree and related hierarchical data structures. Comput. Sci. TR-I 329: University of Maryland; 1983. Meagher, D. Geometric modelling using octree encoding. Com- put. Graph. Imaging Proc. 19:129-147; 1982.

Huijsmans, D.P. Closed 2D-contour algorithms for 3D-reconstruction. In: ten Hagen, P.J.W., ed. Eurographics ‘83. Amster- dam: North Holland; 1983: 157-168.

Christiansen, H.N.; Sederberg, T.W. Conversion of complex contour line definitions into polygonal element mosaics. Comput. Graph. 12:187-192; 1978.

Ekoule, A.B.; Peyrin, F.C.; Odet, C.L. A triangulation algorithm for arbitrary shaped multiple planar contours. ACM Trans. Graph. 10(2):182-199; 1991.

Fuchs, H.; Kedem, Z.M.; Uselton, S.P. Optimal surface reconstruction from planar contours. Comm. ACM 20:693-702; 1977. Ganapathy, S.; Dennchy, T.G. A new general triangulation method for planar contours. Comput. Graph. 16:317-325; 1982. Keppel, E. Approximating complex surfaces by triangulation of contour lines. IBM J. Res. Dev. 19:2-l 1; 1975.

Boissonnat, J.D. Shape reconstruction from planar cross-sections. Comput. Vis. Graph. Imaging Proc. 44: l-29; 1988.

(13)

3D base ??F. J. VERBEEK et nl. 163 24. Gouraud, H. Continuous shading of curved surfaces. IEEE Trans.

25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 31. 38. Comput. 20(6):623-629; 1917.

Phong, B.T. Illumination for computer generated pictures. Comm. ACM 1831 I-317; 1975.

Mulhkin, J.C.; Verbeek, P.W. Surface estimation of digitized planes. BioImaging 1:6- 16; 1993.

Jones, 0. Introduction to the Window X System. New Jersey: Prentice-Hall Inc; 1989.

Howard, T.J.L.; Hewitt, W.T.; Hubbold, R.J.; Wynvas, K.M. A practical introduction to PHIGS and PHIGS+. Reading, MA: Addison-Wesley; I99 1.

tenKate,T.T.;vanBalen,R.;Smeulders,A.W.M.;Groen,F.C.A.: den Boer, G.A. SCIL-Aim: A multi-level interactive image processing environment. Pattern Recognition Lett. 11:429-441;

1990.

Lamers, W.H.; Verbeek, F.J.; Moorman, A.F.M.; Charles, R. The metabolic lobulus, a key to the architecture of the liver. In: Gumucio. J.J., ed. Hepatocyte heterogeneity and liver function RBC, Vol 19. Berlin: Springer International; 1989.

Lamers, W.H.; Wessels, A.; Verbeek, F.J.; Moorman, A.F.M.; Viragh, Sz.; Wenink, A.C.G.; Gittenberger-de Groot, A.C.; An- derson, R.H. New findings concerning ventricular septation in the human heart-their implications for maldevelopment. Cir- culation 86: 1 194- 1205; 1992.

Artzy, E.: Frieder, G.; Herman, G.T. The theory, design, implementation and evaluation of a three-dimensional surface detection algorithm. Comput. Graph. Imaging Proc. 15: l-24; 198 1. Gordon, D.; Udupa, J.K. Fast surface tracking in three-dimensional binary images. Comput. Vis. Graph. Imaging Proc. 45:196-214; 1989.

Cline, H.E.; Lorensen, W.E.; Ludke, S.; Crawford, C.R.; Teeter, B.C. Two algorithms for the three-dimensional reconstruction of tomograms. Med. Phys. 15(3):320-327; 1988.

Kriete, A.: Pepping, T. Volumetric data representations in mi- croscopy: Applications to confocal and NMR-microImaging. In: Kriete, A., ed. Visuahsation in biomedical microscopies. Wein- heim: VHC; 1992: 329-359.

Goldwasser, S.M.: Reynolds, R.A.; Talton, D.A.; Walsh, E.S. Techniques for the rapid display and manipulation of 3D biomedical data. Comput. Med. Imaging Graph. 12( l):l-24;

1988.

Verbeek, F.J.: Lamers, W.H.; Young, LT. Deformation correction by a contour-matching approach. Cytometry (Suppl. 5):22; 1991. Verbeek, F.J. Deformation correction using Euclidean contour distance maps. Proceedings 1 lth International Conference on Pattern Recognition: IEEE Comp. Sot. Press: 347-35 1; 1992.

About the Author-FoNs J. VERBEEK received his M.Sc. in Biology and Computer Science at the Agricultural University Wageningen. A substantial part of his graduation, comprehending 3D reconstruction and 3D image analysis, was performed at the Computer Science Department of the Leiden University, the Netherlands. He is currently working as a Ph.D. student in a collaboration-project between the Delft University of Technology (Pattern Recognition Group) and the

University of Amsterdam (Dpt. Anatomy and Embryology) on 3D image analysis. His main interests are 3D reconstruction, image analysis, image processing, and image warping.

About the Author-MARK0 M. DE GROOT received his M.Sc. in Medical Computer Science at the Leiden University. As a part of his masters thesis he cooperated in the research on the geometrical representation models. He is currently working as a Ph.D. student at the Department of Computer Science of the University of Utrecht. His interests: theory and application of computational geometry and computer graphics.

About the Author-DtoNYstus (NIES) P. HUUSMANS received the M.Sc. and Ph.D. degrees, both in Physics, at the University of Am- sterdam. His Ph.D. research comprehended time variations in sec- ondary cosmic rays. After completion of his Ph.D., Dr. Huijsmans spent 3 years (1982-1985) with the Laboratory of Medical Physics at the University of Amsterdam. During this period he became involved in the research in Computer-Aided 3D Reconstruction. This area of research always retained his interest. In 1985, he joined the Department Computer Science of Leiden University where he holds a position as Assistant Professor in Application Oriented Computer Science. His main interest include computer graphics, image processing and analysis and scientific visuahsation of multi-dimensional datasets.

About the Author-WOUTER H. LAMERS received the M.D. and Ph.D. degrees from the University of Amsterdam in 1974 and 1980, respectively. In 198 I-1982 he was a visiting assistant professor in the Department of Biochemistry at Case Western Reserve University, Cleveland, Ohio. After returning to Amsterdam Dr. Lamers started working on developmental and topographic aspects of the regulation of gene expression, in particular in the heart and liver. In 199 I, Dr. Lamers became Professor of Anatomy and Embryology in the De- partment of Anatomy and Embryology of the University of Am- sterdam. In 1992. he became director of research of the Centre for Liver and Gastrointestinal Research at the Academic Medical Centre of Amsterdam.

About the Author-IAN T. YOUNG received the B.S., M.Sc., and Ph.D. degrees, all in Electrical Engineering, from the Massachusetts Institute of Technology. From 1969 to 1973 he was an Assistant Professor of Electrical Engineering and from 1973 to 1979 an As- sociate Professor Electrical Engineering at MIT. From 1978 to 198 I Dr. Young was a group leader in the Biomedical Sciences Division of Lawrence Livermore National Laboratory. He has been a Visiting Professor in the Electrical Engineering Departments of Delft IJni- versity of Technology, The Netherlands in 1975-1976, Technical University of Linkoping, Sweden, 1976, and the Ecole Polytechnique Federale de Lausanne, 1979-1980. In December 198 I Dr. Young became Professor of Measurement Technology and Instrumentation Science in the Department of Applied Physics, Delft University of Technology, The Netherlands. Over the past decade Dr. Young has been a consultant to a number of companies including MIT’s Lincoln Laboratory, the Coulter Biomedical Research Corporation, and Bec- ton Dickinson. From 1977 to 198 I he was a member of the Cytology Automation Committee of the National Cancer Institute.