#!perl -w

=pod

=head1 NAME

Neuron::Morphometrics::Segment - Information on Segments

=head1 SYNOPSIS

Neuron::Morphometrics::Segment provides different functions to return
morphometric values concerning a neuron's segments, like number of segments,
average branch order etc.

=head1 FUNCTIONS

In all the following functions, the optional parameters I<min> and I<max> define
limits in the parameter specified by I<which>: branch order ("order"), path
distance ("distance") or euclidian distance ("edistance"); for example,
numSegments(\@myneuron, 2, 4, order) would only consider segments of branch
orders 2 through 4, inclusively. If parameters are given at all, both I<min>
and I<max> must be given.

=over 0

=item *
I<$n> = numSegments(I<neuron> [, I<min>, I<max>, I<which>]);

return number of segments.

=item *
I<$n> = contraction(I<neuron> [, I<min>, I<max>, I<which>]);

returns average contraction of segments.

=item *
I<$n> = (min|max|avg)Fragmentation(I<neuron> [, I<min>, I<max>, I<which>]);

returns minimal, maximal or average fragmentation.

=item *
I<$n> = (min|max|avg)SegmentEuclidianLength(I<neuron> [, I<min>, I<max>,
I<which>]);

returns minimal, maximal or average euclidian segment length.

=item *
I<$n> = (min|max|avg)SegmentPathLength(I<neuron> [, I<min>, I<max>,
I<which>]);

returns minimal, maximal or average segment path length.

=item *
I<$n> = (min|max|avg)SegmentTaper(I<neuron> [, I<min>, I<max>, I<which>]);

returns minimal, maximal or average segment taper.

=back

=head1 AUTHOR

Gunther Schmidl <gschmidl@gmx.at>

=cut

package Neuron::Morphometrics::Segment;
use strict;
use warnings;
no warnings 'recursion';

use Neuron::Morphometrics::Tools;

# prepare variable export

use vars qw(@ISA @EXPORT @EXPORT_OK $VERSION);
require Exporter;
@ISA = qw(Exporter);
@EXPORT = qw(numSegments minSegmentTaper avgSegmentTaper maxSegmentTaper
    minFragmentation avgFragmentation maxFragmentation minSegmentPathLength
    maxSegmentPathLength avgSegmentPathLength minSegmentEuclidianLength
    maxSegmentEuclidianLength avgSegmentEuclidianLength contraction
    minSegmentPathDistance maxSegmentPathDistance avgSegmentPathDistance
    minSegmentEuclidianDistance maxSegmentEuclidianDistance
    avgSegmentEuclidianDistance);
@EXPORT_OK = qw();
$VERSION = 0.01;






#==============================================================================
#
# numSegments(neuron [, min, max, which]) - return number of segments.
#
#==============================================================================

sub numSegments
{
    # get parameters
    my @arr = @{$_[0]};
    my $min = (defined($_[1]) ? $_[1] : 0);
    my $max = (defined($_[2]) ? $_[2] : maxBranchOrder(\@arr));
    my $wch = (defined($_[3]) ? $_[3] : "order");

    # only include stem if no limits given
    my $count = ((defined($_[1]) && defined($_[2])) ? 0 : 1);

    foreach my $seg (@SEGLIST)
    {
        # skip this segment if outside of boundaries
        next if $arr[@{$seg}[1]]->{$wch} < $min;
        next if $arr[@{$seg}[1]]->{$wch} > $max;

        $count++;
    }

    return $count;
}






#==============================================================================
#
# generic min/max subroutine for fragmentation
#
# usage: frminmax(array, value, lolimit, uplimit, which)
#
# value       0 for minimum, 1 for maximum, 2 for average
# lolimit     lower branch order limit
# hilimit     higher branch order limit
# which       order, distance, euclidian distance
#
#==============================================================================

sub frminmax
{
    my @arr = @{$_[0]};
    my $num = numSegments(\@arr);
    my $min = (defined($_[2]) ? $_[2] : 0);
    my $max = (defined($_[3]) ? $_[3] : maxBranchOrder(\@arr));
    my $wch = (defined($_[4]) ? $_[4] : "order");

    my $comp;                               # workhorse value
    my $count = 0;                          # for average

    SWITCH:
    {
        $comp = 65000, last SWITCH if ($_[1] == 0);
        $comp = -1,    last SWITCH if ($_[1] == 1);
        $comp = 0,     last SWITCH if ($_[1] == 2);
    }

    for (my $i = 0; $i < $num; $i++)
    {
        my $number = @{$SEGLIST[$i]} - 1;

        if (($arr[$SEGLIST[$i][1]]->{$wch} >= $min) &&
            ($arr[$SEGLIST[$i][1]]->{$wch} <= $max))
        {
            if ($_[1] == 0)                 # minimum
            {
                if ($number < $comp)
                {
                    $comp = $number;
                }
            }
            elsif ($_[1] == 1)              # maximum
            {
                if ($number > $comp)
                {
                    $comp = $number;
                }
            }
            else                            # average
            {
                $count++;
                $comp += $number;
            }
        }

        $SEGLIST[$i][0]->{'fragmentation'} = $number;
    }

    if ($_[1] == 2)
    {
        $comp = (($count == 0) ? 0 : ($comp / $count));
    }

    return $comp;
}






