#!perl -w

=pod

=head1 NAME

Neuron::Morphometrics::Bifurcation - Information on Bifurcations

=head1 SYNOPSIS

Neuron::Morphometrics::Bifurcation provides different functions to return
morphometric values concerning a neuron's birfurcations, like number of
bifurcations, asymmetry 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> = avgAsymmetry(I<neuron> [, I<min>, I<max>, I<which>]);

return average asymmetry of neuron.

=item *
I<$n> = avgDaughterRatio(I<neuron> [, I<min>, I<max>, I<which>]);

return average daughter ratio of neuron.

item*
I<$n> = bifurcationTorqueAngle(I<neuron> [, I<min>, I<max>, I<which>]);

return average bifurcation torque angle of bifurcations.

=item*
I<$n> = compartmentAmplitudeAngle(I<neuron> [, I<min>, I<max>, I<which>]);

return average amplitude angle of bifurcations with respect to their daughter
compartments.

=item*
I<$n> = diameterRatio(I<neuron> [, I<min>, I<max>, I<which>]);

return average parent/child diameter ratio.

=item *
I<$n> = numBifurcations(I<neuron> [, I<min>, I<max>, I<which>]);

return number of bifurcations.

=item*
I<$n> = powerCorrection(I<neuron>, I<rall> [, I<min>, I<max>, I<which>]);

return average power correction for a given Rall's Value.

=item*
I<$n> = rall(I<neuron> [, I<min>, I<max>, I<which>]);

return Rall's Value.

=item*
I<$n> = segmentAmplitudeAngle(I<neuron> [, I<min>, I<max>, I<which>]);

return average amplitude angle of bifurcations with respect to their daughter
segments.

=item*
I<$n> = tilt(I<neuron> [, I<min>, I<max>, I<which>]);

return average tilt between bifurcations.

=back

=head1 AUTHOR

Gunther Schmidl <gschmidl@gmx.at>

=cut

package Neuron::Morphometrics::Bifurcation;
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(numBifurcations avgDaughterRatio rall avgAsymmetry tilt
    compartmentAmplitudeAngle segmentAmplitudeAngle bifurcationTorqueAngle
    diameterRatio powerCorrection);
@EXPORT_OK = qw();
$VERSION = 0.01;






#==============================================================================
#
# numBifurcations(neuron [, min, max, which]) - return number of bifurcations.
#
#==============================================================================

sub numBifurcations
{
    # get parameters
    my @arr = @{$_[0]};
    my $min = (defined($_[1]) ? $_[1] : 0);
    my $max = (defined($_[2]) ? $_[2] : maxBranchOrder(\@arr));
    my $wch = (defined($_[3]) ? $_[3] : "order");

    my $num = numCompartments(\@arr);
    my $count = 0;
    foreach my $f (@FORKLIST)
    {
        my $i = $f->{'number'};
        if (($arr[$i]->{$wch} >= $min) && ($arr[$i]->{$wch} <= $max))
        {
            $count++;
        }
    }
    return $count;
}






#==============================================================================
#
# avgDaughterRatio(neuron [, min, max, which]) - average daughter ratio,
# defined as (diameter of left child) / (diameter of right child).
#
#==============================================================================

sub avgDaughterRatio
{
    # get parameters
    my @arr = @{$_[0]};
    my $min = (defined($_[1]) ? $_[1] : 0);
    my $max = (defined($_[2]) ? $_[2] : maxBranchOrder(\@arr));
    my $wch = (defined($_[3]) ? $_[3] : "order");

    my $num = numCompartments(\@arr);
    my $count = 0;                          # number of relevant bifurcations
    my $sum = 0;                            # sum of daughter ratios

    foreach my $f (@FORKLIST)
    {
        my $i = $f->{'number'};
        if (($arr[$i]->{$wch} >= $min) && ($arr[$i]->{$wch} <= $max))
        {
            # get the diameters of both children, divide them
            my @temp = @{$arr[$i]->{'kids'}};

            # assumption: smaller / larger
            if ($temp[0]->{'d'} > $temp[1]->{'d'})
            {
                ($temp[0], $temp[1]) = ($temp[1], $temp[0]);
            }
            my $tmp = ($temp[0]->{'d'} / $temp[1]->{'d'});
            $sum = $sum + $tmp;
            $f->{'kidratio'} = $tmp;

            $count++;
        }
    }
    return ($count == 0 ? 0 : $sum / $count);
}






#==============================================================================
#
# avgAsymmetry(neuron [, min, max, which]) - average asymmetry value, defined
# as |l-d|/(l+d-2) where l, d are the number of terminals on the left/right
# side of a bifurcation.
#
#==============================================================================

