Designing texture filters with genetic algorithms: An application to ...
SIGNW
PROCESSING Signal Processing
57 (1997) 19-33
Designing texture filters with genetic algorithms: An application to medical images Konstantinos Department
Delibasis, Peter E. Undrill*, George G. Cameron
of BiomedicalPhysics and Bioengineering, Vniversiv of Aberdeen, Foresterhill, Aberdeen, AB25 220, Scotland, UK Received 27 April 1995; revised 14 December
1995, 6 May 1996 and 28 October
1996
Abstract The problem of texture recognition is addressed by studying appropriate descriptors in the spatial frequency domain. During a training phase a filter is configured to determine different classes of texture by the response of its correlation with the Fourier spectrum of training-image templates. This is achieved by genetic algorithm-based optimisation. The technique is tested on standard texture patterns and then applied to magnetic resonance images of the brain to segment the cerebellum from the surrounding white and grey matter. Comparisons with established texture recognition techniques are presented, which show that the proposed method performs as well as, or better than, traditional techniques for the chosen instances of standard and anatomical texture and has the advantage of not having to decide which texture measure to use for a specific image structure. 0 1997 Elsevier Science B.V. Zusammenfassung In diesem Artikel wird das Problem der Texturerkennung dadurch gel&t, dal3 geeignete Abbildungen im RaumFrequenz-Bereich betrachtet werden. Verschiedene Texturklassen sollen mit Hilfe eines Filters erkannt werden, welches in einer Trainingsphase konfiguriert wird. Dabei wird die Korreliertheit der zu bestimmenden Klassen mit dem Fourierspektrum einiger Trainingsschablonen ausgenutzt. Die Optimierung des Filters basiert auf der Anwendung von Genetischen Algorithmen. Diese Technik wird zunlchst an Standardtexturmustern getestet und dann auf Magnetresonanzbilder des Gehirns angewandt, urn das Zerebellum von der umgebenden weigen und grauen Materie zu segmentieren. Vergleiche mit anerkannten Verfahren zur Texturerkennung zeigen, dab die vorgestellte Methode eine gleiche oder sogar bessere Leistungsfahigkeit als die traditionellen Techniken aufweist. Desweiteren besitzt sie den Vorteil der Unabhangigkeit der Megmethode von der speziellen Struktur des zu untersuchenden Bildes. 0 1997 Elsevier Science B.V. R&sum& Le probltme de reconnaissance de texture est abordt en Ctudiant des descripteurs appropribs dans le domaine frequentiel spatial. Un filtre est configure durant une phase d’apprentissage, qui permet de determiner les differentes classes de textures a l’aide de la reponse de leur correlation avec le spectre de Fourier de modeles d’images d’apprentissage. Ceci est effectut par optimisation baste sur un algorithme genetique. La technique est testee sur des motifs de textures standards, puis appliquee a des images par resonnance magnetique du cerveau, afin de segmenter le cerebellum de la mat&e grise et blanche qui l’entoure. Des comparaisons sont faites avec des techniques Ctablies de
*Corresponding
author.
0165-1684/97/$17.00 0 PIISO165-1684(96)00183-l
Tel.: 224 681 818; fax: 224 685 645; e-mail: [email protected].
1997 Elsevier Science B.V. All rights reserved.
20
K. Delibasis et al. / Signal Processing 57 (1997) 19-33
reconnaissance de structures, qui montrent que la methode proposee se comporte aussi bien, voire mieux, que le: techniques traditionelles pour les exemples choisis de textures standard et anatomiques, et qu’elle a l’avantage de ne pas ntcessiter de decision quant a la mesure de texture g utiliser pour une structure d’image spkcifique. 0 1997 Elseviei Science B.V. Keywords: Texture recognition; Image segmentation; Frequency descriptors; Genetic algorithms
1. Introduction
Medical images contain objects with characteristic shape and texture. Whilst several detection algorithms are edge based, detection, classification and determination of regional outlines can be made more robust if specific texture descriptors are established and then matched in an appropriate feature space. The objective of this work is the development of an algorithm which, after training, will be able to discriminate between different textures classes. The system proposed uses the properties of the Fourier spectrum of each pattern to achieve this discrimination,
1. I. Statistical approaches segmentation
to texture-based
Several techniques have been employed for texture-based segmentation, Sonka et al. [42] provides a comprehensive review of them. Most of them derive categories of texture descriptors and then, during a training phase, cluster these descriptors to achieve discrimination. Traditional methods of texture feature extraction are based either on statistical or structural models [13, 18, 381. In the statistical model, texture is defined by a characteristic set of relationships between image elements, and for most practical purposes these are determined from tonal values. Although their investigation is not the primary concern of this work, we will be using one or more of these as comparators and, therefore, they will be briefly mentioned. Grey-level spatial dependency (GLSD) matrices, or co-occurrence matrices, were first introduced by Haralick et al. [19]. Out of 20 descriptors proposed, three of the most powerful are contrast, entropy and energy. In segmentation applications, regions are labelled following an examination of
either first-order [44,45] or second-order [7] stat. istics between the tonal values of pixel pairs with relative orientation d, and displacement Zero Valued Gene
Used Harmonics
i
(a)
\ Complete 2D chromosome
Parents
(W Fig. 4. (a) Coding the Chromosome, segments of the 2D chromosome.
(b) a 2 x 1 crossover
of
25
Offspring 1 receives each segment from parent 1 with probability pcross (and the same segment from parent 2 with probability (1 - ppcross)).The other unused segment goes to the second offspring. Each bit of each segment is copied with probability (1 - p,,,), otherwise a mutation occurs which acts as a logical NOT on the value of the bit. Dependent on the choice of the cut-off frequency m, there will be sections of unused genetic material. Although in principle this might be seen as analogous to the existence of introns in the genetic material [30], and some theoretical specialist work [29] suggests that they can be beneficially used in the evolutionary process, we prohibit these sections of the spectral representation for being affected by simulated evolution. For most genetic processes, a binary alphabet seems intuitive, however non-binary alphabets have been studied where the representation is a continuous variable [23]. The actions of crossover and mutation are no longer unambiguous, the former could be either a geometric or arithmetic mean, or a magnitude adjustment based on a function of absolute difference. Mutation could be a complete replacement with a random value or the arithmetic modification by a random value. Because of these uncertainties the binary representation was retained. 2.3.3. Choice of GA parameters The GA population was limited to 100 and 100 generations of evolution allowed. Crossover probability was typically 0.8 and 2 breakpoints allowed in each dimension. Where 2 or more gene segments were selected, the probability of segment exchange was 0.25. Bit mutation was 0.5%. The quantities correlated were the logarithm of the real part of the spectral magnitude, encoded to 8 bits. We have experimented with changing the dynamic range of each gene, dependent on its distance from the DC point but have observed no improvement in performance, therefore a uniform range mapping was used throughout. Experimental variation of the genetic parameters usually failed to alter the final convergence, although changing the crossover strategy from l-breakpoint to n-breakpoints resulted in a halving of the number of generations needed to achieve a given convergence [9].
26
2.4. Implementation
K. Delibasis et al. J Signal Processing
of the system
(a) The starting point is an image with regions classified into two different classes. One or more training patches (either 16 x 16 or 32 x 32 pixels) are selected as members of each class. (b) The first generation of the GAS is initialised with random values. These are applied up to the cut-off frequency and the rest of the gene matrix set to zero. The design of a discriminating filter now proceeds by GA optimisation using the objective function and the training data. (c) This filter can be used either to find regions of these classes in the same image, or in a different image, containing similar target regions.
2.5. Post-processing segmentation
the results of texture-based
The result of the application of the GA-designed filter to an image is likely to be an incompletely segmented image, where regions having the texture which belongs to class A have high values and regions similar to members of class B have low values. 2.5.1. Extending the classification To improve the segmentation and to apply the method to more than two textures, a maximum likelihood decision rule that minimises the probability of false classification can be used (i) The members of each of n classes are selected as rectangular training samples. (ii) For i = 1 to n _ set class i to class A and the rest of the classes to class B; _ use the GA-based filter configuration system to design a filter fi which will produce high response for texture of class A and low response for class B; - filter the original image usingfi to produce image Ii. (iii) For all pixels in the segmented image _ locate the filtered image II, resulting from applying filter fk on the original image, in which the current pixel has value greater
57 (1997) 19-33
than the value of the same pixel in any of the remaining (n - 1) filtered images; _ set the current pixel in the segmented image to value k The image resulting from the above segmenta tion algorithm will contain values from 1 to n indicating the class of texture for each pixel. 2.5.2. Improving the classijcation The output of the segmentation algorithm can bc improved by applying a median filter to the seg mented image. The size of the filter found to bc useful was in the range from 5 x 5 to 11 x 11 pixels The median filter ensures that small islands o pixels are not classified as belonging to a textura region different from their surroundings.
2.6. Derivation of an enclosing contour For the MR image, there is the anatomical con straint of seeking a single feature, i.e. the enclosing contour. A ‘blob-colouring’ technique [2] labels al separate regions, so that those which are ana tomically small can be rejected, any remaining arc outlined. The candidate boundary that has been obtainer so far may not accurately enclose the region o contrasting texture because of the filter’s inheren spatial resolution, but it is likely that it will follow it; boundary and have a similar shape. It can be furthe refined if there is strong-edge information to bc exploited. The brain image is an example where the distinctly different texture of the cerebellum may bc enclosed in part by a strong boundary, informatiot that is not utilised by a texture segmentation tech nique. We therefore assist the performance of thl GA by an active contour model, or ‘snake’ [24,27] This snake will be initialised with the boundar: resulting from the texture segmentation and wil migrate, via energy minimisation, towards a bound ary suggested by edge-strength evidence in the image
2.7. Practical realisation The above algorithms were implemented in Pas cal on a SUN SparcStation 10-41. The design o
K. Delibasis et al. / Signal Processing 57 (1997) 19-33
a 16 x 16 pixel filter, using classes with three members each, takes 3-5 min and the application of the filter to a 256 x 256 pixel image 7 min. The overall approach therefore combines an initial segmentation focus based exclusively on texture characteristics with further refinements utilising regional and edge information. The use of a snake as a final refinement was employed only for the cerebellum segmentation.
3. Results 3. I. Standard texture patterns Our initial objective is to compare the operation of the GA-based texture classification with conventional methods on controlled texture. Fig. 3(a) contains three sections of real, but standardised, Brodatz texture. Textures 1, 2 and 3 represent the upper left, upper right and lower right quadrant of each 256 x 256 pixel image. The algorithm used to segment them was based on that described earlier: - four 32 x 32 pixel patches were extracted from
Type
Original Texture
GA method
27
each texture type, to serve as class members; _ the GA-based filter design algorithm was executed three times, producing 3 filters: 6) the first identifies texture 1, whereas it produces a low response for textures 2 and 3; (ii) the second filter identifies texture 2 and suppresses textures 1 and 3; (iii) the third filter identifies texture 3 and suppresses textures 2 and 1. _ the original image was filtered with all the above filters; _ the segmented image was produced by comparing the three filtered images on a pixel by pixel basis and setting the current pixel of the segmented image to a value equal to the rank of the image whose pixel had the maximum value compared to the other two images. Fig. 5 shows the result of the GA-based texture filter, using 32 x 32 pixel training segments, compared with the co-histogram method of Unser and the Laws R5R5 matrix. Each method is applied in turn to 3 montages of texture type, B1, B2 and B3 and smoothed using a 9 x 9 median filter prior to
Co-histogram method
Fig. 5. Comparison of classification of Brodatz texture types.
Laws method
28
K. Delibasis et al. J Signal Processing 57 (1997) 19-33
unsupervised classification, as described by Heckbert [20]. For Br, the GA method is visually superior, for Bz all methods are good, although co-occurrence histograms give more uniform regions and for B3 the Laws method is the best, followed by the GA method. Whilst it could be argued that any rigorous comparison necessitates the application of supervised classification with defined cluster centres, the objective here is to confirm the effectiveness of our GA method, which achieves classification through regional training rather than rigorous supervision, before moving to a practical medical imaging task.
3.2. Magnetic resonance brain images A more ambitious application was to design a filter able to discriminate between the characteristic texture of the cerebellum and surrounding tissue. Although conventional edge based approaches might be used (since an edge exists separating at least part of the cerebellum from the rest of the brain), the complexity of grey-scale variation within the cerebellum itself does not allow successful exploitation of such methods, and the shape is insufficiently characteristic of template matching. We used trans-axial slices from an MRI data set (236 x 176 pixels at 1 mm resolution). Fig. 6 shows the essential nature of the problem. The lobes containing white and grey matter are relatively uniform, connected regions with randomly folded outlines, whilst the brain stem is an enclosed uniform object. The texture of the cerebellum has a characteristic structure with longitudinal folds which give rise to a tonal variation having spatially dependent orientation, albeit in a generally circular manner. For the training, class A consists of four 16 x 16 pixel patches marked as black squares within the cerebellum and class B consists of four external patches, marked, in Fig. 7(a), as white squares. We designed a filter to produce high correlation response over class A texture (cerebellum) and low response over class B texture (non-cerebellum). The resulting image can be seen in Fig. 7(b,). Following this, classes A and B were swapped and the procedure repeated. The resultant image is shown in Fig. 7(b,). A pixel by pixel comparison of
Fig. 6. Outline of medical image segmentation task.
the two filtered images assigned each pixel to one of the two classes (cerebellum or non-cerebellum) and the segmented image is shown in Fig. 7(c). 3.2.1. A comparison
of texture segmentation
The effectiveness of our technique needs to be tested against other well-established texture classification methods. We use three descriptors derived from the GLSD matrices [19], namely the contrast, entropy and energy, together with the co-occurrence histograms and the fractal dimension. Fig. 8 shows the result of determining the segmentation boundaries from the same section of image as in Fig. 7(a) using (a) the fractal dimension and (b) the co-occurrence histograms method for which two training regions (32 x 32 pixels) were used for each texture class. The maximum likelihood decision rule was applied in each case. The fractal descriptor performed well, identifying most of the cerebellum tissue, but also produced areas of false identification. Co-occurrence histograms produced more extended regions of false identification. The GLSD descriptors (contrast, entropy and energy) are displayed as tonal maps in Figs. 8(c-e), with entropy the most characteristic of the three methods.
3.3. Boundary refinement Fig. 9(a) shows the image after 7 x 7 median filtering and boundary following. The primary contour lies on the cerebellum, albeit with a smaller, but similar, shape. This underestimation is caused by the filter “seeing” the image using a 16 x 16 pixel
K. Delibasis et al. / Signal Processing 57 (1997) 19-33
/
Locating texture A
/
Ia) . .
29
\ 4
Locating texture 6
Texture classification
Fig. 7. Segmentation of cerebellum using the GA-designed filter.
window and therefore an error in the boundary of a few pixels is expected. The large white rectangle indicates the extent of the image that has been examined by the filters, the much smaller contour is rejected by a magnitude criterion. The next step is to use the resulting contour to initialise the active snake, resulting in Fig. 9(b). The position of the cerebellum has been located and its shape accurately captured. Finally, we applied filters derived from one MRI trans-axial section to another taken at a different axial position where the texture may be slightly different due to biological change. The outcome is shown in Fig 9(c). Other regions as well as the cerebellum are marked on the brainstem and in the
surrounding brain tissue. The refined boundary is shown in Fig. 9(d), indicating that filters derived using one image can detect similar texture in another.
4. Conclusion A system for texture discrimination, based on the spectral frequency properties is described and results produced using images containing standard textures and magnetic resonance images of the human brain. The system exploits well-established Fourier spectral properties. The difference from other texture classification methods lies in the fact
30
K. Delibasis et al. / Signal Processing 57 (1991) 19-33
(a)
WI
le)
Fig. 8. Identification of cerebellum region by statistical methods.
that this approach searches the frequency spectrum of given textures, extracting characteristic relationships which are difficult to discover analytically. In an earlier investigation [9J we developed a GA-based method for designing texture filters for regular patterns with rotational invariance, angular discrimination and sensitivity to scale and intensity. This work extends the GA approach to real texture and a practical task of image classification and segmentation. In this respect the work has much in common, at least in terms of end-point objective, with the work of You and Cohen [47] in which the Laws matrix approach [ZS] is the starting point for the determination of a 5 x 5 mask with a more robust performance in classifying rotated and scaled textures. The GA approach with which we have experimented differs in two major features: (a) You and Cohen correctly identify the difficulties of gradient descent optimisation and propose a ‘guided search procedure’ using a vector of image dependent variables as an objective function. The guided search procedure optimises the magnitude
of this vector by proposing random changes on its components, taking the best resultant at each iteration. The GA approach also starts off with a randomised set of variables, and evolves towards a solution using a rigorously defined objective function, based on a stringent analysis of structural texture, designed to separate classes, through a process with similarities to natural selection. This allows the continuous adjustment of all components in the objective function (via the chromosome) rather than a single step-by-step approach. (b) The output from You’s work is a Laws-type discriminatory mask, 5 x 5 in size, with the sum of its terms equal to zero but not necessarily symmetric. The mask produced by the GA approach does not have this constraint, which may be desirable, but more importantly our mask is at least 16 x 16 pixels, having a greater discriminatory power for problems with considerable variation within the texture types, such as the medical imaging problem, albeit with additional computational complexity and the possible need
PROCESSING Signal Processing
57 (1997) 19-33
Designing texture filters with genetic algorithms: An application to medical images Konstantinos Department
Delibasis, Peter E. Undrill*, George G. Cameron
of BiomedicalPhysics and Bioengineering, Vniversiv of Aberdeen, Foresterhill, Aberdeen, AB25 220, Scotland, UK Received 27 April 1995; revised 14 December
1995, 6 May 1996 and 28 October
1996
Abstract The problem of texture recognition is addressed by studying appropriate descriptors in the spatial frequency domain. During a training phase a filter is configured to determine different classes of texture by the response of its correlation with the Fourier spectrum of training-image templates. This is achieved by genetic algorithm-based optimisation. The technique is tested on standard texture patterns and then applied to magnetic resonance images of the brain to segment the cerebellum from the surrounding white and grey matter. Comparisons with established texture recognition techniques are presented, which show that the proposed method performs as well as, or better than, traditional techniques for the chosen instances of standard and anatomical texture and has the advantage of not having to decide which texture measure to use for a specific image structure. 0 1997 Elsevier Science B.V. Zusammenfassung In diesem Artikel wird das Problem der Texturerkennung dadurch gel&t, dal3 geeignete Abbildungen im RaumFrequenz-Bereich betrachtet werden. Verschiedene Texturklassen sollen mit Hilfe eines Filters erkannt werden, welches in einer Trainingsphase konfiguriert wird. Dabei wird die Korreliertheit der zu bestimmenden Klassen mit dem Fourierspektrum einiger Trainingsschablonen ausgenutzt. Die Optimierung des Filters basiert auf der Anwendung von Genetischen Algorithmen. Diese Technik wird zunlchst an Standardtexturmustern getestet und dann auf Magnetresonanzbilder des Gehirns angewandt, urn das Zerebellum von der umgebenden weigen und grauen Materie zu segmentieren. Vergleiche mit anerkannten Verfahren zur Texturerkennung zeigen, dab die vorgestellte Methode eine gleiche oder sogar bessere Leistungsfahigkeit als die traditionellen Techniken aufweist. Desweiteren besitzt sie den Vorteil der Unabhangigkeit der Megmethode von der speziellen Struktur des zu untersuchenden Bildes. 0 1997 Elsevier Science B.V. R&sum& Le probltme de reconnaissance de texture est abordt en Ctudiant des descripteurs appropribs dans le domaine frequentiel spatial. Un filtre est configure durant une phase d’apprentissage, qui permet de determiner les differentes classes de textures a l’aide de la reponse de leur correlation avec le spectre de Fourier de modeles d’images d’apprentissage. Ceci est effectut par optimisation baste sur un algorithme genetique. La technique est testee sur des motifs de textures standards, puis appliquee a des images par resonnance magnetique du cerveau, afin de segmenter le cerebellum de la mat&e grise et blanche qui l’entoure. Des comparaisons sont faites avec des techniques Ctablies de
*Corresponding
author.
0165-1684/97/$17.00 0 PIISO165-1684(96)00183-l
Tel.: 224 681 818; fax: 224 685 645; e-mail: [email protected].
1997 Elsevier Science B.V. All rights reserved.
20
K. Delibasis et al. / Signal Processing 57 (1997) 19-33
reconnaissance de structures, qui montrent que la methode proposee se comporte aussi bien, voire mieux, que le: techniques traditionelles pour les exemples choisis de textures standard et anatomiques, et qu’elle a l’avantage de ne pas ntcessiter de decision quant a la mesure de texture g utiliser pour une structure d’image spkcifique. 0 1997 Elseviei Science B.V. Keywords: Texture recognition; Image segmentation; Frequency descriptors; Genetic algorithms
1. Introduction
Medical images contain objects with characteristic shape and texture. Whilst several detection algorithms are edge based, detection, classification and determination of regional outlines can be made more robust if specific texture descriptors are established and then matched in an appropriate feature space. The objective of this work is the development of an algorithm which, after training, will be able to discriminate between different textures classes. The system proposed uses the properties of the Fourier spectrum of each pattern to achieve this discrimination,
1. I. Statistical approaches segmentation
to texture-based
Several techniques have been employed for texture-based segmentation, Sonka et al. [42] provides a comprehensive review of them. Most of them derive categories of texture descriptors and then, during a training phase, cluster these descriptors to achieve discrimination. Traditional methods of texture feature extraction are based either on statistical or structural models [13, 18, 381. In the statistical model, texture is defined by a characteristic set of relationships between image elements, and for most practical purposes these are determined from tonal values. Although their investigation is not the primary concern of this work, we will be using one or more of these as comparators and, therefore, they will be briefly mentioned. Grey-level spatial dependency (GLSD) matrices, or co-occurrence matrices, were first introduced by Haralick et al. [19]. Out of 20 descriptors proposed, three of the most powerful are contrast, entropy and energy. In segmentation applications, regions are labelled following an examination of
either first-order [44,45] or second-order [7] stat. istics between the tonal values of pixel pairs with relative orientation d, and displacement Zero Valued Gene
Used Harmonics
i
(a)
\ Complete 2D chromosome
Parents
(W Fig. 4. (a) Coding the Chromosome, segments of the 2D chromosome.
(b) a 2 x 1 crossover
of
25
Offspring 1 receives each segment from parent 1 with probability pcross (and the same segment from parent 2 with probability (1 - ppcross)).The other unused segment goes to the second offspring. Each bit of each segment is copied with probability (1 - p,,,), otherwise a mutation occurs which acts as a logical NOT on the value of the bit. Dependent on the choice of the cut-off frequency m, there will be sections of unused genetic material. Although in principle this might be seen as analogous to the existence of introns in the genetic material [30], and some theoretical specialist work [29] suggests that they can be beneficially used in the evolutionary process, we prohibit these sections of the spectral representation for being affected by simulated evolution. For most genetic processes, a binary alphabet seems intuitive, however non-binary alphabets have been studied where the representation is a continuous variable [23]. The actions of crossover and mutation are no longer unambiguous, the former could be either a geometric or arithmetic mean, or a magnitude adjustment based on a function of absolute difference. Mutation could be a complete replacement with a random value or the arithmetic modification by a random value. Because of these uncertainties the binary representation was retained. 2.3.3. Choice of GA parameters The GA population was limited to 100 and 100 generations of evolution allowed. Crossover probability was typically 0.8 and 2 breakpoints allowed in each dimension. Where 2 or more gene segments were selected, the probability of segment exchange was 0.25. Bit mutation was 0.5%. The quantities correlated were the logarithm of the real part of the spectral magnitude, encoded to 8 bits. We have experimented with changing the dynamic range of each gene, dependent on its distance from the DC point but have observed no improvement in performance, therefore a uniform range mapping was used throughout. Experimental variation of the genetic parameters usually failed to alter the final convergence, although changing the crossover strategy from l-breakpoint to n-breakpoints resulted in a halving of the number of generations needed to achieve a given convergence [9].
26
2.4. Implementation
K. Delibasis et al. J Signal Processing
of the system
(a) The starting point is an image with regions classified into two different classes. One or more training patches (either 16 x 16 or 32 x 32 pixels) are selected as members of each class. (b) The first generation of the GAS is initialised with random values. These are applied up to the cut-off frequency and the rest of the gene matrix set to zero. The design of a discriminating filter now proceeds by GA optimisation using the objective function and the training data. (c) This filter can be used either to find regions of these classes in the same image, or in a different image, containing similar target regions.
2.5. Post-processing segmentation
the results of texture-based
The result of the application of the GA-designed filter to an image is likely to be an incompletely segmented image, where regions having the texture which belongs to class A have high values and regions similar to members of class B have low values. 2.5.1. Extending the classification To improve the segmentation and to apply the method to more than two textures, a maximum likelihood decision rule that minimises the probability of false classification can be used (i) The members of each of n classes are selected as rectangular training samples. (ii) For i = 1 to n _ set class i to class A and the rest of the classes to class B; _ use the GA-based filter configuration system to design a filter fi which will produce high response for texture of class A and low response for class B; - filter the original image usingfi to produce image Ii. (iii) For all pixels in the segmented image _ locate the filtered image II, resulting from applying filter fk on the original image, in which the current pixel has value greater
57 (1997) 19-33
than the value of the same pixel in any of the remaining (n - 1) filtered images; _ set the current pixel in the segmented image to value k The image resulting from the above segmenta tion algorithm will contain values from 1 to n indicating the class of texture for each pixel. 2.5.2. Improving the classijcation The output of the segmentation algorithm can bc improved by applying a median filter to the seg mented image. The size of the filter found to bc useful was in the range from 5 x 5 to 11 x 11 pixels The median filter ensures that small islands o pixels are not classified as belonging to a textura region different from their surroundings.
2.6. Derivation of an enclosing contour For the MR image, there is the anatomical con straint of seeking a single feature, i.e. the enclosing contour. A ‘blob-colouring’ technique [2] labels al separate regions, so that those which are ana tomically small can be rejected, any remaining arc outlined. The candidate boundary that has been obtainer so far may not accurately enclose the region o contrasting texture because of the filter’s inheren spatial resolution, but it is likely that it will follow it; boundary and have a similar shape. It can be furthe refined if there is strong-edge information to bc exploited. The brain image is an example where the distinctly different texture of the cerebellum may bc enclosed in part by a strong boundary, informatiot that is not utilised by a texture segmentation tech nique. We therefore assist the performance of thl GA by an active contour model, or ‘snake’ [24,27] This snake will be initialised with the boundar: resulting from the texture segmentation and wil migrate, via energy minimisation, towards a bound ary suggested by edge-strength evidence in the image
2.7. Practical realisation The above algorithms were implemented in Pas cal on a SUN SparcStation 10-41. The design o
K. Delibasis et al. / Signal Processing 57 (1997) 19-33
a 16 x 16 pixel filter, using classes with three members each, takes 3-5 min and the application of the filter to a 256 x 256 pixel image 7 min. The overall approach therefore combines an initial segmentation focus based exclusively on texture characteristics with further refinements utilising regional and edge information. The use of a snake as a final refinement was employed only for the cerebellum segmentation.
3. Results 3. I. Standard texture patterns Our initial objective is to compare the operation of the GA-based texture classification with conventional methods on controlled texture. Fig. 3(a) contains three sections of real, but standardised, Brodatz texture. Textures 1, 2 and 3 represent the upper left, upper right and lower right quadrant of each 256 x 256 pixel image. The algorithm used to segment them was based on that described earlier: - four 32 x 32 pixel patches were extracted from
Type
Original Texture
GA method
27
each texture type, to serve as class members; _ the GA-based filter design algorithm was executed three times, producing 3 filters: 6) the first identifies texture 1, whereas it produces a low response for textures 2 and 3; (ii) the second filter identifies texture 2 and suppresses textures 1 and 3; (iii) the third filter identifies texture 3 and suppresses textures 2 and 1. _ the original image was filtered with all the above filters; _ the segmented image was produced by comparing the three filtered images on a pixel by pixel basis and setting the current pixel of the segmented image to a value equal to the rank of the image whose pixel had the maximum value compared to the other two images. Fig. 5 shows the result of the GA-based texture filter, using 32 x 32 pixel training segments, compared with the co-histogram method of Unser and the Laws R5R5 matrix. Each method is applied in turn to 3 montages of texture type, B1, B2 and B3 and smoothed using a 9 x 9 median filter prior to
Co-histogram method
Fig. 5. Comparison of classification of Brodatz texture types.
Laws method
28
K. Delibasis et al. J Signal Processing 57 (1997) 19-33
unsupervised classification, as described by Heckbert [20]. For Br, the GA method is visually superior, for Bz all methods are good, although co-occurrence histograms give more uniform regions and for B3 the Laws method is the best, followed by the GA method. Whilst it could be argued that any rigorous comparison necessitates the application of supervised classification with defined cluster centres, the objective here is to confirm the effectiveness of our GA method, which achieves classification through regional training rather than rigorous supervision, before moving to a practical medical imaging task.
3.2. Magnetic resonance brain images A more ambitious application was to design a filter able to discriminate between the characteristic texture of the cerebellum and surrounding tissue. Although conventional edge based approaches might be used (since an edge exists separating at least part of the cerebellum from the rest of the brain), the complexity of grey-scale variation within the cerebellum itself does not allow successful exploitation of such methods, and the shape is insufficiently characteristic of template matching. We used trans-axial slices from an MRI data set (236 x 176 pixels at 1 mm resolution). Fig. 6 shows the essential nature of the problem. The lobes containing white and grey matter are relatively uniform, connected regions with randomly folded outlines, whilst the brain stem is an enclosed uniform object. The texture of the cerebellum has a characteristic structure with longitudinal folds which give rise to a tonal variation having spatially dependent orientation, albeit in a generally circular manner. For the training, class A consists of four 16 x 16 pixel patches marked as black squares within the cerebellum and class B consists of four external patches, marked, in Fig. 7(a), as white squares. We designed a filter to produce high correlation response over class A texture (cerebellum) and low response over class B texture (non-cerebellum). The resulting image can be seen in Fig. 7(b,). Following this, classes A and B were swapped and the procedure repeated. The resultant image is shown in Fig. 7(b,). A pixel by pixel comparison of
Fig. 6. Outline of medical image segmentation task.
the two filtered images assigned each pixel to one of the two classes (cerebellum or non-cerebellum) and the segmented image is shown in Fig. 7(c). 3.2.1. A comparison
of texture segmentation
The effectiveness of our technique needs to be tested against other well-established texture classification methods. We use three descriptors derived from the GLSD matrices [19], namely the contrast, entropy and energy, together with the co-occurrence histograms and the fractal dimension. Fig. 8 shows the result of determining the segmentation boundaries from the same section of image as in Fig. 7(a) using (a) the fractal dimension and (b) the co-occurrence histograms method for which two training regions (32 x 32 pixels) were used for each texture class. The maximum likelihood decision rule was applied in each case. The fractal descriptor performed well, identifying most of the cerebellum tissue, but also produced areas of false identification. Co-occurrence histograms produced more extended regions of false identification. The GLSD descriptors (contrast, entropy and energy) are displayed as tonal maps in Figs. 8(c-e), with entropy the most characteristic of the three methods.
3.3. Boundary refinement Fig. 9(a) shows the image after 7 x 7 median filtering and boundary following. The primary contour lies on the cerebellum, albeit with a smaller, but similar, shape. This underestimation is caused by the filter “seeing” the image using a 16 x 16 pixel
K. Delibasis et al. / Signal Processing 57 (1997) 19-33
/
Locating texture A
/
Ia) . .
29
\ 4
Locating texture 6
Texture classification
Fig. 7. Segmentation of cerebellum using the GA-designed filter.
window and therefore an error in the boundary of a few pixels is expected. The large white rectangle indicates the extent of the image that has been examined by the filters, the much smaller contour is rejected by a magnitude criterion. The next step is to use the resulting contour to initialise the active snake, resulting in Fig. 9(b). The position of the cerebellum has been located and its shape accurately captured. Finally, we applied filters derived from one MRI trans-axial section to another taken at a different axial position where the texture may be slightly different due to biological change. The outcome is shown in Fig 9(c). Other regions as well as the cerebellum are marked on the brainstem and in the
surrounding brain tissue. The refined boundary is shown in Fig. 9(d), indicating that filters derived using one image can detect similar texture in another.
4. Conclusion A system for texture discrimination, based on the spectral frequency properties is described and results produced using images containing standard textures and magnetic resonance images of the human brain. The system exploits well-established Fourier spectral properties. The difference from other texture classification methods lies in the fact
30
K. Delibasis et al. / Signal Processing 57 (1991) 19-33
(a)
WI
le)
Fig. 8. Identification of cerebellum region by statistical methods.
that this approach searches the frequency spectrum of given textures, extracting characteristic relationships which are difficult to discover analytically. In an earlier investigation [9J we developed a GA-based method for designing texture filters for regular patterns with rotational invariance, angular discrimination and sensitivity to scale and intensity. This work extends the GA approach to real texture and a practical task of image classification and segmentation. In this respect the work has much in common, at least in terms of end-point objective, with the work of You and Cohen [47] in which the Laws matrix approach [ZS] is the starting point for the determination of a 5 x 5 mask with a more robust performance in classifying rotated and scaled textures. The GA approach with which we have experimented differs in two major features: (a) You and Cohen correctly identify the difficulties of gradient descent optimisation and propose a ‘guided search procedure’ using a vector of image dependent variables as an objective function. The guided search procedure optimises the magnitude
of this vector by proposing random changes on its components, taking the best resultant at each iteration. The GA approach also starts off with a randomised set of variables, and evolves towards a solution using a rigorously defined objective function, based on a stringent analysis of structural texture, designed to separate classes, through a process with similarities to natural selection. This allows the continuous adjustment of all components in the objective function (via the chromosome) rather than a single step-by-step approach. (b) The output from You’s work is a Laws-type discriminatory mask, 5 x 5 in size, with the sum of its terms equal to zero but not necessarily symmetric. The mask produced by the GA approach does not have this constraint, which may be desirable, but more importantly our mask is at least 16 x 16 pixels, having a greater discriminatory power for problems with considerable variation within the texture types, such as the medical imaging problem, albeit with additional computational complexity and the possible need
