Math::Project3D - Project functions of multiple parameters

  use Math::Project3D;

  my $projection = Math::Project3D->new(
    plane_basis_vector => [0,  0, 0],
    plane_direction1   => [.4, 1, 0],
    plane_direction2   => [.4, 0, 1],
    projection_vector  => [1,  1, 1], # defaults to normal of the plane
  );

  $projection->new_function(
    'u,v', 'sin($u)', 'cos($v)', '$u' 
  );

  # Rotate the points before projecting them.
  # Rotate every point the same way we need to rotate the
  # z-axis to get the x-axis.
  $projection->rotate([1,0,0]);

  # Nah, changed my mind
  $projection->unrotate();

  ($plane_coeff1, $plane_coeff2, $distance_coeff) =
     $projection->project( $u, $v );

  basis point of the plane + i2 * d + i3 * e

   n1*i1 + d1*i2 + e1*i3 = p1 + x(t)
   n2*i1 + d2*i2 + e2*i3 = p2 + y(t)
   n3*i1 + d3*i2 + e3*i3 = p3 + z(t)

# See the POD documentation at the end of this
# document for detailed copyright information.
# (c) 2002-2006 Steffen Mueller, all rights reserved.

package Math::Project3D;

use strict;
use warnings;

use 5.006;

use vars qw/$VERSION/;

$VERSION = 1.02;

use Carp;

use Math::MatrixReal;
use Math::Project3D::Function;

    sub _max {
        my $max = 0;
        for (@_) {
            $max = $_ if $_ > $max;
        }
        $max
    }
    sub _new_from_rows {
        my $class = 'Math::MatrixReal';
        my $rows = shift;
        my @rows;
        foreach (@$rows) {
            push @rows,
              (ref($_) eq 'Math::MatrixReal' ? [@{$_->[0][0]}] : [@{$_}]);
        }
        return bless [
            \@rows,
            scalar(@rows),
            _max(
                map {
                    ref($_) eq 'Math::MatrixReal' ? $_->[2] : scalar(@$_)
                }
                @rows
            ),
            undef,
            undef,
            undef
        ] => $class;
    }

    sub _new_from_cols {
        my $class = 'Math::MatrixReal';
        my $cols = shift;
        my $num_rows = _max(
            map {
                ref($_) eq 'Math::MatrixReal' ? $_->[1] : scalar(@$_)
            }
            @$cols
        );
        my $num_cols = @$cols;
        
        my @rows;
        my $cn = 0;
        foreach (@$cols) {
            my $col = $_;
            if (ref($col) eq 'Math::MatrixReal') {
                $col = $col->[0];
                $rows[$_][$cn] = $col->[$_][0] for 0..$num_rows-1;
            }
            else {
                $rows[$_][$cn] = $col->[$_] for 0..$num_rows-1;
            }
            $cn++;
        }
        return bless [
            \@rows,
            $num_rows,
            $num_cols,
            undef,
            undef,
            undef
        ] => $class;
    }


# class and object method new_function
# 
# Wrapper around Math::Project3D->new()
# Returns compiled function (anon sub)
# As a class method, it does not have side effects.
# As an object method, it assigns the returned function to
# the object's function attribute.

sub new_function {
   @_ > 1 or croak "No arguments supplied to 'new_function' method.";

   my $self = shift;

   my @components = @_;

   my $function = Math::Project3D::Function->new(@components);

   if (ref $self eq __PACKAGE__) {
      $self->{function} = $function;
   }

   return $function;
}


# Class and object method/constructor new
# 
# Arguments are used in the object's anon hash.
# Creates new object.
# Returns MP3D instance.