sub avgAsymmetry
{
    # get parameters
    my @arr = @{$_[0]};
    my $min = (defined($_[1]) ? $_[1] : 0);
    my $max = (defined($_[2]) ? $_[2] : maxBranchOrder(\@arr));
    my $wch = (defined($_[3]) ? $_[3] : "order");

    my $asym = 0;                           # total asymmetry
    my $count = 0;
    my @work;                               # workhorse array
    my %done;                               # avoid duplicate work

    # create list of all terminals
    foreach my $seg (@SEGLIST)
    {
        if (not($arr[@{$seg}[-1]]->{'kids'}))
        {
            push @work, $seg;
        }
    }

    while (@work > 0)
    {
        my $current = shift @work;

        next if exists $done{@{$current}[1]};

        if (not($arr[@{$current}[-1]]->{'kids'}))
        {
            # it's a terminal, so it counts 1 towards total asymmetry, which is
            # arbitrarily assigned to the terminal's lambda
            $arr[@{$current}[1]]->{'lambda'} = 1;
            $arr[@{$current}[1]]->{'delta'} = 0;
            $done{@{$current}[1]} = 1;
        }
        else
        {
            # it's a bifurcation, which means we add all the terminals on the
            # left side to get the new lambda and all the terminals on the
            # right side to get the new delta

            # if one of the values is not yet defined, push the segment to
            # the back of the list for a later try
            if (!defined(@{$arr[@{$current}[-1]]->{'kids'}}[0]->{'lambda'}) ||
                !defined(@{$arr[@{$current}[-1]]->{'kids'}}[1]->{'lambda'}) ||
                !defined(@{$arr[@{$current}[-1]]->{'kids'}}[0]->{'delta'}) ||
                !defined(@{$arr[@{$current}[-1]]->{'kids'}}[1]->{'delta'}))
            {
                push @work, $current;
            }
            else
            {
                my $l = @{$arr[@{$current}[-1]]->{'kids'}}[0]->{'lambda'};
                $l += @{$arr[@{$current}[-1]]->{'kids'}}[0]->{'delta'};

                my $r = @{$arr[@{$current}[-1]]->{'kids'}}[1]->{'lambda'};
                $r += @{$arr[@{$current}[-1]]->{'kids'}}[1]->{'delta'};

                $arr[@{$current}[1]]->{'lambda'} = $l;
                $arr[@{$current}[1]]->{'delta'} = $r;

                # count only if it's inside the correct limits
                if (($arr[@{$current}[1]]->{$wch} >= $min) &&
                    ($arr[@{$current}[1]]->{$wch} <= $max))
                {
                    $asym += (abs($l - $r) / ($l + $r - 2))
                        unless (($l + $r - 2) == 0);

                    $count++;
                }

                $done{@{$current}[0]} = 1;
            }
        }

        # finally, add the parent segment to the list, unless it would be the
        # stem
        foreach my $seg (@SEGLIST)
        {
            if ($arr[@{$current}[1]]->{'parent'} != 1 &&
            $arr[@{$seg}[-1]]->{'number'} == $arr[@{$current}[1]]->{'parent'})
            {
                push @work, $seg;
            }
        }
    }
    return (($count == 0) ? 0 : ($asym / $count));
}






#==============================================================================
#
# acos(i) - return inverse cosine of number.
#
#==============================================================================

sub acos
{
    return atan2(sqrt(1 - ($_[0] * $_[0])), $_[0]);
}






#==============================================================================
#
# compartmentAmplitudeAngle(neuron [, min, max, which]) - average compartment
# amplitude angle of the neuron. The angle is defined as
#
#             A x B
#  cos phi = -------, where A and B are the respective vectors.
#            |A|*|B|
#
#==============================================================================