#==============================================================================
#
# (min|max|avg)Fragmentation(neuron [, min, max, which]) - return
# min/max/average fragmentation of segment.
#
#==============================================================================

sub minFragmentation
{
    return frminmax(\@{$_[0]}, 0, $_[1], $_[2], $_[3]);
}

sub maxFragmentation
{
    return frminmax(\@{$_[0]}, 1, $_[1], $_[2], $_[3]);
}

sub avgFragmentation
{
    return frminmax(\@{$_[0]}, 2, $_[1], $_[2], $_[3]);
}






#==============================================================================
#
# generic min/max subroutine for segment taper and euclidian/path length
#
# usage: stminmax(array, parameter, value, lolimit, uplimit, which)
#
# value       0 for minimum, 1 for maximum, 2 for average
# lolimit     lower branch order limit
# hilimit     higher branch order limit
# which       order, distance, euclidian distance
#
#==============================================================================

sub stminmax
{
    my @arr = @{$_[0]};
    my $num = numSegments(\@arr);
    my $min = (defined($_[3]) ? $_[3] : 0);
    my $max = (defined($_[4]) ? $_[4] : maxBranchOrder(\@arr));
    my $wch = (defined($_[5]) ? $_[5] : "order");

    my $comp;                               # workhorse value
    my $count = 0;                          # for average

    SWITCH:
    {
        $comp = 65000, last SWITCH if ($_[2] == 0);
        $comp = -1,    last SWITCH if ($_[2] == 1);
        $comp = 0,     last SWITCH if ($_[2] == 2);
    }

    for (my $i = 0; $i < $num; $i++)
    {
        my $first = $SEGLIST[$i][1];
        my $last = $SEGLIST[$i][-1];
        my $work;

        if (($arr[$first]->{$wch} >= $min) &&
            ($arr[$first]->{$wch} <= $max))
        {
            if ($_[1] eq "d")
            {
                $work = $arr[$last]->{"$_[1]"} / $arr[$first]->{"$_[1]"};
            }
            else
            {
                $work = abs($arr[$last]->{"$_[1]"} - $arr[$first]->{"$_[1]"});
            }

            if ($_[2] == 0)                 # minimum
            {
                if ($work < $comp)
                {
                    $comp = $work;
                }
            }
            elsif ($_[2] == 1)              # maximum
            {
                if ($work > $comp)
                {
                    $comp = $work;
                }
            }
            else                            # average
            {
                $count++;
                $comp += $work;
            }
        }

        if ($_[1] eq "d")            { $SEGLIST[$i][0]->{'taper'} = $work; }
        elsif ($_[1] eq "distance")  { $SEGLIST[$i][0]->{'length'} = $work; }
        elsif ($_[1] eq "edistance") { $SEGLIST[$i][0]->{'elength'} = $work; }
    }

    if ($_[2] == 2)
    {
        $comp = (($count == 0) ? 0 : ($comp / $count));
    }

    return $comp;
}






#==============================================================================
#
# (min|max|avg)SegmentTaper(neuron [, min, max, which]) - return 
# min/max/average segment taper.
#
#==============================================================================

sub minSegmentTaper
{
    return stminmax(\@{$_[0]}, "d", 0, $_[1], $_[2], $_[3]);
}

sub maxSegmentTaper
{
    return stminmax(\@{$_[0]}, "d", 1, $_[1], $_[2], $_[3]);
}

sub avgSegmentTaper
{
    return stminmax(\@{$_[0]}, "d", 2, $_[1], $_[2], $_[3]);
}






#==============================================================================
#
# (min|max|avg)SegmentPathLength(neuron [, min, max, which]) - return
# min/max/average segment path length.
#
#==============================================================================

sub minSegmentPathLength
{
    return stminmax(\@{$_[0]}, "distance", 0, $_[1], $_[2], $_[3]);
}

sub maxSegmentPathLength
{
    return stminmax(\@{$_[0]}, "distance", 1, $_[1], $_[2], $_[3]);
}

sub avgSegmentPathLength
{
    return stminmax(\@{$_[0]}, "distance", 2, $_[1], $_[2], $_[3]);
}






#==============================================================================
#
# (min|max|avg)SegmentEuclidianLength(neuron [, min, max, which]) - return
# min/max/average segment path length.
#
#==============================================================================

sub minSegmentEuclidianLength
{
    return stminmax(\@{$_[0]}, "edistance", 0, $_[1], $_[2], $_[3]);
}