sub new {
   my $proto = shift;
   my $class = ref $proto || $proto;

   # !!!FIXME!!! Insert argument checking here
   my $self  = {
     function => undef,
     @_
   };

   bless $self => $class;

   # Some attributes are required.
   my $missing_attribute = $self->_require_attributes(qw(
                                   plane_basis_vector
                                   plane_direction1
                                   plane_direction2
                           ));
   croak "Required attribute '$missing_attribute' missing."
     if defined $missing_attribute;

   # Transform the vector/matrix specifications ((nested) array refs)
   # into Math::MatrixReal's.
   foreach ( qw(
        plane_basis_vector
        plane_direction1
        plane_direction2
      ) ) {

      # Argument checking, we either want an array ref or a MMR object.
      croak "Invalid argument '$_' passed to new. Should be an array reference."
        if not exists $self->{$_};
      
      next if ref $self->{$_} eq 'Math::MatrixReal';
      croak "Invalid argument '$_' passed to new. Should be an array reference."
        if ref $self->{$_} ne 'ARRAY';

      # Transform into an MMR object.
      $self->{$_} = $self->_make_matrix($self->{$_});
   }

   # Projection vector defaults to the normal vector of the plane,
   # but may also be specified explicitly as array ref or MMR object.
   if ( defined $self->{projection_vector} and
        ref $self->{projection_vector} eq 'ARRAY' ) {

      $self->{projection_vector} = $self->_make_matrix( $self->{projection_vector} );

   } elsif ( not defined $self->{projection_vector} or
             ref $self->{projection_vector} ne 'Math::MatrixReal' ) {

      # Defaults to the normal of the plane.
      $self->{projection_vector} = $self->{plane_direction1}->vector_product(
                                     $self->{plane_direction2}
                                   );
   }

   # Now generate the linear equation system that has to be solved by
   # MMR in order to do project a point onto the plane.
   # 
   # MMR does all the dirty work. (Thanks!) So we just create a matrix.
   # (d, e be the directional vectors, n the vector we project along,
   #  p the basis point. Let i be the solution)
   # 
   # x(t) + n1*i1 = d1*i2 + e1*i3 + p1
   # y(t) + n2*i1 = d2*i2 + e2*i3 + p2
   # z(t) + n3*i1 = d3*i2 + e3*i3 + p3
   # 
   # This is the intersection of the plane and the corresponding orthogonal
   # line through s(t). Another form:
   # 
   # n1*i1 + d1*i2 + e1*i3 = p1 + x(t)
   # n2*i1 + d2*i2 + e2*i3 = p2 + y(t)
   # n3*i1 + d3*i2 + e3*i3 = p3 + z(t)
   # 
   # linear_eq_system will be the matrix of n,d,e.
   # result_vector will be p+s(t) (later)

   my $linear_eq_system = _new_from_cols(
                         [
                           $self->{projection_vector},
                           $self->{plane_direction1},
                           $self->{plane_direction2},
                         ]
   );

   # Now we do the first step towards solving the equation system.
   # This does not have to be repeated for every projection.
   $self->{lr_matrix} = $linear_eq_system->decompose_LR();

   return $self;
}


# Method project
# 
# Does the projection calculations for a single set of function
# parameters.
# Takes the function parameters as argument.
# Returns undef on failure (plane and P + lambda*projection_vector do
# not intersect!). On success:
# Returns the coefficients of the point on the plane, 
# the distance of the source point from the plane in lengths
# of the projection vector.

sub project {
   my $self  = shift;

   croak "You need to assign a function before projecting it."
     if not defined $self->get_function();

   # Apply function
   my $point = _new_from_cols(
                 [
                   [ $self->{function}->(@_) ],
                 ]
   );

   # Generate result_vector
   my $result_vector = _new_from_cols(
                         [
                           $self->{plane_basis_vector} + $point
                         ]
   );

   # Solve the l_e_s.
   my ($dimension, $p_vector, undef) = $self->{lr_matrix}->solve_LR($result_vector);
   
   # Did we find a solution?
   return undef if not defined $dimension;

   # $dimension == 0 => one solution
   # $dimension == 1 => straight line
   # $dimension == 2 => plane (impossible :) )
   # ...
   # $dimension == 1 is possible (point and projection_vector part
   # of the plane). Hence: !!!FIXME!!!

   return $p_vector->element(2,1), # coefficient 1
          $p_vector->element(3,1), # coefficient 2
          $p_vector->element(1,1); # distance in lengths of the projection vector
}