sub compartmentAmplitudeAngle
{
    # get parameters
    my @arr = @{$_[0]};
    my $min = (defined($_[1]) ? $_[1] : 0);
    my $max = (defined($_[2]) ? $_[2] : maxBranchOrder(\@arr));
    my $wch = (defined($_[3]) ? $_[3] : "order");

    my $phi = 0;
    my $count = 0;

    foreach my $f (@FORKLIST)
    {
        my $fork = $f->{'number'};
        if (($arr[$fork]->{$wch} >= $min) &&
            ($arr[$fork]->{$wch} <= $max))
        {
            # get vectors
            my $ax = $arr[$fork]->{'kids'}[0]->{'x'} - $arr[$fork]->{'x'};
            my $ay = $arr[$fork]->{'kids'}[0]->{'y'} - $arr[$fork]->{'y'};
            my $az = $arr[$fork]->{'kids'}[0]->{'z'} - $arr[$fork]->{'z'};

            my $bx = $arr[$fork]->{'kids'}[1]->{'x'} - $arr[$fork]->{'x'};
            my $by = $arr[$fork]->{'kids'}[1]->{'y'} - $arr[$fork]->{'y'};
            my $bz = $arr[$fork]->{'kids'}[1]->{'z'} - $arr[$fork]->{'z'};

            my $dividend = $ax * $bx + $ay * $by + $az * $bz;
            my $divisor = sqrt($ax ** 2 + $ay ** 2 + $az ** 2) *
                sqrt($bx ** 2 + $by ** 2 + $bz ** 2);
            if ($divisor != 0)
            {
                my $temp = (acos($dividend / $divisor) / $PI) * 180;
                $phi += $temp;
                $f->{'caangle'} = $temp;
                $count++;
            }
        }
    }

    return (($count == 0) ? 0 : ($phi / $count));
}






#==============================================================================
#
# segmentAmplitudeAngle(neuron [, min, max, which]) - average segment
# amplitude angle of the neuron. Definition as above.
#
#==============================================================================

sub segmentAmplitudeAngle
{
    # get parameters
    my @arr = @{$_[0]};
    my $min = (defined($_[1]) ? $_[1] : 0);
    my $max = (defined($_[2]) ? $_[2] : maxBranchOrder(\@arr));
    my $wch = (defined($_[3]) ? $_[3] : "order");

    my $phi = 0;
    my $count = 0;

    foreach my $f (@FORKLIST)
    {
        my $fork = $f->{'number'};
        if (($arr[$fork]->{$wch} >= $min) &&
            ($arr[$fork]->{$wch} <= $max))
        {
            # get vectors
            my $lseg = $arr[$fork]->{'kids'}[0]->{'segnum'};
            my $rseg = $arr[$fork]->{'kids'}[1]->{'segnum'};

            my $left = $SEGLIST[$lseg][-1];
            my $right = $SEGLIST[$rseg][-1];

            my $ax = $arr[$left]->{'x'} - $arr[$fork]->{'x'};
            my $ay = $arr[$left]->{'y'} - $arr[$fork]->{'y'};
            my $az = $arr[$left]->{'z'} - $arr[$fork]->{'z'};

            my $bx = $arr[$right]->{'x'} - $arr[$fork]->{'x'};
            my $by = $arr[$right]->{'y'} - $arr[$fork]->{'y'};
            my $bz = $arr[$right]->{'z'} - $arr[$fork]->{'z'};

            my $dividend = $ax * $bx + $ay * $by + $az * $bz;
            my $divisor = sqrt($ax ** 2 + $ay ** 2 + $az ** 2) *
                sqrt($bx ** 2 + $by ** 2 + $bz ** 2);

            if ($divisor != 0)
            {
            	my $temp = (acos($dividend / $divisor) / $PI) * 180;
                $phi += $temp;
                $f->{'saangle'} = $temp;
                $count++;
            }
        }
    }

    return (($count == 0) ? 0 : ($phi / $count));
}






#==============================================================================
#
# bifurcationTorqueAngle(neuron [, min, max, which]) - average bifurcation
# torque angle of the neuron. Defined as the angle between the vectors
# (bifurcation - parent) and (left child + right child - bifurcation).
#
#==============================================================================