sub maxSegmentEuclidianLength
{
    return stminmax(\@{$_[0]}, "edistance", 1, $_[1], $_[2], $_[3]);
}

sub avgSegmentEuclidianLength
{
    return stminmax(\@{$_[0]}, "edistance", 2, $_[1], $_[2], $_[3]);
}






#==============================================================================
#
# contraction(neuron [, min, max, which]) - return euclidian length / path 
# length.
#
#==============================================================================

sub contraction
{
    my @arr = @{$_[0]};
    my $num = numSegments(\@arr);
    my $min = (defined($_[1]) ? $_[1] : 0);
    my $max = (defined($_[2]) ? $_[2] : maxBranchOrder(\@arr));
    my $wch = (defined($_[3]) ? $_[3] : "order");

    my $comp;                               # workhorse value
    my $count = 0;                          # for average

    for (my $i = 0; $i < $num; $i++)
    {
        my $first = $SEGLIST[$i][1];
        my $last = $SEGLIST[$i][-1];

        my $e = 0;                          # euclidian distance
        my $p = 0;                          # path distance

        if (($arr[$first]->{$wch} >= $min) &&
            ($arr[$first]->{$wch} <= $max))
        {
            $e = abs($arr[$last]->{'edistance'} - $arr[$first]->{'edistance'});
            $p = abs($arr[$last]->{'distance'} - $arr[$first]->{'distance'});

            if ($p == 0) { $p = $arr[$first]->{'length'}; }
            if ($e == 0) { $e = $arr[$first]->{'length'}; }

            my $temp = ($e / $p);

            $SEGLIST[$i][0]->{'contraction'} = $temp;
            $comp += $temp;
            $count++;
        }
    }

    return (($count == 0) ? 0 : ($comp / $count));
}







#==============================================================================
#
# generic min/max subroutine for euclidian/path distance
#
# usage: diminmax(array, parameter, value, lolimit, uplimit, which)
#
# value       0 for minimum, 1 for maximum, 2 for average
# lolimit     lower branch order limit
# hilimit     higher branch order limit
# which       order, distance, euclidian distance
#
#==============================================================================

sub diminmax
{
    my @arr = @{$_[0]};
    my $num = numSegments(\@arr);
    my $min = (defined($_[3]) ? $_[3] : 0);
    my $max = (defined($_[4]) ? $_[4] : maxBranchOrder(\@arr));
    my $wch = (defined($_[5]) ? $_[5] : "order");

    my $comp;                               # workhorse value
    my $count = 0;                          # for average

    SWITCH:
    {
        $comp = 65000, last SWITCH if ($_[2] == 0);
        $comp = -1,    last SWITCH if ($_[2] == 1);
        $comp = 0,     last SWITCH if ($_[2] == 2);
    }

    for (my $i = 0; $i < $num; $i++)
    {
        my $work = 0;

        if (($arr[$SEGLIST[$i][1]]->{$wch} >= $min) &&
            ($arr[$SEGLIST[$i][1]]->{$wch} <= $max))
        {
            $work = $arr[$SEGLIST[$i][1]]->{"$_[1]"};

            if ($_[2] == 0)                 # minimum
            {
                if ($work < $comp)
                {
                    $comp = $work;
                }
            }
            elsif ($_[2] == 1)              # maximum
            {
                if ($work > $comp)
                {
                    $comp = $work;
                }
            }
            else                            # average
            {
                $count++;
                $comp += $work;
            }
        }
    }

    if ($_[2] == 2)
    {
        $comp = (($count == 0) ? 0 : ($comp / $count));
    }

    return $comp;
}






#==============================================================================
#
# (min|max|avg)SegmentPathDistance(neuron [, min, max, which]) - return
# min/max/average segment path distance.
#
#==============================================================================

sub minSegmentPathDistance
{
    return diminmax(\@{$_[0]}, "distance", 0, $_[1], $_[2], $_[3]);
}

sub maxSegmentPathDistance
{
    return diminmax(\@{$_[0]}, "distance", 1, $_[1], $_[2], $_[3]);
}

sub avgSegmentPathDistance
{
    return diminmax(\@{$_[0]}, "distance", 2, $_[1], $_[2], $_[3]);
}






#==============================================================================
#
# (min|max|avg)SegmentEuclidianDistance(neuron [, min, max, which]) - return
# min/max/average euclidian segment distance.
#
#==============================================================================

sub minSegmentEuclidianDistance
{
    return diminmax(\@{$_[0]}, "edistance", 0, $_[1], $_[2], $_[3]);
}

sub maxSegmentEuclidianDistance
{
    return diminmax(\@{$_[0]}, "edistance", 1, $_[1], $_[2], $_[3]);
}

sub avgSegmentEuclidianDistance
{
    return diminmax(\@{$_[0]}, "edistance", 2, $_[1], $_[2], $_[3]);
}