# Method project_range_callback
# 
# Does the projection calculations for ranges of function parameters.
# Takes an code reference and a number array refs of ranges as argument.
# Ranges are specified using three values:
# [lower_bound, upper_bound, increment].
# Alternatively, one may specify only one element in any of the array
# references. This will then be a static parameter.
# The number of ranges (array references) corresponds to the number of
# parameters to the vectorial function.
# The callback is called for every set of result values with an array
# reference of parameters and the list of the three result values.

sub project_range_callback {
   my  $self     = shift;

   croak "You need to assign a function before projecting it."
     if not defined $self->get_function();

   # Argument checking
   my  $callback = shift;
   ref $callback eq 'CODE' or
     croak "Invalid code reference passed to project_range_callback.";

   my @ranges    = @_;

   croak "Invalid range parameters passed to project_range_callback"
     if grep { ref $_ ne 'ARRAY' } @ranges;


   my $function = $self->get_function();

   # Replace all ranges that consist of a single number with ranges
   # that iterate from that number to itself.
   foreach (@ranges) {
      $_ = [$_->[0], $_->[0], 1] if @$_ == 1;
   }

   # Calculate every range's length and store it in @lengths.
   my @lengths   = map {
                         #      upper  - lower    / increment
                         int( ($_->[1] - $_->[0]) / $_->[2] ),
                       } @ranges;

   # Prepare counters for every range.
   my @counters  = (0) x scalar(@ranges);

   # Calculate the number if iterations needed.
   # It is $_+1 and not $_ because the lengths are the lengths we
   # need for the comparisons inside the long for loop. We save one
   # op in there that way.
   my $iterations = 1;
   $iterations *= ( $_ + 1 ) for @lengths;
   
   # For all possible combinations of parameters...
   for (my $i = 1; $i <= $iterations; $i++) {
      
      # Get current function parameters
      my @params;

      # Get one parameter for every range
      for (my $range_no = 0; $range_no < @ranges; $range_no++) {
         # lower + increment * current_count
         push @params, $ranges[$range_no][0] +
                       $ranges[$range_no][2] * $counters[$range_no];
      }
      
      # Increment outermost range by one. If it got out of bounds,
      # make it 0 and increment the next range, etc.
      my $j = 0;
      while (defined $ranges[$j] and ++$counters[$j] > $lengths[$j]) {
         $counters[$j] = 0;
         $j++;
      }

      # Apply function
      my $point = _new_from_cols(
                    [
                      [ $function->(@params) ],
                    ]
      );

      # Generate result_vector
      my $result_vector = _new_from_cols(
                            [
                              $self->{plane_basis_vector} + $point
                            ]
      );

      # Solve the l_e_s.
      my ($dimension, $p_vector, undef) = $self->{lr_matrix}->solve_LR($result_vector);
   
      # Did we find a solution?
      croak "Could not project $result_vector."
        if not defined $dimension;

      # $dimension == 0 => one solution
      # $dimension == 1 => straight line
      # $dimension == 2 => plane (impossible :) )
      # ...
      # $dimension == 1 is possible (point and projection_vector part
      # of the plane). Hence: !!!FIXME!!!
      
      $callback->(
             $p_vector->element(2,1), # coefficient 1
             $p_vector->element(3,1), # coefficient 2
             $p_vector->element(1,1), # distance in lengths of the projection vector
             $j,                      # how many ranges did we increment?
      );
   }
   
   return();
}


# Method project_list
# 
# Wrapper around project(), therefore slow.
# Takes a list of array refs as argument. The array refs
# are to contain sets of function parameters.
# Calculates every set of parameters' projection and stores the
# three associated values (coefficients and distance coefficient)
# in an n*3 matrix (MMR obj) wheren is the number of points
# projected.
# Currently croaks if any of the points cannot be projected onto
# the plane using the given projection vector. (In R3 -> R2, try
# using the plane's normal vector which is guaranteed not to be
# parallel to the plane.)
# Returns that MMR object.

sub project_list {
   my $self = shift;
   croak "No arguments passed to project_list()."
     if @_ == 0;

   croak "You need to assign a function before projecting it."
     if not defined $self->get_function();

   # Create result matrix to hold individual results.
   my $result_matrix = Math::MatrixReal->new(scalar(@_), 3);

   # Keep track of the matrix row
   my $result_no = 1;

   foreach my $array_ref (@_) {
      my ($coeff1, $coeff2, $dist) = $self->project(@$array_ref);

      croak "Vector $result_no cannot be projected."
        if not defined $coeff1;

      # Assign results
      $result_matrix->assign($result_no, 1, $coeff1);
      $result_matrix->assign($result_no, 2, $coeff2);
      $result_matrix->assign($result_no, 3, $dist);
      $result_no++;
   }

   return $result_matrix;
}


