- Email: [email protected]
Detection and three-dimensional reconstruction of a vascular network from serial sections
Pattern Recognition Letters 1 (1983) 497-505 North-Holland
July 1983
Detection and three-dimensional reconstruction o f a vascular network from serial sections H. H.-S. I P * Image Processing Group, Department of Physics and Astronomy, University College London, London WCIE 6BT, U.K. Received 30 March 1983
Abstract: The process of three-dimensional reconstruction from serial sections includes aligning adjacent sections, segmenting the desired objects and constructing a computer internal model of the reconstructed object. Computational methodologies taking advantage of the parallel processing facilities of CLIP4 are presented for automating these tasks.
Key words: 3-D reconstruction, serial sections, array processors.
1. Introduction One way to gain a complete description of the interior as well as the exterior structure of a threedimensional object is by cutting it into slices whilst carefuly preserving the structural relationships of the components of the object. In biological terms, this means anatomical dissection, or sectioning the biological specimen on a microtome and reconstructing the object of interest from the set of serial sections. The process of three-dimensional reconstruction from serial sections may be divided into the following subtasks: (a) Alignment of the adjacent sections. (b) Segmentation of the regions corresponding to the intersections of the object of interest with a slicing plane from each of the section images. (c) Construction of a compact model suitable for subsequent morphological analysis and display of the 3-D object. Computational methodologies which take advantage of the parallel processing facilities provided by an array processor, CLIP4 (Duff (1976), Fountain (1982)), have been developed with the * Supported by the University of London.
aim of automating some of the tedious and timeconsuming elements of the reconstruction process. We are interested in the detailed structure of the carotid body which is a highly vascularized organ with the largest blood flow rate of any tissue in the body (Biscoe (1971), Seidl (1975), Liibbers et al. (1977), Clarke and Daly (1982)). It plays an important role in monitoring the chemical composition o f arterial blood (Po2, Pco2, pH). The aim of the present investigation is to reconstruct the total vasculature of the organ and to make an analytical study of the geometrical configuration of its vessels.
2. System configuration The present hardware configuration consists of a light microscope with a mechanical specimen stage. A rotating prism is attached between a TV camera and the barrel of the microscope along a common optical axis. The Cellular Logic Image Processor, CLIP4, consists of 96 × 96 binary processors each having 32 bits of local storage. Each of the processors can communicate with its eight nearest neighbours. Since these communication channels are in-
0167-8655/83/$3.00 © 1983, Elsevier Science Publishers B.V. (North-Holland)
497
Volume 1, Numbers 5,6
PATTERN RECOGNITION LETTERS
dividually gated, the user can tailor the shape o f the neighbourhood to suit the various local-context sensitive operations. In addition, a signal can propagate across the array of processors at the maximum clockrate o f the machine during a single instruction, giving an order of magnitude speed improvement over local operations executed iteratively. One utilization of this global propagation operation is object labelling (see Otto (1982)), which forms the basic building block of a 3-D object isolation algorithm outlined in a later section. Data and instructions are passed between CLIP4 and the PDP11/34 minicomputer via the UNIBUS1; most o f the application software is currently written in IPC (see Reynolds and Otto (1981)), an extension of the C language.
3. Aligning sections A set of serial sections, 2/am thick, of the Carotid Body of a rat were provided by St. Bartholomew's Hospital Medical College for the initial investigation (see Clarke and Daly (1982)). An area corresponding to 125/amx 125/am on a tissue section is digitized into a 96 x 96, 6-bit pixel image with a resulting pixel resolution of 1.3/am. The process of aligning two successive sections aims at approximating the spatial structure of the object in the third dimension which is lost in sectioning. In the absence of any fiducial marks, some kind of alignment criterion for two section images has to be defined. In our case, the design of the alignment procedure is based on the assumption that the corresponding vessel regions in adjacent section images, when superimposed, are close to each other and therefore only small translations and small angular rotations are needed to align them. The assumption can be met by overlaying the two section images on a TV monitor and approximately positioning one with the other under the microscope by manually adjusting the mechanical stage and the rotating prism. The computer then determines the exact movements required to bring the images into register. Local distortions within a section as result of the slicing and sample l UNIBUS is a trademark of Digital Equipment Corp. 498
July 1983
preparation process, which are found to be small, are not corrected since that would require the evaluation of a distortion map for each section. The alignment problem can be considered as one o f change detection in motion analysis, e.g. Yalamanchili et al. (1982). One image is transformed with respect to the other in such a way as to minimize the observed changes between the two images. In our case, the major vessel regions within a section image provide the necessary landmarks. They are extracted by thresholding the grey section images and comparing the resultant images. Initially, the centres of mass of the two aligning images are made to coincide; subsequent alignments proceed iteratively. During each iteration, the movements of the image are guided by a difference measure. The set difference, n, as defined by equation (I), is minimized, through successive iterations, by hill-climbing in three parameters, i.e. horizontal shifts, vertical shifts, and rotations in predetermined increments.
n = ]~ Pl(i,j) @P2(i,j)
(1)
i,j
where p l ( i , j ) and p2(i,j) are the pixel values of the two thresholded images at location (i, j ) , and @ denotes the Boolean exclusive OR operator. Experience has shown that whilst it is relatively easy to align two images roughly, immense concentration and patience is needed to bring them into exact registration manually. The advantage of the hill-climbing approach is that shifting and rotating images, though computationally expensive in a serial computer, are relatively trivial to compute in an array processor like CLIP4 which makes it a feasible and practical technique for the search of a minimum in the difference measure. Shifting an image by one pixel can be accomplished in a single instruction. Similarly, image rotation can be achieved by a controlled sequence o f image shifts in various directions (see Clarke and Ip (1982)), and does not involve matrix multiplications. A sequence of aligned section images can be obtained by repeating the process.
4. Detection of blood vessels The serial sections are prepared and stained in
Volume 1, Numbers 5,6
PATTERN RECOGNITION LETTERS
July 1983
Fig. 1. (a) A sequence of aligned tissue section images. (b) Result of vessel segmentation. (c) A sequence showing an isolated branching vascular structure. 'x' denotes the initial 'seed'. 499
Volume 1, Numbers 5,6
PATTERN RECOGNITION LETTERS
such a way that the vessels appear bright against a dark and textural background of cells and connective tissues. An example of the typical appearance of the tissue sections is shown in Figure l(a). In order to facilitate measurement and automatic tracking of vessel connectivities, it is usually helpful to reduce the grey section image into a binary picture in which ones (white) denote vessel interiors and zeros (black) denote background. Since we are dealing with three-dimensional structures and have, at our disposal, essentially three-dimensional data, it is natural that both the knowledge of the vessel region's image characteristics (grey values, shapes, etc.), and the assumption that they are continuous in the third dimension are exploited in the vessel segmentation process.
4.1. Threshold selection One method of selecting the 'optimal' threshold which minimizes the error of rejecting pixels belonging to the objects of interest is the p-tile method (see Rosenfeld and Kak (1976), p. 263). Since we do not usually know in advance the area of the image occupied by the objects, this method is rarely used for general image processing. However, in an image sequence, with a sufficiently large between-frame sampling frequency, the desired threshold for the next image in the sequence can be estimated by correlating the area occupied by the objects, in our case vessels, in the images already processed. Such a p-tile correlation method can be summarized as follows: For the nth section image, (a) calculate the percentage-area, kn_l, occupied by the vessel regions extracted from the ( n - 1)th section image, and (b) choose a threshold value, t, such that [kt -kn-11 is minimized, where k t is the percentage-area occupied pixels having grey values greater than t in the current section image. Alternatively, the selected threshold can be expressed as follows:
kt= 100× ~ Pt(i,j)/N i.j
(2)
and pick a t' such that
I k r - k , _ l l = M i n lkL-tn_ll
(3)
t
where pt(i,j) is the pixel value, i.e. 1 or 0, at 500
July 1983
location (i, j) of the section image thresholded at t; t' is the selected threshold for the current image; and N is the total number of pixels in the image.
4.2. Refinements Threshold selection based on a correlation with the segmentation result obtained previously, if used alone, would suffer the risk of allowing propagation of errors if the previous images were not segmented properly. Therefore, an 'optimal' threshold value having been selected using the method described above, it is modified according to the current image grey-level statistics. A sample of pixels is selected from the image resulting from thresholding the current section image at the 'optimal' value. This is done by the object labelling operation, i.e. regions in the thresholded image qualify as samples if they overlap one or more vessel regions in the previous image when superimposed. Using these sample regions as masks, the grey-value characteristics, e.g. mean, variance, etc., of the objects (pixels of the original grey image lying inside the masks), and background (pixels lying outside the masks) are evaluated and a better estimate for the final threshold value is obtained from these statistical parameters. Usually, no single threshold value will be sufficient to separate the objects from the background exactly. The thresholded regions are further subjected to a classification procedure during which they are assigned to one of three classes, namely: high-confidence vessel regions, query regions and artefacts. The design of the decision tree is based on knowledge of the shapes and sizes of vessel cross-sections and the gradient characteristics around the perimeters of vessel regions. Therefore, the information used for discrimination is essentially two-dimensional. The structure of the decision tree is shown schematically in Figure 2.
4.3. Three-dimensional referencing Regions labelled as artefacts are discarded immediately, whilst high-confidence ones are accepted as valid vessel regions and stored in the result image for the current section. The query regions usually correspond to tissue-tears, lightly stained intercellular gaps and small vessels. The
Volume 1, Numbers'5,6
PATTERN RECOGNITION LETTERS
occurrences of the first two types of structure are expected to be of a r a n d o m nature and, consequently, are not continuous in the third dimension. These query regions are held temporarily in a reference image which provides a means of looking ahead for corresponding vessel regions in the next few consecutive sections. Each region in the reference image is labelled with an 'age' denoting from how far back in the image sequence it originally came. Any query regions for which matches are not found within a pre-set number of following sections are rejected as artefacts. Otherwise, they are accepted as valid vessel regions and transferred f r o m the reference image to the corresponding result images. This referencing strategy aims to minimize the rejection of valid vessel regions which share m a n y c o m m o n features with other artefacts in a section image. It should be noted that a h u m a n expert, when identifying vessels on sections, frequently refers to the next few consecutive sections before or after the current section in the sequence. This
July 1983
referencing scheme is designed to model such a procedure. In essence, both two- and three-dimensional information inherent in an image sequence is used to extract the desired objects reliably from the available data. All the processes described above are implemented in parallel in the sense that all regions in an image, whenever necessary, are processed simultaneously. Results of the segmentations on an image sequence are shown in Figure
l(b).
5. Three-dimensional object isolation The result of the segmentation process is a sequence of binary images. The question to ask at this stage is which of the regions in a section are connected to other regions to form a connected 3-D object (structure). The solution of this question enables the isolation of particular 3-D objects, counting objects, and subsequently, completely
SET OF POSSIBLE VESSEL REGIONS ( THRESHOLDEOIMAGE)
SIZE FILTERING
J LARGE REGIONS (>.5 PIXELS IN DIAMETER)
SHAPE MEASURES
I SMALL REGIONS
BOUNDARY GRADIENT CHARACTERISTICS
It__ LARGE VESSELS (HIGH CONFIDENCE REGIONS)
l POSSIBLETISSUE TEARS (QUERY)
3-D CORRELATION REQUIRED
POSSIBLE SMALL VESSELS (QUERY)
CONNECTIVE TISSUES etc. (ARTEFACTS)
3-0 CORRELATION REQUIRED
Fig. 2. Classification scheme. 501
Volume 1, Numbers 5,6
PATTERN RECOGNITION LETTERS
describing the morphology of any or all of the 3-D objects in the sequence. The main problem in computing such a solution is usually the sheer amount of computation involved. A 3-D object is typically digitized into 10 5 voxels (volume elements, a 3-D analogy of pixels); to isolate or detect the object implies the process of classifying each o f these voxels as a member of the object or otherwise. In these situations, the processing power of a parallel array processor may be exploited. Using a control search strategy and the object labelling operation, once one or more 'seeds' have been defined to denote the desired object, the process searches up and down the sequence of images for voxels 'belonging' to the same object until the
July 1983
entire object is completely isolated. Depending on the complexity of the 3-D object, this process may need to be repeated several times. The object labelling operation usually involves two binary images. Regions (or points) in one o f the images, say B~, act as labels for the objects in the other image, say B 2. Propagation signals initiated by the labels in B 1 are passed to the relevant pixels in image B2, the extent o f the propagation being determined by the result of a Boolean function between the two images involved. The outcome of the operation is that objects in image B2 are retained if they overlap one or more of the labels in image B 1 when superimposed. The rest of the regions in image B 2 are discarded. This makes essential use of the global propagation functions
N-1
P.1 BLOCK
2
K*|
- -
- - - -
BLOCK BOUNDARY
BLOCK 1 SECTION
f
=
'SEED'
e
a)
K
!
! 0
b)
K+I
c)
d)
Fig. 3. (a) Original set of 3-D objects, partitioned into 2 separate blocks. (b) Result after one search of Block 1. Region in the kth plane becomes the label for the (k+ 1)th plane of Block 2. (c) Result after one search of Block 2. Regions in the (k+ 1)th plane become the labels for the kth plane of Block 1. (d) Re-entering Block 1. Result after second search of Block 1. Note that, if the block boundary is located between planes p and (p+ 1), only one search through each of the resulting blocks is needed to isolate the desired object. 502
Volume 1, N u m b e r s 5,6
July 1983
PATTERN RECOGNITION LETTERS
I I
i
SON POINTER RECORD
ii
J -I
L
PARENT POINTER
J
SON POINTERS
~ RECORD PARENT POINTER
SON POINTER RECORD PARENT POINTER
VESSEL REGION I. D
q
" 1 J
CENTROID AREA CHAIN CODE POINTERS
RECORD
PARENTPOINTER
CHAIN CODES
T a)
b) Fig. 4. (a) A n example of the internal model of a branching vessel. (b) Vessel region (node) record.
unique to C L I P 4 to speed up the signal propagation process. For large objects, due to the limited local storage within each processor, some kind of efficient search strategy is needed to minimize the search time. The problem is similar to that encountered in trying to process a large amount of data in a serial machine. The core can only hold part of the data set and the remainder is stored in some form of mass-storage device, e.g. magnetic disk or tape, which has a relatively slow access time. To this end, the set of serial section images is partitioned into several blocks 2 of images. The partial structures of the desired object are initially isolated f r o m the others within each block, and subsequently combined to form the entire object. Figure l(c) shows a sequence of the isolated vessels. Block re-entries may be needed if the object is highly convoluted. This method makes effective use of the fast memories available within 2 The n u m b e r of images in each block does not have to be the same. With careful partitioning of the data set, the search time can be further reduced.
the computer, and is particularly suitable for isolating convoluted 3-D structures, e.g. vessels having lots of recurrent branches. Figure 3 depicts diagrammatically the 3-D object isolation process.
6. Internal models It is desirable to store only the essential information required for subsequent morphological analysis and displays of the reconstructed structures. The reasons are (a) to economize on the storage space, and (b) if the vascular structures were to be analysed in a serial machine, pixel (voxel) representation or lists of contour coordinates are hardly the appropriate data structures to use. 6.1. Data structure
A connected graph seems to be a natural representation for a vascular network reconstructed f r o m serial sections. Each node of the graph represents a vessel region and the inter503
Volume 1, Numbers 5,6
PATTERN RECOGNITION LETTERS
plane connectivities of various vessel regions are represented by the arcs joining the corresponding nodes. Within a node, a set of features, e.g. centroid, area, Freeman chaincodes of the boundaries, etc., of the corresponding vessel region in the section image is recorded (see Figure 4). The result is a compact model containing quantitative and structural information of the vasculature. Current implementation using an array of the 'struct' data structure in the C language (similar to the 'record' in Pascal) typically reduces the storage requirement of a vascular structure by half. Furthermore, most of the existing graph-traversal and search algorithms, e.g. breadth-first search, depth-first search (e.g. Pavlidis (1977), p. 51), etc., can readily be adapted to provide efficient morphological analysis of the vasculature. For example, variations of the vessel diameters along all possible pathways starting from any particular point in the vascular network can be studied using a modified depth-first search method to transverse the desired routes in the graph representing the vasculature. 6.2. Presentation
The results of the reconstruction may be presented to experts for interpretation in various ways. Currently, there are three modes of display available: (a) The 'skeleton' of the vascular network can
be displayed with each 'bone' colour-coded according to the size of the corresponding vessel (see Figure 5). (b) The outlines (boundaries) of the vessel crosssections in each section are stacked and displayed with hidden lines removed (see Figure 5). (c) Stereo-pairs of the vascular structure provide a 'feel' for the complex spatial relationships of the vaious branches in the network. Algorithms for approximating the surface of a 3-D object from a series of planar contours (e.g. Fuchs et al. (1977), Christianson and Sederberg (1978), Cook et al. (1981)) can also be applied to produce shadedsurface displays of the vascular network from its graph representation.
7. Conclusion
We have described the approach taken to relieve human experts of the tedious and time-consuming tasks of reconstructing a complex vascular network from a set of serial sections. In a serial machine, the computational complexity of such reconstruction process was proportional to the number of voxels representing the objects. The use of parallel processing in CLIP4 reduces this complexity so that it becomes typically proportional to the number of serial sections. With the advent of non-invasive imaging techni-
a Fig. 5. Two modes of displaying the reconstructed vasculature.
504
July 1983
Volume 1, Numbers 5,6
PATTERN RECOGNITION LETTERS
ques, such as X-ray tomography, NMR tomography, and ultrasound tomography, planar data becomes, increasingly, a common form of data for three-dimensional objects. It is expected that the methodologies reported in this paper may be extended to these forms of planar data.
Acknowledgements The author wishes to thank Dr. M.J.B. Duff for his encouragement and support, his colleagues in the Image Processing Group for their discussions. Dr. J.A. Clarke and Prof. M. de Burgh Daly (Departments of Anatomy and Physiology, St. Bartholomew's Hospital Medical College, London) for their special biological knowledge and the supply of sample tissue sections, and Mr. S. Jones of the Department of Histopathology, St. Bartholomew's Hospital, for expert assistance in the preparation of the material.
References Biscoe, T.J. (1971). Carotid body: Structure and function. Physiological Reviews 51, 437-481. Christianson, H.N. and T.W. Sederberg (1978). Conversion of complex contour line definitions into polygonal element mosaics. Computer Graphics 12, 187-192. Clarke, J.A. and M. de B. Daly (1982). Distribution of carotid body type I ceils and periadventitial type I cells in the carotid bifurcation regions of the dog. Acta Anat. 113, 352-370.
July 1983
Clarke, K.A. and H. H.-S. Ip (1982). A parallel implementation of geometric transformations. Pattern Recognition Letters 1, 51-53. Cook, P.N. et al. (1981). Three-dimensional reconstruction from serial sections for medical applications. SPIE 3-D Machine Perception 283, 98-105. Duff, M.J.B. (1976). CLIP4: A large scale integrated circuit array parallel processor. Proc. 3rd Int. Joint Conference Pattern Recognition, November 1976, 728-733. Fountain, T.J. (1982). The CLIP4 image processing system. General Engineer 93, 267-281. Fuchs, H., Z.M. Kedem and S.P. Uselton (1977). Optimal surface reconstruction from planar contours. Comm. A C M 20, 693-702. Liibbers, D.W., L. Teckhaus and E. Seidl (1977). Capillary distances and oxygen supply to the specific tissue of the carotid body. In: H. Acker et al., eds., Chemoreceptor in the Carotid Body, 62-68. Springer, Berlin. Otto, G.P. (1982). Counting and labelling objects on CLIP 4. Report 82/11, Image Processing Group, Dept. of Physics and Astronomy, University College London. Pavlidis, T. (1977). Structural Pattern Recognition. Springer, Berlin. Reynolds, D.E. and G.P. Otto (1981). IPC user manual. Report 82/4, Image Processing Group, Dept. of Physics and Astronomy, University College London. Rosenfeld, A. and A.C. Kak (1976). Digital Picture Processing. Academic Press, New York. Seidl, E. (1975). On the morphology of the vascular system of the carotid body of cat and rabbit and its relation to the glomus type I cells. In: M.J. Purves, ed., The Peripheral Arterial Chemoreceptors. Cambridge University Press, Great Britain. Yalamanchili, S., W.N. Martin and J.K. Aggarwal (1982). Extraction of moving object descriptions via differencing. Computer Graphics and Image Processing 18, 188-201.
505
July 1983
Detection and three-dimensional reconstruction o f a vascular network from serial sections H. H.-S. I P * Image Processing Group, Department of Physics and Astronomy, University College London, London WCIE 6BT, U.K. Received 30 March 1983
Abstract: The process of three-dimensional reconstruction from serial sections includes aligning adjacent sections, segmenting the desired objects and constructing a computer internal model of the reconstructed object. Computational methodologies taking advantage of the parallel processing facilities of CLIP4 are presented for automating these tasks.
Key words: 3-D reconstruction, serial sections, array processors.
1. Introduction One way to gain a complete description of the interior as well as the exterior structure of a threedimensional object is by cutting it into slices whilst carefuly preserving the structural relationships of the components of the object. In biological terms, this means anatomical dissection, or sectioning the biological specimen on a microtome and reconstructing the object of interest from the set of serial sections. The process of three-dimensional reconstruction from serial sections may be divided into the following subtasks: (a) Alignment of the adjacent sections. (b) Segmentation of the regions corresponding to the intersections of the object of interest with a slicing plane from each of the section images. (c) Construction of a compact model suitable for subsequent morphological analysis and display of the 3-D object. Computational methodologies which take advantage of the parallel processing facilities provided by an array processor, CLIP4 (Duff (1976), Fountain (1982)), have been developed with the * Supported by the University of London.
aim of automating some of the tedious and timeconsuming elements of the reconstruction process. We are interested in the detailed structure of the carotid body which is a highly vascularized organ with the largest blood flow rate of any tissue in the body (Biscoe (1971), Seidl (1975), Liibbers et al. (1977), Clarke and Daly (1982)). It plays an important role in monitoring the chemical composition o f arterial blood (Po2, Pco2, pH). The aim of the present investigation is to reconstruct the total vasculature of the organ and to make an analytical study of the geometrical configuration of its vessels.
2. System configuration The present hardware configuration consists of a light microscope with a mechanical specimen stage. A rotating prism is attached between a TV camera and the barrel of the microscope along a common optical axis. The Cellular Logic Image Processor, CLIP4, consists of 96 × 96 binary processors each having 32 bits of local storage. Each of the processors can communicate with its eight nearest neighbours. Since these communication channels are in-
0167-8655/83/$3.00 © 1983, Elsevier Science Publishers B.V. (North-Holland)
497
Volume 1, Numbers 5,6
PATTERN RECOGNITION LETTERS
dividually gated, the user can tailor the shape o f the neighbourhood to suit the various local-context sensitive operations. In addition, a signal can propagate across the array of processors at the maximum clockrate o f the machine during a single instruction, giving an order of magnitude speed improvement over local operations executed iteratively. One utilization of this global propagation operation is object labelling (see Otto (1982)), which forms the basic building block of a 3-D object isolation algorithm outlined in a later section. Data and instructions are passed between CLIP4 and the PDP11/34 minicomputer via the UNIBUS1; most o f the application software is currently written in IPC (see Reynolds and Otto (1981)), an extension of the C language.
3. Aligning sections A set of serial sections, 2/am thick, of the Carotid Body of a rat were provided by St. Bartholomew's Hospital Medical College for the initial investigation (see Clarke and Daly (1982)). An area corresponding to 125/amx 125/am on a tissue section is digitized into a 96 x 96, 6-bit pixel image with a resulting pixel resolution of 1.3/am. The process of aligning two successive sections aims at approximating the spatial structure of the object in the third dimension which is lost in sectioning. In the absence of any fiducial marks, some kind of alignment criterion for two section images has to be defined. In our case, the design of the alignment procedure is based on the assumption that the corresponding vessel regions in adjacent section images, when superimposed, are close to each other and therefore only small translations and small angular rotations are needed to align them. The assumption can be met by overlaying the two section images on a TV monitor and approximately positioning one with the other under the microscope by manually adjusting the mechanical stage and the rotating prism. The computer then determines the exact movements required to bring the images into register. Local distortions within a section as result of the slicing and sample l UNIBUS is a trademark of Digital Equipment Corp. 498
July 1983
preparation process, which are found to be small, are not corrected since that would require the evaluation of a distortion map for each section. The alignment problem can be considered as one o f change detection in motion analysis, e.g. Yalamanchili et al. (1982). One image is transformed with respect to the other in such a way as to minimize the observed changes between the two images. In our case, the major vessel regions within a section image provide the necessary landmarks. They are extracted by thresholding the grey section images and comparing the resultant images. Initially, the centres of mass of the two aligning images are made to coincide; subsequent alignments proceed iteratively. During each iteration, the movements of the image are guided by a difference measure. The set difference, n, as defined by equation (I), is minimized, through successive iterations, by hill-climbing in three parameters, i.e. horizontal shifts, vertical shifts, and rotations in predetermined increments.
n = ]~ Pl(i,j) @P2(i,j)
(1)
i,j
where p l ( i , j ) and p2(i,j) are the pixel values of the two thresholded images at location (i, j ) , and @ denotes the Boolean exclusive OR operator. Experience has shown that whilst it is relatively easy to align two images roughly, immense concentration and patience is needed to bring them into exact registration manually. The advantage of the hill-climbing approach is that shifting and rotating images, though computationally expensive in a serial computer, are relatively trivial to compute in an array processor like CLIP4 which makes it a feasible and practical technique for the search of a minimum in the difference measure. Shifting an image by one pixel can be accomplished in a single instruction. Similarly, image rotation can be achieved by a controlled sequence o f image shifts in various directions (see Clarke and Ip (1982)), and does not involve matrix multiplications. A sequence of aligned section images can be obtained by repeating the process.
4. Detection of blood vessels The serial sections are prepared and stained in
Volume 1, Numbers 5,6
PATTERN RECOGNITION LETTERS
July 1983
Fig. 1. (a) A sequence of aligned tissue section images. (b) Result of vessel segmentation. (c) A sequence showing an isolated branching vascular structure. 'x' denotes the initial 'seed'. 499
Volume 1, Numbers 5,6
PATTERN RECOGNITION LETTERS
such a way that the vessels appear bright against a dark and textural background of cells and connective tissues. An example of the typical appearance of the tissue sections is shown in Figure l(a). In order to facilitate measurement and automatic tracking of vessel connectivities, it is usually helpful to reduce the grey section image into a binary picture in which ones (white) denote vessel interiors and zeros (black) denote background. Since we are dealing with three-dimensional structures and have, at our disposal, essentially three-dimensional data, it is natural that both the knowledge of the vessel region's image characteristics (grey values, shapes, etc.), and the assumption that they are continuous in the third dimension are exploited in the vessel segmentation process.
4.1. Threshold selection One method of selecting the 'optimal' threshold which minimizes the error of rejecting pixels belonging to the objects of interest is the p-tile method (see Rosenfeld and Kak (1976), p. 263). Since we do not usually know in advance the area of the image occupied by the objects, this method is rarely used for general image processing. However, in an image sequence, with a sufficiently large between-frame sampling frequency, the desired threshold for the next image in the sequence can be estimated by correlating the area occupied by the objects, in our case vessels, in the images already processed. Such a p-tile correlation method can be summarized as follows: For the nth section image, (a) calculate the percentage-area, kn_l, occupied by the vessel regions extracted from the ( n - 1)th section image, and (b) choose a threshold value, t, such that [kt -kn-11 is minimized, where k t is the percentage-area occupied pixels having grey values greater than t in the current section image. Alternatively, the selected threshold can be expressed as follows:
kt= 100× ~ Pt(i,j)/N i.j
(2)
and pick a t' such that
I k r - k , _ l l = M i n lkL-tn_ll
(3)
t
where pt(i,j) is the pixel value, i.e. 1 or 0, at 500
July 1983
location (i, j) of the section image thresholded at t; t' is the selected threshold for the current image; and N is the total number of pixels in the image.
4.2. Refinements Threshold selection based on a correlation with the segmentation result obtained previously, if used alone, would suffer the risk of allowing propagation of errors if the previous images were not segmented properly. Therefore, an 'optimal' threshold value having been selected using the method described above, it is modified according to the current image grey-level statistics. A sample of pixels is selected from the image resulting from thresholding the current section image at the 'optimal' value. This is done by the object labelling operation, i.e. regions in the thresholded image qualify as samples if they overlap one or more vessel regions in the previous image when superimposed. Using these sample regions as masks, the grey-value characteristics, e.g. mean, variance, etc., of the objects (pixels of the original grey image lying inside the masks), and background (pixels lying outside the masks) are evaluated and a better estimate for the final threshold value is obtained from these statistical parameters. Usually, no single threshold value will be sufficient to separate the objects from the background exactly. The thresholded regions are further subjected to a classification procedure during which they are assigned to one of three classes, namely: high-confidence vessel regions, query regions and artefacts. The design of the decision tree is based on knowledge of the shapes and sizes of vessel cross-sections and the gradient characteristics around the perimeters of vessel regions. Therefore, the information used for discrimination is essentially two-dimensional. The structure of the decision tree is shown schematically in Figure 2.
4.3. Three-dimensional referencing Regions labelled as artefacts are discarded immediately, whilst high-confidence ones are accepted as valid vessel regions and stored in the result image for the current section. The query regions usually correspond to tissue-tears, lightly stained intercellular gaps and small vessels. The
Volume 1, Numbers'5,6
PATTERN RECOGNITION LETTERS
occurrences of the first two types of structure are expected to be of a r a n d o m nature and, consequently, are not continuous in the third dimension. These query regions are held temporarily in a reference image which provides a means of looking ahead for corresponding vessel regions in the next few consecutive sections. Each region in the reference image is labelled with an 'age' denoting from how far back in the image sequence it originally came. Any query regions for which matches are not found within a pre-set number of following sections are rejected as artefacts. Otherwise, they are accepted as valid vessel regions and transferred f r o m the reference image to the corresponding result images. This referencing strategy aims to minimize the rejection of valid vessel regions which share m a n y c o m m o n features with other artefacts in a section image. It should be noted that a h u m a n expert, when identifying vessels on sections, frequently refers to the next few consecutive sections before or after the current section in the sequence. This
July 1983
referencing scheme is designed to model such a procedure. In essence, both two- and three-dimensional information inherent in an image sequence is used to extract the desired objects reliably from the available data. All the processes described above are implemented in parallel in the sense that all regions in an image, whenever necessary, are processed simultaneously. Results of the segmentations on an image sequence are shown in Figure
l(b).
5. Three-dimensional object isolation The result of the segmentation process is a sequence of binary images. The question to ask at this stage is which of the regions in a section are connected to other regions to form a connected 3-D object (structure). The solution of this question enables the isolation of particular 3-D objects, counting objects, and subsequently, completely
SET OF POSSIBLE VESSEL REGIONS ( THRESHOLDEOIMAGE)
SIZE FILTERING
J LARGE REGIONS (>.5 PIXELS IN DIAMETER)
SHAPE MEASURES
I SMALL REGIONS
BOUNDARY GRADIENT CHARACTERISTICS
It__ LARGE VESSELS (HIGH CONFIDENCE REGIONS)
l POSSIBLETISSUE TEARS (QUERY)
3-D CORRELATION REQUIRED
POSSIBLE SMALL VESSELS (QUERY)
CONNECTIVE TISSUES etc. (ARTEFACTS)
3-0 CORRELATION REQUIRED
Fig. 2. Classification scheme. 501
Volume 1, Numbers 5,6
PATTERN RECOGNITION LETTERS
describing the morphology of any or all of the 3-D objects in the sequence. The main problem in computing such a solution is usually the sheer amount of computation involved. A 3-D object is typically digitized into 10 5 voxels (volume elements, a 3-D analogy of pixels); to isolate or detect the object implies the process of classifying each o f these voxels as a member of the object or otherwise. In these situations, the processing power of a parallel array processor may be exploited. Using a control search strategy and the object labelling operation, once one or more 'seeds' have been defined to denote the desired object, the process searches up and down the sequence of images for voxels 'belonging' to the same object until the
July 1983
entire object is completely isolated. Depending on the complexity of the 3-D object, this process may need to be repeated several times. The object labelling operation usually involves two binary images. Regions (or points) in one o f the images, say B~, act as labels for the objects in the other image, say B 2. Propagation signals initiated by the labels in B 1 are passed to the relevant pixels in image B2, the extent o f the propagation being determined by the result of a Boolean function between the two images involved. The outcome of the operation is that objects in image B2 are retained if they overlap one or more of the labels in image B 1 when superimposed. The rest of the regions in image B 2 are discarded. This makes essential use of the global propagation functions
N-1
P.1 BLOCK
2
K*|
- -
- - - -
BLOCK BOUNDARY
BLOCK 1 SECTION
f
=
'SEED'
e
a)
K
!
! 0
b)
K+I
c)
d)
Fig. 3. (a) Original set of 3-D objects, partitioned into 2 separate blocks. (b) Result after one search of Block 1. Region in the kth plane becomes the label for the (k+ 1)th plane of Block 2. (c) Result after one search of Block 2. Regions in the (k+ 1)th plane become the labels for the kth plane of Block 1. (d) Re-entering Block 1. Result after second search of Block 1. Note that, if the block boundary is located between planes p and (p+ 1), only one search through each of the resulting blocks is needed to isolate the desired object. 502
Volume 1, N u m b e r s 5,6
July 1983
PATTERN RECOGNITION LETTERS
I I
i
SON POINTER RECORD
ii
J -I
L
PARENT POINTER
J
SON POINTERS
~ RECORD PARENT POINTER
SON POINTER RECORD PARENT POINTER
VESSEL REGION I. D
q
" 1 J
CENTROID AREA CHAIN CODE POINTERS
RECORD
PARENTPOINTER
CHAIN CODES
T a)
b) Fig. 4. (a) A n example of the internal model of a branching vessel. (b) Vessel region (node) record.
unique to C L I P 4 to speed up the signal propagation process. For large objects, due to the limited local storage within each processor, some kind of efficient search strategy is needed to minimize the search time. The problem is similar to that encountered in trying to process a large amount of data in a serial machine. The core can only hold part of the data set and the remainder is stored in some form of mass-storage device, e.g. magnetic disk or tape, which has a relatively slow access time. To this end, the set of serial section images is partitioned into several blocks 2 of images. The partial structures of the desired object are initially isolated f r o m the others within each block, and subsequently combined to form the entire object. Figure l(c) shows a sequence of the isolated vessels. Block re-entries may be needed if the object is highly convoluted. This method makes effective use of the fast memories available within 2 The n u m b e r of images in each block does not have to be the same. With careful partitioning of the data set, the search time can be further reduced.
the computer, and is particularly suitable for isolating convoluted 3-D structures, e.g. vessels having lots of recurrent branches. Figure 3 depicts diagrammatically the 3-D object isolation process.
6. Internal models It is desirable to store only the essential information required for subsequent morphological analysis and displays of the reconstructed structures. The reasons are (a) to economize on the storage space, and (b) if the vascular structures were to be analysed in a serial machine, pixel (voxel) representation or lists of contour coordinates are hardly the appropriate data structures to use. 6.1. Data structure
A connected graph seems to be a natural representation for a vascular network reconstructed f r o m serial sections. Each node of the graph represents a vessel region and the inter503
Volume 1, Numbers 5,6
PATTERN RECOGNITION LETTERS
plane connectivities of various vessel regions are represented by the arcs joining the corresponding nodes. Within a node, a set of features, e.g. centroid, area, Freeman chaincodes of the boundaries, etc., of the corresponding vessel region in the section image is recorded (see Figure 4). The result is a compact model containing quantitative and structural information of the vasculature. Current implementation using an array of the 'struct' data structure in the C language (similar to the 'record' in Pascal) typically reduces the storage requirement of a vascular structure by half. Furthermore, most of the existing graph-traversal and search algorithms, e.g. breadth-first search, depth-first search (e.g. Pavlidis (1977), p. 51), etc., can readily be adapted to provide efficient morphological analysis of the vasculature. For example, variations of the vessel diameters along all possible pathways starting from any particular point in the vascular network can be studied using a modified depth-first search method to transverse the desired routes in the graph representing the vasculature. 6.2. Presentation
The results of the reconstruction may be presented to experts for interpretation in various ways. Currently, there are three modes of display available: (a) The 'skeleton' of the vascular network can
be displayed with each 'bone' colour-coded according to the size of the corresponding vessel (see Figure 5). (b) The outlines (boundaries) of the vessel crosssections in each section are stacked and displayed with hidden lines removed (see Figure 5). (c) Stereo-pairs of the vascular structure provide a 'feel' for the complex spatial relationships of the vaious branches in the network. Algorithms for approximating the surface of a 3-D object from a series of planar contours (e.g. Fuchs et al. (1977), Christianson and Sederberg (1978), Cook et al. (1981)) can also be applied to produce shadedsurface displays of the vascular network from its graph representation.
7. Conclusion
We have described the approach taken to relieve human experts of the tedious and time-consuming tasks of reconstructing a complex vascular network from a set of serial sections. In a serial machine, the computational complexity of such reconstruction process was proportional to the number of voxels representing the objects. The use of parallel processing in CLIP4 reduces this complexity so that it becomes typically proportional to the number of serial sections. With the advent of non-invasive imaging techni-
a Fig. 5. Two modes of displaying the reconstructed vasculature.
504
July 1983
Volume 1, Numbers 5,6
PATTERN RECOGNITION LETTERS
ques, such as X-ray tomography, NMR tomography, and ultrasound tomography, planar data becomes, increasingly, a common form of data for three-dimensional objects. It is expected that the methodologies reported in this paper may be extended to these forms of planar data.
Acknowledgements The author wishes to thank Dr. M.J.B. Duff for his encouragement and support, his colleagues in the Image Processing Group for their discussions. Dr. J.A. Clarke and Prof. M. de Burgh Daly (Departments of Anatomy and Physiology, St. Bartholomew's Hospital Medical College, London) for their special biological knowledge and the supply of sample tissue sections, and Mr. S. Jones of the Department of Histopathology, St. Bartholomew's Hospital, for expert assistance in the preparation of the material.
References Biscoe, T.J. (1971). Carotid body: Structure and function. Physiological Reviews 51, 437-481. Christianson, H.N. and T.W. Sederberg (1978). Conversion of complex contour line definitions into polygonal element mosaics. Computer Graphics 12, 187-192. Clarke, J.A. and M. de B. Daly (1982). Distribution of carotid body type I ceils and periadventitial type I cells in the carotid bifurcation regions of the dog. Acta Anat. 113, 352-370.
July 1983
Clarke, K.A. and H. H.-S. Ip (1982). A parallel implementation of geometric transformations. Pattern Recognition Letters 1, 51-53. Cook, P.N. et al. (1981). Three-dimensional reconstruction from serial sections for medical applications. SPIE 3-D Machine Perception 283, 98-105. Duff, M.J.B. (1976). CLIP4: A large scale integrated circuit array parallel processor. Proc. 3rd Int. Joint Conference Pattern Recognition, November 1976, 728-733. Fountain, T.J. (1982). The CLIP4 image processing system. General Engineer 93, 267-281. Fuchs, H., Z.M. Kedem and S.P. Uselton (1977). Optimal surface reconstruction from planar contours. Comm. A C M 20, 693-702. Liibbers, D.W., L. Teckhaus and E. Seidl (1977). Capillary distances and oxygen supply to the specific tissue of the carotid body. In: H. Acker et al., eds., Chemoreceptor in the Carotid Body, 62-68. Springer, Berlin. Otto, G.P. (1982). Counting and labelling objects on CLIP 4. Report 82/11, Image Processing Group, Dept. of Physics and Astronomy, University College London. Pavlidis, T. (1977). Structural Pattern Recognition. Springer, Berlin. Reynolds, D.E. and G.P. Otto (1981). IPC user manual. Report 82/4, Image Processing Group, Dept. of Physics and Astronomy, University College London. Rosenfeld, A. and A.C. Kak (1976). Digital Picture Processing. Academic Press, New York. Seidl, E. (1975). On the morphology of the vascular system of the carotid body of cat and rabbit and its relation to the glomus type I cells. In: M.J. Purves, ed., The Peripheral Arterial Chemoreceptors. Cambridge University Press, Great Britain. Yalamanchili, S., W.N. Martin and J.K. Aggarwal (1982). Extraction of moving object descriptions via differencing. Computer Graphics and Image Processing 18, 188-201.
505