perlmeditation
lin0
<p>Greetings Fellow Monks,</p>
<p>For an [http://en.wikipedia.org/wiki/Image_analysis|image analysis] application I am writing using the [http://pdl.perl.org/|PDL], I needed to compute some [http://en.wikipedia.org/wiki/Texture|texture] measures. After some research, I decided to go with the measures proposed by [http://en.wikipedia.org/wiki/R_M_Haralick|Rober Haralick] based on the [http://en.wikipedia.org/wiki/Co-occurrence_matrix|Gray Level Co-occurrence Matrix (GLCM)]. To make a long story short, I found a [http://www.fp.ucalgary.ca/mhallbey/tutorial.htm|nice tutorial on the GLCM] and started implementing the code for computing the GLCM and the texture measures following the equations presented in the tutorial. This has shown me that I don't know PDL nearly as well as I would like. However, here is my first take to computing the GLCM and some of the texture measures:</p>
<readmore>
<code>
#!/usr/bin/perl
use warnings;
use strict;
use PDL;
use PDL::NiceSlice;
# ================================
# cooccurrence:
#
# $glcm = cooccurrence( $pdl, $dir, $dist, $symmetric )
#
# computes the grey level coocurrence coocurrence
# matrix of piddle $pdl for a given direction and
# distance
#
# Inputs:
# $pdl
# $dir: direction of evaluation
# $dir angle
# 0 +90
# 1 +45
# 2 0
# 3 -45
# 4 -90
# $dist: distance between pixels
# $symmetric: 0 => non-symmetric $glcm
#
# ================================
sub cooccurrence {
my ( $pdl, $dir, $dist, $symmetric ) = @_;
my $min_quantization_level = int( min( $pdl ) );
my $max_quantization_level = int( max( $pdl ) );
my $glcm = zeroes( $max_quantization_level
- $min_quantization_level + 1
, $max_quantization_level
- $min_quantization_level + 1 );
my ($dir_x, $dir_y);
if ( $dir == 0 ){
$dir_x = 0;
$dir_y = 1;
} elsif ( $dir == 1 ){
$dir_x = 1;
$dir_y = 1;
} elsif ( $dir == 2 ){
$dir_x = 1;
$dir_y = 0;
} elsif ( $dir == 3 ){
$dir_x = 1;
$dir_y = -1;
} elsif ( $dir == 4 ){
$dir_x = 0;
$dir_y = -1;
} else {
$dir_x = 0;
$dir_y = 0;
}
$dir_x *= $dist;
$dir_y *= $dist;
my $glcm_ind_x = 0;
my $glcm_ind_y = 0;
foreach my $grey_level_1 ( $min_quantization_level .. $max_quantization_level ){
my ( $ind_x_1, $ind_y_1 )
= whichND( $pdl == $grey_level_1 );
$ind_x_1 += $dir_x;
$ind_y_1 += $dir_y;
foreach my $grey_level_2 ( $min_quantization_level .. $max_quantization_level ){
my ( $ind_x_2, $ind_y_2 )
= whichND( $pdl == $grey_level_2 );
my $count = 0;
foreach my $i (0..$ind_x_1->getdim(0) - 1) {
foreach my $j (0..$ind_x_2->getdim(0) - 1) {
if ( ($ind_x_1($i) == $ind_x_2($j))
and ($ind_y_1($i) == $ind_y_2($j)) ) {
$count++;
}
}
}
$glcm( $glcm_ind_x, $glcm_ind_y ) .= $count;
$glcm_ind_y++;
}
$glcm_ind_y = 0;
$glcm_ind_x++;
}
if ( $symmetric ) {
$glcm += transpose( $glcm );
}
$glcm /= sum( $glcm );
return $glcm;
}
# ================================
# texture_descriptors:
#
# ( $contrast, $dissimilarity, $homogeneity
# , $inverse_difference, $asm, $energy )
# = texture_descriptors( $glcm );
#
# computes a set of texture descriptors
# associated with the GLCM $glcm
#
# $contrast:
# Range = [0 .. ($glcm->getdim(0)-1)^2]
# $contrast = 0 for a constant image.
# $homogeneity:
# Measures the closeness of the distribution
# of elements in the GLCM to the GLCM diagonal.
# Range = [0 1]
# $homogeneity is 1 for a diagonal GLCM.
# ================================
sub texture_descriptors{
my $glcm = pdl( @_ );
my $n = $glcm->getdim(0);
my $i = sequence( $n );
my $j = sequence( $n );
my $diff = $i->dummy(0, $n) - $j->dummy(1, $n);
my $contrast = sum( $glcm * ($diff ** 2) );
my $dissimilarity = sum( $glcm * abs( $diff ) );
my $homogeneity = sum( $glcm / ( 1 + $diff ** 2) );
my $inverse_difference = sum( $glcm / ( 1 + abs( $diff ) ) );
my $asm = sum( $glcm ** 2 );
my $energy = sqrt( $asm );
return ( $contrast, $dissimilarity, $homogeneity
, $inverse_difference, $asm, $energy );
}
my $pdl = pdl([0,0,1,1],[0,0,1,1],[0,2,2,2],[2,2,3,3]);
my $glcm = cooccurrence( $pdl, 2, 1, 1 );
print "glcm: $glcm\n";
my ( $contrast, $dissimilarity, $homogeneity
, $inverse_difference, $asm, $energy )
= texture_descriptors( $glcm );
print "contrast: $contrast\tdissimilarity: $dissimilarity\n";
print "homogeneity: $homogeneity\t";
print "inverse difference: $inverse_difference\n";
print "ASM: $asm\tenergy: $energy\n";
</code>
</readmore>
<p>I still have a lot more work to do, but feel free to go ahead and share with us your approach for finding these or the remaining texture measures.</p>
<p>Cheers,</p>
[lin0]