#!perl -w

=pod

=head1 NAME

Neuron::Morphometrics::Tools - Common functions for all morphometric modules

=head1 SYNOPSIS

Neuron::Morphometrics::Tools provides a number of functions used by more than
one of the Neuron::Morphometrics modules to avoid code duplication.

=head1 FUNCTIONS

=over 0

=item *
I<@list> = initSegments(I<neuron>);

initialises segments and puts a list in the exported variable @SEGLIST;
automatically called from Neuron::SWC::Read.

=item *
I<$n> = numCompartments(I<neuron>);

returns the number of compartments in the neuron.

=item *
I<$n> = (min|max|avg)BranchOrder(I<neuron>);

returns minimal, maximal or average branch order in the neuron.

=back

=head1 AUTHOR

Gunther Schmidl <gschmidl@gmx.at>

=cut

package Neuron::Morphometrics::Tools;
use strict;
use warnings;
no warnings 'recursion';

use vars qw($NEURON_VERSION);
$NEURON_VERSION = "0.1";

# prepare variable export

use vars qw(@ISA @EXPORT @EXPORT_OK $VERSION);
require Exporter;
@ISA = qw(Exporter);
@EXPORT = qw(numCompartments minmax minBranchOrder maxBranchOrder
    avgBranchOrder @SEGLIST @FORKLIST $PI $NEURON_VERSION);
@EXPORT_OK = qw(initSegments);
$VERSION = $NEURON_VERSION;






#==============================================================================
#
# initSegments(neuron) - initialise hash by depth-first traversal of tree and
# associating an (arbitrary) segment number and (relevant) branching depth with
# each compartment, as well as each compartment's area, volume, path length,
# path distance and euclidean distance from the soma, and Burke's taper rate.
# It also sets up @SEGLIST, a list of compartments belonging to each segment,
# and @FORKLIST, a list of bifurcations (and trifurcations, should such exist).
#
#==============================================================================

use vars qw($seg @SEGLIST @FORKLIST @arr $PI);

$PI = 3.1415926535897932384626433832795;

sub recurse
{
    my $n   = $_[0];                        # neuron node
    my $dep = $_[1];                        # current depth

    # put segment number and branch order in node
    $n->{'segnum'} = $seg;
    $n->{'order'}  = $dep;

    # area, volume, compartment length, path distance, euclidian distance
    if ($n->{'parent'} == -1)
    {
        $n->{'area'} = 0;
        $n->{'volume'} = 0;
        $n->{'length'} = 0;
        $n->{'distance'} = 0;
        $n->{'edistance'} = 0;
        $n->{'taper'} = 0;
    }
    else
    {
        my $parent = $arr[$n->{'parent'}];

        # length of current compartment
        my $height = sqrt((abs($parent->{'x'} - $n->{'x'}) ** 2) +
                          (abs($parent->{'y'} - $n->{'y'}) ** 2) +
                          (abs($parent->{'z'} - $n->{'z'}) ** 2));

        # euclidian distance of current compartment
        my $euclid = sqrt((abs($arr[1]->{'x'} - $n->{'x'}) ** 2) +
                          (abs($arr[1]->{'y'} - $n->{'y'}) ** 2) +
                          (abs($arr[1]->{'z'} - $n->{'z'}) ** 2));

        # Burke's taper rate
        my $taper = (($parent->{'length'} == 0) ? 0 :
           (($parent->{'d'} - $n->{'d'}) / ($height + $parent->{'length'})));

        # assign all the values
        $n->{'area'} = $n->{'d'} * $PI * $height;
        $n->{'volume'} = (($n->{'d'} / 2) ** 2) * $PI * $height;
        $n->{'length'} = $height;
        $n->{'distance'} = $parent->{'distance'} + $height;
        $n->{'edistance'} = $euclid;
        $n->{'taper'} = $taper;
    }

    # add this compartment to list of segments
    do
    {
        $SEGLIST[$seg] = [];

        # the zeroth element of the segment contains additional segment info
        my $temp = {};
        $temp->{'number'} = -1;
        $temp->{'segnum'} = $n->{'segnum'};
        $temp->{'parent'} = $n->{'parent'};
        $temp->{'order'} = $n->{'order'};
        $temp->{'distance'} = $n->{'distance'};
        $temp->{'edistance'} = $n->{'edistance'};
        $temp->{'taper'} = $n->{'taper'};
        $temp->{'initiald'} = $n->{'d'};

        push @{$SEGLIST[$seg]}, $temp;
    }
    unless exists($SEGLIST[$seg]);

    push @{$SEGLIST[$seg]}, $n->{'number'};

    if ($n->{'kids'})                       # not a bachelor node
    {
        my @kids = @{ $n->{'kids'} };
        my $numkids = @kids;

        # check if it's a bifurcation (stem doesn't count)
        if ($n->{'parent'} != -1 && $numkids > 1)
        {
            my $f = {};
            $f->{'number'} = $n->{'number'};  
            $f->{'order'} = $n->{'order'};
            $f->{'distance'} = $n->{'distance'};
            $f->{'edistance'} = $n->{'edistance'};            
            push @FORKLIST, $f;
        }

        # check for trifurcations (not at stem)
        if (($n->{'parent'} != -1) && ($numkids > 2))
        {
            print "warning: trifurcation at compartment $n->{'number'}\n";
        }

        # begin recursion
        if ($numkids > 1) { $dep++; }       # one order deeper at bifurcation

        foreach my $kid (@kids)
        {
            if ($numkids > 1) { $seg++; }
            recurse($kid, $dep);
        }

        if ($numkids > 1) { $dep--; }       # one order up at end of recursion
    }
}

sub initSegments(@)
{
    # initialise values for recursion
    @arr = @{$_[0]};
    $seg = 0;

    recurse($arr[1], -1);                   # recurse tree

    # set up additional per-segment statistics
    foreach my $s (@SEGLIST)
    {
        @{$s}[0]->{'endd'} = $arr[@{$s}[-1]]->{'d'};
    }

    return 0;
}






#==============================================================================
#
# numCompartments(neuron) - return number of compartments.
#
#==============================================================================

sub numCompartments(@)
{
    my $temp = @{$_[0]};                    # simply get length of array
    return $temp - 1;                       # minus one, as @neuron[0] is empty
}






#==============================================================================
#
# generic min/max subroutine
#
# usage: minmax(array, variable, value)
#
# variable    the hash key to look up
# value       0 for minimum, 1 for maximum, 2 for average
#
#==============================================================================

sub minmax
{
    my @arr = @{$_[0]};
    my $num = numCompartments(\@arr);

    my $comp;                               # workhorse value

    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 = 1; $i <= $num; $i++)
    {
        if ($_[2] == 0)                     # minimum
        {
            if ($arr[$i]->{"$_[1]"} < $comp)
            {
                $comp = $arr[$i]->{"$_[1]"};
            }
        }
        elsif ($_[2] == 1)                  # maximum
        {
            if ($arr[$i]->{"$_[1]"} > $comp)
            {
                $comp = $arr[$i]->{"$_[1]"};
            }
        }
        else                                # average
        {
            $comp = $comp + $arr[$i]->{"$_[1]"};
        }
    }

    if ($_[2] == 2)
    {
        $comp = $comp / $num;
    }

    return $comp;
}






#==============================================================================
#
# (min|max|avg)BranchOrder(neuron) - return min/max/average branch order.
#
#==============================================================================

sub minBranchOrder(@) { return 0; }

sub maxBranchOrder(@) { return minmax(\@{$_[0]}, "order", 1); }

sub avgBranchOrder(@) { return minmax(\@{$_[0]}, "order", 2); }