# Accessor get_function
# 
# No parameters.
# Returns the current function code ref.

sub get_function {
   my $self = shift;
   return $self->{function};
}


# Accessor set_function
# 
# Takes a code ref as argument.
# Sets the object's function to that code ref.
# Returns the code ref.

sub set_function {
   my $self     = shift;
   my $function = shift;

   ref $function eq 'CODE' or croak "Argument to set_function must be code reference.";

   $self->{function} = $function;

   return $self->{function};
}


# Private method _require_attributes
# 
# Arguments must be a list of attribute names (strings).
# Tests for the existance of those attributes.
# Returns the missing attribute on failure, undef on success.

sub _require_attributes {
   my $self = shift;
   exists $self->{$_} or return $_ foreach @_;
   return undef;
}


# Private method _make_matrix
# 
# Takes a list of array refs as arguments.
# Creates a Math::MatrixReal object from the arrays.
# Returns the Math::MatrixReal object.

sub _make_matrix {
   my $self = shift;
   @_ or croak "No arguments passed to _make_matrix.";

   my $matrix = Math::MatrixReal->new_from_cols( [ @_ ] );

   return $matrix;
}


# Evil method rotate
# 
# Takes a vector as argument. (Either an MMR object or an array ref
# of an array containing three components.)
# The method will replace the function with a wrapper around the original
# function that rotates the result of the original function by the same
# angles that the z-axis (e3) needs to be turned to become the passed vector.
# Example: Passed [1,0,0] means that e3 will be rotated to become e1.
# Hence all points will be rotated by 90 degrees and orthogonally to e2.
# Returns the wrapper function and the old function as code references.
# 
# Important note: You can apply this function multiple times. That means
# if you rotated the function once, you may do so again and the original
# rotated function will be wrapped again so that you effectively rotate
# it twice. Note, however, that you will sacrifice performance on the
# altar of recursion that way.
# 
# !!!FIXME!!! This is slow. Sloooooow. Conceptually slow, but the
# implementation sucks badly, too. It is especially unnerving to me that
# the whole lexical scope of a myriad of intermediate result matrices and
# vectors is kept because we use closures. Closures are great, but not
# in bloated lexical scopes. How fix that without introducing either
#  - an ugly additional method
#  - ugly additional blocks to keep the scope clean. (Yuck!)