sub bifurcationTorqueAngle
{
    # get parameters
    my @arr = @{$_[0]};
    my $min = (defined($_[1]) ? $_[1] : 0);
    my $max = (defined($_[2]) ? $_[2] : maxBranchOrder(\@arr));
    my $wch = (defined($_[3]) ? $_[3] : "order");

    my $phi = 0;
    my $count = 0;

    foreach my $f (@FORKLIST)
    {
        my $fork = $f->{'number'};
        if (($arr[$fork]->{$wch} >= $min) &&
            ($arr[$fork]->{$wch} <= $max))
        {
            # span a parallelogram with the two daughter vectors
            my $tx = $arr[$fork]->{'kids'}[0]->{'x'} - $arr[$fork]->{'x'};
            my $ty = $arr[$fork]->{'kids'}[0]->{'y'} - $arr[$fork]->{'y'};
            my $tz = $arr[$fork]->{'kids'}[0]->{'z'} - $arr[$fork]->{'z'};

            my $qx = $arr[$fork]->{'kids'}[1]->{'x'} - $arr[$fork]->{'x'};
            my $qy = $arr[$fork]->{'kids'}[1]->{'y'} - $arr[$fork]->{'y'};
            my $qz = $arr[$fork]->{'kids'}[1]->{'z'} - $arr[$fork]->{'z'};

            # get end point
            $tx += $qx; $ty += $qy; $tz += $qz;

            # get the vector between that new end point and the bifurcation
            my $ax = $tx - $arr[$fork]->{'x'};
            my $ay = $ty - $arr[$fork]->{'y'};
            my $az = $tz - $arr[$fork]->{'z'};

            # get the vector between the bifurcation and its parent
            my $bx = $arr[$fork]->{'x'} - $arr[$arr[$fork]->{'parent'}]->{'x'};
            my $by = $arr[$fork]->{'y'} - $arr[$arr[$fork]->{'parent'}]->{'y'};
            my $bz = $arr[$fork]->{'z'} - $arr[$arr[$fork]->{'parent'}]->{'z'};

            my $dividend = $ax * $bx + $ay * $by + $az * $bz;
            my $divisor = sqrt($ax ** 2 + $ay ** 2 + $az ** 2) *
                sqrt($bx ** 2 + $by ** 2 + $bz ** 2);
            
            if ($divisor != 0)
            {
                my $temp = (acos($dividend / $divisor) / $PI) * 180;
                $phi += $temp;
                $f->{'btangle'} = $temp;
                $count++;
            }
        }
    }

    return (($count == 0) ? 0 : ($phi / $count));
}






#==============================================================================
#
# tilt(neuron [, min, max, which]) - average tilt of the neuron.
#
#==============================================================================

sub tilt
{
    # get parameters
    my @arr = @{$_[0]};
    my $min = (defined($_[1]) ? $_[1] : 0);
    my $max = (defined($_[2]) ? $_[2] : maxBranchOrder(\@arr));
    my $wch = (defined($_[3]) ? $_[3] : "order");


    my $phi = 0;
    my $count = 0;

    foreach my $f (@FORKLIST)
    {
        my $fork = $f->{'number'};

        if (($arr[$fork]->{$wch} >= $min) &&
            ($arr[$fork]->{$wch} <= $max))
        {
            # get vectors of plane spanned by vectors of current fork
            my $lseg = $arr[$fork]->{'kids'}[0]->{'segnum'};
            my $rseg = $arr[$fork]->{'kids'}[1]->{'segnum'};

            my $left = $SEGLIST[$lseg][-1];
            my $right = $SEGLIST[$rseg][-1];

            my $cax = $arr[$left]->{'x'} - $arr[$fork]->{'x'};
            my $cay = $arr[$left]->{'y'} - $arr[$fork]->{'y'};
            my $caz = $arr[$left]->{'z'} - $arr[$fork]->{'z'};

            my $cbx = $arr[$right]->{'x'} - $arr[$fork]->{'x'};
            my $cby = $arr[$right]->{'y'} - $arr[$fork]->{'y'};
            my $cbz = $arr[$right]->{'z'} - $arr[$fork]->{'z'};

            # get vectors of plane spanned by vectors of previous fork
            my $t = $arr[$fork]->{'parent'};
            my $prefork = $arr[$SEGLIST[$arr[$t]->{'segnum'}][1]]->{'parent'};

            # this doesn't work if the previous fork is the stem, so skip
            $f->{'tilt'} = -1 if ($prefork == -1);
            next if ($prefork == -1);

            # otherwise, go on processing
            $lseg = $arr[$prefork]->{'kids'}[0]->{'segnum'};
            $rseg = $arr[$prefork]->{'kids'}[1]->{'segnum'};
            $left = $SEGLIST[$lseg][-1];
            $right = $SEGLIST[$rseg][-1];

            my $pax = $arr[$left]->{'x'} - $arr[$prefork]->{'x'};
            my $pay = $arr[$left]->{'y'} - $arr[$prefork]->{'y'};
            my $paz = $arr[$left]->{'z'} - $arr[$prefork]->{'z'};

            my $pbx = $arr[$right]->{'x'} - $arr[$prefork]->{'x'};
            my $pby = $arr[$right]->{'y'} - $arr[$prefork]->{'y'};
            my $pbz = $arr[$right]->{'z'} - $arr[$prefork]->{'z'};

            # get cross products of both pairs of vectors (normal to planes)
            my $cnx = ($cay * $cbz) - ($caz * $cby);
            my $cny = ($caz * $cbx) - ($cax * $cbz);
            my $cnz = ($cax * $cby) - ($cay * $cbx);

            my $pnx = ($pay * $pbz) - ($paz * $pby);
            my $pny = ($paz * $pbx) - ($pax * $pbz);
            my $pnz = ($pax * $pby) - ($pay * $pbx);

            # and now the arccos (A x B / |A| * |B|) formula as before
            my $dividend = $cnx * $pnx + $cny * $pny + $cnz * $pnz;
            my $divisor = sqrt($cnx ** 2 + $cny ** 2 + $cnz ** 2) *
                sqrt($pnx ** 2 + $pny ** 2 + $pnz ** 2);
            if ($divisor != 0)
            {
                my $temp = (acos($dividend / $divisor) / $PI) * 180;
                $phi += $temp;
                $f->{'tilt'} = $temp;
                $count++;
            }
        }
    }

    return (($count == 0) ? 0 : ($phi / $count));
}






