科学文献

header for SPIE usePoint-based warping with optimized weighting factorsof displacement vectorsRainer Pielot a, Michael Scholz b, Klaus Obermayer b, Eckart D. Gundelfinger a and Andreas Hess aa Leibniz Institute for Neurobiology, Brenneckestr. 6,D-39118 Magdeburg, Germanyb TU Berlin, Department of Computer Science, Sekr. FR 2-1, Franklinstr. 28/29,D-10587 Berlin, GermanyABSTRACTThe accurate comparison of inter-individual 3D image brain datasets requires non-affine transformation techniques (warping) to reduce geometric variations. Constrained by the biological prerequisites we use in this study a landmark-based warping method with weighted sums of displacement vectors, which is enhanced by an optimization process. Furthermore, we investigate fast automatic procedures for determining landmarks to improve the practicability of 3D warping. This combined approach was tested on 3D autoradiographs of Gerbil brains. The autoradiographs were obtained after injecting a non-metabolizable radioactive glucose derivative into the Gerbil thereby visualizing neuronal activity in the brain. Afterwards the brain was processed with standard autoradiographical methods. The landmark-generator computes corresponding reference points simultaneously within a given number of datasets by Monte-Carlo-techniques. The warping function is a distance weighted exponential function with a landmark-specific weighting factor. These weighting factors are optimized by a computational evolution strategy. The warping quality is quantified by several coefficients (correlation coefficient, overlap-index, and registration error). The described approach combines a highly suitable procedure to automatically detect landmarks in autoradiographical brain images and an enhanced point-based warping technique, optimizing the local weighting factors. This optimization process significantly improves the similarity between the warped and the target dataset.Keywords:Warping, Brain Imaging, Registration, Landmark Detection1.INTRODUCTIONGeometric transformation techniques play an important role for matching information of individual brain image datasets against a reference system to reduce inter-individual variations1. Such transformations consist in the simplest affine case of repositioning by translation and rotation of the image dataset relative to the reference system. Additional transforming the shape of the brain requires non-affine transformations (warping)1-3. Implementations of warping algorithms can be based on low-level information, such as intensity distributions, or on high-level information like structural information, such as corresponding points or contours of important anatomic elements. Intensity-based warping algorithms maximize local intensity similarities between the datasets; the definition of the spatial relationships can be based e.g. local cross correlation4 or squared differences of pixel intensities5. Since our interest is the inter-individual local spatial distribution of grayvalue differences instead of subject-specific shape variations intensity-based warping is not appropriate for our kind of biological analysis. Another possibility to match individual image datasets is to find corresponding reference points (landmarks) and touse them for calculation of the transformation (point-based warping). In this study, we use a distance-weighted method of warping with weighted sums of displacement vectors6-7, which is extended by an iterative optimization process8.The manual setting of such reference points between large 3D image datasets depends on prior knowledge of the expert and is very time-consuming. To overcome these problems several approaches for an automatic extraction of landmarks by incorporating knowledge about anatomic features were developed9-10. These algorithms utilize anatomic structures such as sulcal-ventricular curves and gyral crests1 as matching criteria. We approach this problem by investigating fast automatic procedures for determing landmarks, which are independent of prior anatomic knowledge in our model.To evaluate the results, three different similarity functions were applied to the datasets before and after matching against a reference template. The mean linear cross correlation coefficient measures the intensity similarity between all datasets according to correspondent grayvalue information. The overlap index describes the geometric shape correspondence of the datasets. The registration error utilizes structural information in terms of characteristic brain description points to measure the similarity.The combined approach was tested on reconstructed 3D autoradiographs of brains of Mongolian Gerbils (Meriones unguiculatus).This paper is organized as follows: The next section describes the preparation of the brains and the digitization of the autoradiographs. The preprocessing of the image-stacks (reconstruction) by image alignment is outlined in Section 3. Sections 5 and 6 describe the detection of the landmarks and the subsequent warping procedure. Section 6 describes the quality measures, which evaluate the warping results. The results are presented in Section 7 and discussed in Section 8.2.TISSUE-PREPARATIONThe autoradiographs were obtained after intraperitoneal injection of radioactive non-metabolizable 2-fluoro-2-deoxyglucose11 (2FDG) into gerbils12. This method visualizes brain activity due to locally increased excitatory neuronal metabolic activity. The gerbils were acoustically stimulated (pure tones at 70 db) with a stimulating frequency of 1 and 8 kHz. After 45 min the animals were sacrificed and the brains were removed. Then the brains were sectioned on a cryostat and after drying at 60°C autoradiographed with a KODAK NMB X-ray film. After 2 weeks exposure, the films were developed and digitized with a CCD camera (768*512 pixel, 8 bit/pixel). All animal experiments were performed in accordance with the guidelines and regulations defined by the German Animal Welfare Act. Permission was obtained from the Regierungspraesidium Dessau.3.IMAGE PREPROCESSINGReconstructions of image-stacks of physically sectioned slices are performed in a first step by image alignment. For this purpose, each slice is rotated and translated in order to maximize the spatial relations between adjacent slices. In this study we use a combinational reconstruction method13, consisting of principal axes alignment, which is fast automatic technique and a cross-correlation-based method, resulting in a high reconstruction quality14.The intensity of autoradiographs depends on various parameters like the position of the brain slice on the X-ray film and exposure time. Therefore, the images within one image stack may have a different grayvalue distribution and thus resulting in an inhomogeneous data volume. To solve this problem, the datasets were z-transformed to normalize the distribution of grayvalues:σµ-X=Z(1) X and Z are grayvalues of the voxel, µ the average grayvalue and σ the standard deviation.NDMARK DETECTIONIn this study we use an approach to define corresponding reference points, which is based on Monte-Carlo-techniques to detect edges of anatomical structures. In our model, an edge is defined by the grayvalue transition of adjacent voxels, which exceeds a certain threshold value. This value strongly depends on the level of noise and contrast in the image. A three-dimensional grid is placed on each dataset to define non-overlapping subvolumes. Within each subvolume random movements of a single point perform the search for edges. This procedure stops within each a certain subvolume, if an edge is found. Subvolumes in which no edge is found after a maximum number of iterations do not contribute a reference point. Only if a given subvolume contributes a reference point in all datasets, a reference point is generated for this subvolume in the template by averaging the positions of all corresponding reference points of the individual datasets. Consequently, the spatial difference of an individual reference point and the corresponding one in the reference template defines a displacement vector for the warping to be performed.Monte-Carlo-techniques are generally based on random numbers. In this study we use a pseudo random number generator, which is initialized with a certain seed number. As an approach to minimize the influence of the seed number, the final position of each reference point is the average of 10 runs. The pseudo random number generator, used in this study, is taken from the tested standard libraries (MIPSpro C-compiler 7.1). The quality of the pseudo random number generator is tested negatively for autocorrelation.In order to detect different displacement sizes, we use a multiscaling procedure. Two differently sized grids are used successively: First a coarse grid with large first-order subvolumes (1920) is used to detect large distortions (Fig. 1a). Afterwards, each large subvolume is divided into four smaller second-order subvolumes (7680) (Fig. 1b). If the search procedure within the four second-order subvolumes detects more than one surface in the corresponding first-order subvolume, the displacement vector in the latter one is eliminated to avoid non-corresponding surfaces, otherwise both are taken.The Monte-Carlo search was performed in each dataset separately. As mentioned above, the reference points of the template for a given subvolume are defined as the mean positions of the reference points of each individual dataset. Therefore each individual landmark couple is build up by a reference point of the individual dataset and the corresponding one of the reference template. The subsequent warping procedure matches each dataset with the information of all individual landmarks to the template.a bFigure 1: Grids of non-overlapping subvolumes. At the beginning of the Monte-Carlo search each point is located in the center of selected subvolumes. a) Coarse grid (1920 subvolumes). b) Fine grid (7680 subvolumes). For a better illustration four different anatomical structures were manually segmented (septum (SE), striatum (ST), thalamus (TH), andhippocampus (H)) and then surface rendered including the outline of the brain (Visualization software: A MIRA (Indeed GmbH)).5.WARPINGAs described above, the presented warping method is based on displacement vectors, which are defined by the spatial difference of corresponding reference points (landmarks) in the individual datasets and in the reference template. The transformation function determines the displacement of each voxel by the weighted sum of all displacement vectors.The distribution of the generated landmarks across the 3D dataset is very irregular. Simple distance-weighted warping methods determine the transformation only dependent on the absolute distance of the reference points to a given voxel. Therefore landmark-free regions are transformed with the information of the nearest available landmark, which can be localized in another structure. On the other side, there may exist regions with many landmarks lying close together describing even different transformation directions. The transformation of such regions is more or less determined by an average of these displacement vectors. One approach to overcome this problem is to determine each displacement not only by the spatial difference of a given voxel to each landmark but also by a landmark-specific weighting factor8. The weighting factor determines the influence of that specific landmark on the transformation. At the beginning all weighting factors have identical values, so that the displacement of a given voxel is mainly influenced by the closest landmarks.The following transformation function T(P) describes the transformation of the voxel P=(x,y,z) based on the reference points {i1L}={(x i,y i,z i)|0<=i<=M-1} in the source dataset and the corresponding points }{0i L={(u i,v i,w i)|0<=i<=M-1} in the reference template:T(P)=P+wi(P)(Lii=0M-1∑−L i1)wi(P)i=0M-1∑(2)The weighting function w i(P) consists of the landmark-specific weighting factor βi and the distance function d i(P) of the point P to the reference point {i1L}:wi=e−βi d i(P)(3) The distance function d i(P) bases on the city-block metric to reduce computational efforts by avoiding the computation of the square root:d i(P)=|(x-x i)|+ |(y-y i)|+ |(z-z i)|(4) Calculation of the Euclidian distance instead of the city-block distance does not improve the warping results (data not shown).As outlined above it is essential to select the weighting factors optimally due to their strong influence on the warping result. To optimize these weighting factors, we used a computational evolution strategy, which consists of the following steps: At the beginning, a population of identical sets of weighting factors is generated. The fitness of each set is determined by warping of the dataset to the template and then by calculating a registration error to quantify the similarity. The registration error is the sum of the distance values between all so-called registration points of the warped dataset and the corresponding ones of the reference template. Registration points are a small number of manually defined points at characteristic positions in all datasets. In this study, the registration points are defined by one expert. Here, we use 28 registration points, which are located on structures which can be easily identified (Fig. 2).Figure 2: Registration points. For a better illustration only one slice with registration points, for which the z-coordinates correspond ± 1 slice to the z-level of this slice, is shown here (16 out of 28 registration points in the 3D dataset are placed here).For the reference template the registration points are generated by averaging the coordinates of the registration points of the individual datasets around a given position. This registration error has to be minimized by the optimization procedure.Therefore after each generation, the set with the lowest registration error is duplicated and replaces the worst set. Then the weighting factors of all sets other than the fittest set are mutated randomly with normal distributed step sizes. This procedure prevents the optimization to run into deadlocks. These steps are iterated for a constant number of generations,which was selected to obtain a stable fitness. In this study, the optimization was realized with a population of 500 sets of local weighting factors and ended after 1,000 generation cycles.6. SIMILARITY MEASURESThe quality of the warping is determined by comparison of the mean linear cross correlation coefficient, the overlap index and the registration error between the datasets before, after warping and after optimized warping to the reference template.The mean linear cross correlation coefficient measures the intensity similarity between two datasets according to correspondent grayvalue information. The maximum value of 1.0 is reached, when the highest grayvalues exist at the same location in the two 3D datasets. A normalized image cross correlation c is defined as 14:()()∑∑∑∑∑∑∑∑∑=x y z z y x T x y z z y x W x y zz y x T z y x W I I I I c 2),,(2),,(),,(),,((5)I W(x,y,z) and I T(x,y,z) describe the grayvalues of the voxels of the datasets W and T at the position (x,y,z ). The summation is only performed in the area of overlap between the datasets.The overlap index o describes the extent of the volume of overlap, which quantifies the geometric correspondence of the two datasets to be compared. The overlap index is defined as:T W WTN N N o +=*2(6)The overlap volume between the datasets W and T is described by N WT . N W , N T are the overlap volumes of the datasets W and T , respectively. Values of o are between 0 (no overlap) and 1.0 (perfect overlap).The determination of the registration error is already described in the previous section. The registration error reaches 0when the positions of the registration points in the transformed datasets become identical (perfect match). As mentioned before, the optimization procedure is based on this similarity function.7. RESULTS Figure 3 shows as an example a surface view of an individual brain dataset together with defined landmarks (M=1363).Figure 3: Surface view of an individual brain dataset together with 1363 generated landmarks. The small lines at the points describe the displacement vectors. As in Fig. 1, four different regions (septum (SE), striatum (ST), thalamus (TH), and hippocampus (H)) were segmented manually and surface rendered together with the outer surface of the gerbil brain by the visualization software A MIRA (Indeed GmbH). The cerebellum (CE) was not segmented.Most of the landmarks are located on the surface of the structures. Some landmarks are defined within the region of the cerebellum, which is not shown as a segmented structure in this figure. As described above, the positions of the individual reference points are averaged in order to generate a reference template. Then all datasets are warped to this template. Figure 4 shows as an example a slice of a dataset before transformation (Fig. 4a) and after transformation (Fig. 4b). The contour of the corresponding slice of another dataset is shown to illustrate the increase in overlap after warping.a bFigure 4: a) Slice of a dataset before transformation b) same slice after transformation. Black contour: Outline of the corresponding slice of another dataset. The arrows point to positions where a better overlap is achieved after warping. Tab. 1 shows the averaged results of the similarity measure (3 datasets).Table 1: Results of similarity measures (averaged values of all brain datasets).All calculations were performed on a SGI Origin 200 (MIPS R10000/180 MHz). The landmark generation required 7,228 s / dataset computation time, the warping procedure including the optimization required 13,027 s / dataset computation time (averaged values, single CPU).8.DISCUSSIONIn our analysis the first step of an accurate comparison of inter-individual brain datasets is to define individual reference points. Then the reference template is generated by averaging the positions of the corresponding individual reference points. The landmarks provide the geometrical information for the subsequent warping of the individual brain datasets to thistemplate. Finally the increase of similarity between the transformed objects is quantified. The presented approach combines a fast fully automatic procedure to detect landmarks and a point-based warping technique, which optimizes the influence of displacement vectors.The presented landmark generator is a highly suitable tool for detecting landmarks simultaneously in large 3D brain datasets without incorporating any prior knowledge. The average positions of the landmarks are successfully used to generate a reference template.The optimization of the weighting factors improves the similarity of the warped datasets and hence decreases the inter-individual differences in all cases tested. Similarity measures between the warped objects show the increase of similarity after warping. Optimizing the local influence of each landmark results in an accurate warping with very good local adjustments, which justifies the relatively high computational effort of the optimization procedure.For biological questions demanding quantitative evaluation of grayvalues the above described algorithms proved to be a useful tool to make complex biological structures inter-individually more comparable and facilitate a quantitative group comparison.ACKNOWLEDGEMENTSThis work was supported by the Bundesministerium für Bildung und Forschung (BMBF) of the Federal Republic of Germany (AZ 10963, Forschungsverbund “Virtual Brain”)REFERENCES1. A.W. Toga, Brain Warping, Academic Press, San Diego, 19992.G. Wolberg, Digital Image Warping, IEEE Computer Society Press, 19903. A.W. Toga, “Visualization and warping of multimodality brain imagery,” in Functional Neuroimaging, R.W. Thatcher,M. Hallet, T. Zeffiro, E.R. John, and M. Huerta, eds., pp. 171-180. Academic Press, San Diego, 19944.R. Bajcsy, S. Kovacic, “Multiresolution elastic matching,” Computer Vision, Graphics and Image Processing, 46, pp.1-21, 19895.G.E. Christensen, R.D. Rabbit and M.I. Miller, “Deformable templates using large deformation kinematics,” IEEETransactions on Image Processing, 5(10), pp. 1435-1447, 19966. D. Shepard, “A two-dimensional interpolation function for irregularly-spaced data,” Proc. 1968 ACM NationalConference, pp. 517-524, 19687.R. Franke and G. Nielson, “Smooth interpolation of large sets of scattered data,” Int J for Numerical Methods inEngineering, 15, pp.1691-1704, 19808.R. Pielot, M. Scholz, K. Obermayer, E.D. Gundelfinger and A. Hess, “Optimiertes Warping durch gewichtete Summenvon Verschiebungsvektoren – eine neue Methode zur Reduktion von interindividuellen Variabilitäten von Hirndaten,”in Bildverarbeitung für die Medizin 1999, H. Evers, G. Glombitza, T. Lehmann and H.-P. Meinzer, eds., pp. 417-421, 19999.J. Declerck, G. Subsol, J.-P. Thirion and N. Ayache, “Automatic retrieval of anatomical structures in three-dimensionalmedical images,” INRIA Technical Report, p. 2485, 199510.C. Davatzikos, “Spatial normalization of three-dimensional brain images using deformable models,” J Comp AssistedTomography, 20(4) pp. 656-665, 199611.L. Sokoloff, M. Reivich, C. Kennedy, M.H. Des Rosiers, C.S. Patlak, K.D. Pettigrew, O. Sakurada and M. Shinohara,“The (14C)deoxyglucose method for the measurement of local cerebral glucose utilization: Theory, procedure, and normal values in the concious and anaesthetized albino rat,” J Neurochem, 28, pp. 897-916, 197712.H. Scheich, P Heil, G. Langner, “Functional organization of auditory cortex in the mongolian gerbil (Merionesunguiculatus) II. Tonotopic 2-deoxyglucose,” European J of Neuroscience, 5, pp. 898-914, 199313.L.S. Hibbard and R.A. Hawkins, “Objective image alignment for three-dimensional reconstruction of digitalautoradiograms,” J Neurosci Methods, 26, pp. 55-74, 198814.A. Hess, K. Lohmann, E.D. Gundelfinger and H. Scheich, “A new method for reliable and efficient reconstruction of 3-dimensional images from autoradiographs of brain sections,” J Neurosci Methods, 84, pp. 77-86, 199815.R.C Gonzalez, Digital Image Processing, Reading, MA:Addison Wesley, 1992。