#!/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_quantiz
+ation_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_qua
+ntization_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";
[download]