sub rotate {
   my $self = shift;

   croak "You need to assign a function before rotating it."
     if not defined $self->get_function();

   # We want to rotate everything the same way we would
   # have to rotate e3 (z unit vector) to become the unit
   # vector parallel to $e3_.
   my $e3_ = shift;

   # Make sure we work with a MMR vector
   if (ref $e3_ ne 'Math::MatrixReal') {
      croak "Invalid vector passed to rotate()."
        if ref $e3_ ne 'ARRAY';

      $e3_ = _new_from_cols([$e3_]);
   }

   $e3_ *= 1 / $e3_->length();

   my $e3 = _new_from_cols([[0,0,1]]);

   # The axis we want to rotate around
   my $axis = $e3->vector_product($e3_);

   # The angle we want to rotate by.
   my $angle = acos($e3->scalar_product($e3_));

   # Rotationsmatrix um Achse (Vektor) a = (a,b,c):
   #                           [a*a a*b a*c]
   # M(rot) = (1-cos(alpha)) * [a*b b*b b*c]  +
   #                           [a*c b*c c*c]
   #                           [1   0   0  ]
   #          cos(alpha)     * [0   1   0  ]  +
   #                           [0   0   1  ]
   #                           [0   -c  b  ]
   #          sin(alpha)     * [c   0   -a ]
   #                           [-b  a   0  ]
   # (from: Merziger, Wirth: "Repetitorium der Hoeheren Mathematik"
   #        Binomi Verlag. ISBN 3-923923-33-3)

   my ($a, $b, $c) = ($axis->element(1,1), $axis->element(2,1), $axis->element(3,1));

   my $matrix1 = _new_from_cols(
      [
        [$a*$a, $a*$b, $a*$c],
        [$b*$a, $b*$b, $b*$c],
        [$c*$a, $c*$b, $c*$c],
      ],
   );

   my $matrix2 = Math::MatrixReal->new_diag([1,1,1]);

   my $matrix3 = _new_from_cols(
      [
        [0,   -$c, $b ],
        [$c,  0,   -$a],
        [-$b, $a,  0  ],
      ],
   );

   # What a painful birth, but here we are:
   my $rot_matrix = (1-cos($angle)) * $matrix1 + cos($angle) * $matrix2 +
                    sin($angle)     * $matrix3;

   # For (Rotated vector r=(r1,r2,r3), source vector x=(x1,x2,x3)):
   # (eq1) r = $rot_matrix * x
   # Hence:
   # (eq2) x = transposed($rot_matrix) * r
   # 
   # x and $rot_matrix (and therefore also transposed($rot_matrix)
   # are know. Hance, (eq2) can be solved as a linear equation system
   # of the form: ($rot_matrix be all elements a(ij))
   # 
   # a11*r1 + a12*r2 + a13*r3 = x1
   # a21*r1 + a22*r2 + a23*r3 = x2
   # a31*r1 + a32*r2 + a33*r3 = x3

   # Transpose the rotation matrix and then decompose_LR
   $rot_matrix->transpose($rot_matrix);
   $rot_matrix = $rot_matrix->decompose_LR();

   # Save old function
   my $old_function = $self->get_function;

   my $rotator = sub {

      # If the first argument is 'restore', we restore the old
      # function from the lexical $old_function and return.
      # This is needed for the unrotate() method.
      $self->{function} = $old_function, return if $_[0] eq 'restore';

      # We apply the old function as usual.
      my $source_vec = _new_from_cols(
         [ [ $old_function->(@_) ] ],
      );

      # But then we rotate the results.
      my ($dimension, $result_vector, undef) = $rot_matrix->solve_LR($source_vec);

      # There should be a solution to the equation system.
      croak "Dimension of rotation result is '$dimension' but should be 0."
        if $dimension;

      # return the rotated vector
      return (
        $result_vector->element(1,1),
        $result_vector->element(2,1),
        $result_vector->element(3,1),
      );
   };
   
   # Set the new function wrapper to become the actual function
   # that will be used for projection calculations.
   $self->set_function($rotator);

   # The level of nested rotation has been increased.
   $self->{rotated}++;

   return ($rotator, $old_function);
}


# method unrotate
# 
# Removes the evil hack of rotation using an evil hack.
# Takes an optional integer as argument. Removes
# [integer] levels of rotation. 'Level of rotation'
# means: one rotation wrapper of the original function
# as wrapped using rotate().
# If no integer was passed, defaults to the total number
# of rotations that were made, effectively removing
# any rotation.
# Returns the number of wrappers (rotations) removed.

sub unrotate {
   my $self = shift;

   croak "You need to assign a function before rotating it."
     if not defined $self->get_function();

   # How many levels of rotation to we want to remove?
   my $level = shift;

   # Boundary checking for the level. Default to the number
   # of rotation wrappers.
   if (not defined $level or $level <= 0 or $level > $self->{rotated}) {
      $level = $self->{rotated};
   }

   $level = int $level; # I don't trust users.

   # If we're rotated at all.
   if ($self->{rotated}) {

      # Count the number of wrappers removed.
      my $no_unrotated = 0;

      # While we still want to unrotate and while we still can
      while ($level != 0 and $self->{rotated} != 0) {

         # Let the function restore its parent.
         $self->{function}->('restore');

         $self->{rotated}--;
         $level--;
         $no_unrotated++;
      }

      return $no_unrotated;

   } else {
      # We aren't even rotated.
      return 0;

   }
}


# Function acos (arc cosine)
# Not a method.
sub acos { atan2( sqrt( 1 - $_[0] * $_[0] ), $_[0] ) }

1;

__END__