#==============================================================================
#
# diameterRatio(neuron [, min, max, which]) - average ratio between parent and
# child diameters.
#
#==============================================================================

sub diameterRatio
{
    # get parameters
    my @arr = @{$_[0]};
    my $min = (defined($_[1]) ? $_[1] : 0);
    my $max = (defined($_[2]) ? $_[2] : maxBranchOrder(\@arr));
    my $wch = (defined($_[3]) ? $_[3] : "order");

    my $work = 0;
    my $count = 0;

    foreach my $f (@FORKLIST)
    {
        my $fork = $f->{'number'};
        if (($arr[$fork]->{$wch} >= $min) &&
            ($arr[$fork]->{$wch} <= $max))
        {
            # get diameters
            my $dp = $arr[$fork]->{'d'};
            my $da = $arr[$fork]->{'kids'}[0]->{'d'};
            my $db = $arr[$fork]->{'kids'}[1]->{'d'};

            $work += $dp / $da; $count++;
            $work += $dp / $db; $count++;
            $f->{'dratio'} = (($dp / $da) + ($dp / $db)) / 2;
        }
    }

    return (($count == 0) ? 0 : ($work / $count));
}






#==============================================================================
#
# powerCorrection(neuron, rall, [, min, max, which]) - power correction.
#
#==============================================================================

sub powerCorrection
{
    # get parameters
    my @arr = @{$_[0]};
    my $rll = (defined($_[1]) ? $_[1] : rall(\@arr));
    my $min = (defined($_[2]) ? $_[2] : 0);
    my $max = (defined($_[3]) ? $_[3] : maxBranchOrder(\@arr));
    my $wch = (defined($_[4]) ? $_[4] : "order");

    my $pc = 0;
    my $count = 0;

    foreach my $f (@FORKLIST)
    {
        my $fork = $f->{'number'};
        if (($arr[$fork]->{$wch} >= $min) &&
            ($arr[$fork]->{$wch} <= $max))
        {
            # get diameters
            my $dp = $arr[$fork]->{'d'};
            my $da = $arr[$fork]->{'kids'}[0]->{'d'};
            my $db = $arr[$fork]->{'kids'}[1]->{'d'};

            my $temp = ($dp ** $rll) / ($da ** $rll + $db ** $rll);
            $pc += $temp;
            $f->{'pc'} = $temp;
            $count++;
        }
    }

    return (($count == 0) ? 0 : ($pc / $count));
}






#==============================================================================
#
# rall(neuron [, min, max, which]) - average Rall value. (brute-force: SLOW)
#
#==============================================================================

sub rall
{
    # get parameters
    my @arr = @{$_[0]};
    my $min = (defined($_[1]) ? $_[1] : 0);
    my $max = (defined($_[2]) ? $_[2] : maxBranchOrder(\@arr));
    my $wch = (defined($_[3]) ? $_[3] : "order");

    my $rll = 0;
    my $count = 0;

    foreach my $f (@FORKLIST)
    {
        my $fork = $f->{'number'};
        if (($arr[$fork]->{$wch} >= $min) &&
            ($arr[$fork]->{$wch} <= $max))
        {
            # get diameters
            my $dp = $arr[$fork]->{'d'};
            my $da = $arr[$fork]->{'kids'}[0]->{'d'};
            my $db = $arr[$fork]->{'kids'}[1]->{'d'};

            # try to find minimum of (Dp^v = Da^v + Db^v) between 0 and 5
            my $minimum = 65000;

            for (my $v = 1; $v <= 25000; $v++)
            {
                my $div = $v / 5000;
                if ((($dp ** $div) - ($da ** $div) - ($db ** $div)) < $minimum)
                {
                    $minimum = $div;
                }
            }
            if ($minimum < 5)
            {
                $rll += $minimum;
                $f->{'rall'} = $minimum;
                $count++;
            }
            else
            {
                $f->{'rall'} = -1;
            }
        }
    }

    return (($count == 0) ? 0 : ($rll / $count));
}