[system] / trunk / pg / lib / Hermite.pm Repository:
ViewVC logotype

View of /trunk/pg/lib/Hermite.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1079 - (download) (as text) (annotate)
Mon Jun 9 17:36:12 2003 UTC (15 years, 6 months ago) by apizer
File size: 11526 byte(s)
removed unneccesary shebang lines

    1 
    2 
    3 use strict;
    4 use Carp;
    5 
    6 =head1 NAME
    7 
    8   Hermite.pm
    9 
   10 =head1 SYNPOSIS
   11 
   12       Usage:
   13                 $obj = new Hermit(\@x_values, \y_valuses \@yp_values);
   14 
   15     #get and set methods
   16                 $ra_x_values = $obj -> ra_x(\@x_values);
   17     $ra_y_values = $obj -> ra_y;
   18     $ra_yp_values = $obj -> ra_yp;
   19 
   20     $obj -> initialize;           # calculates the approximation
   21 
   22     #get methods
   23     $rf_function                  = $obj -> rf_f;
   24     $rf_function_derivative       = $obj -> rf_fp;
   25     $rf_function_2nd_derivative   = $obj -> rf_fpp;
   26 
   27     $rh_critical_points           =$obj -> rh_critical_points
   28     $rh_inflection_points         =$obj -> rh_inflection_points
   29 
   30 
   31 
   32 
   33 =head1 DESCRIPTION
   34 
   35 This module defines an object containing a Hermite spline approximation to a function.
   36   The approximation
   37 consists of a piecewise cubic polynomial which agrees with the original
   38 function and its first derivative at
   39 the node points.
   40 
   41 This is useful for creating on the fly graphics.  Care must be taken to use a small
   42 number of points spaced reasonably far apart, preferably
   43 points with alternating slopes, since this will minimize
   44 the number of extraneous critical points
   45 introduced.  Too many points will introduce many small variations and
   46 a large number of extraneous critical points.
   47 
   48 There are even more extraneous inflections points.
   49 This parameter is probably not very useful.   A different
   50 approximation scheme needs to be used.
   51 
   52 =cut
   53 
   54 
   55 # an object that contains a cubic hermite spline
   56 
   57 package Hermite;
   58 
   59 
   60 sub new {
   61     my $class = shift;
   62     my ($ra_x,$ra_y,$ra_yp) = @_;
   63     my $self = {
   64            ra_x       =>  [],
   65      ra_y       =>  [],
   66      ra_yp        =>  [],
   67      rf_f       =>  0,   # refers to a subroutine for the function
   68      rf_fp        =>  0,   # refers to a subroutine for the derivative of the function
   69      rf_fpp       =>  0,   # refers to a subroutine for the second derivative of the function
   70      rh_critical_points   =>  {},
   71      rh_inflection_points   =>  {},
   72      rh_maximum_points    =>  {},
   73      rh_minimum_points    =>  {},
   74     };
   75     bless $self, $class;
   76     $self->define($ra_x,$ra_y,$ra_yp) if ref($ra_x) eq 'ARRAY' and ref($ra_y) eq 'ARRAY' and ref($ra_yp) eq 'ARRAY';
   77     $self;
   78 }
   79 
   80 sub define{
   81     my $self   = shift;
   82     my $ra_x  = shift;
   83     my $ra_y  = shift;
   84     my $ra_yp = shift;
   85     $self->ra_x($ra_x);
   86     $self->ra_y($ra_y);
   87     $self->ra_yp($ra_yp);
   88 
   89     $self -> initialize();
   90     ($self->{x_ref},$self->{y_ref}, $self->{yp_ref} )
   91 }
   92 
   93 sub initialize {
   94     my $self = shift;
   95     # define the functions rf_f
   96 
   97     $self -> {rf_f} = hermite_spline($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} );
   98 #     # define the function  rf_fp
   99      $self -> {rf_fp} = hermite_spline_p($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} );
  100 #     # define the function  rf_fp
  101      $self -> {rf_fpp} = hermite_spline_pp($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} );
  102 #     # define the critical points
  103      %{$self -> {rh_critical_points}} = critical_points($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp},
  104                                                      $self -> { rf_f },
  105                  );
  106 #     # define the inflection_points
  107      %{$self -> {rh_inflection_points}} = inflection_points($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} ,
  108                                                      $self -> { rf_f },
  109                  );
  110 #     # define the maximum points
  111 
  112     # define the minimum points
  113 
  114 }
  115 
  116 # fix up these accesses to do more checking
  117 sub ra_x {
  118     my $self = shift;
  119     my $ra_x  = shift;
  120     @{$self->{ra_x}} = @$ra_x if ref($ra_x) eq 'ARRAY';
  121     $self->{ra_x};
  122 
  123 }
  124 
  125 sub ra_y {
  126     my $self = shift;
  127     my $ra_y  = shift;
  128     @{$self->{ra_y}} = @$ra_y if ref($ra_y) eq 'ARRAY';
  129     $self->{ra_y};
  130 }
  131 
  132 sub ra_yp {
  133     my $self = shift;
  134     my $ra_yp  = shift;
  135     @{$self->{ra_yp}} = @$ra_yp if ref($ra_yp) eq 'ARRAY';
  136     $self->{ra_yp};
  137 }
  138 
  139 sub rf_f {
  140     my $self = shift;
  141     my $rf_f =shift;
  142     $self ->{rf_f} = $rf_f if defined( $rf_f );
  143     $self ->{rf_f};
  144 }
  145 sub rf_fp {
  146     my $self = shift;
  147     my $rf_fp =shift;
  148     $self ->{rf_fp} = $rf_fp if defined( $rf_fp );
  149     $self ->{rf_fp};
  150 }
  151 sub rf_fpp {
  152     my $self = shift;
  153     my $rf_fpp =shift;
  154     $self ->{rf_fpp} = $rf_fpp if defined( $rf_fpp );
  155     $self ->{rf_fpp};
  156 }
  157 
  158 sub rh_critical_points {
  159     my $self = shift;
  160     $self -> { rh_critical_points};
  161 }
  162 
  163 sub rh_inflection_points {
  164     my $self = shift;
  165     $self -> { rh_inflection_points};
  166 }
  167 
  168 ##Internal subroutines
  169 
  170 
  171 sub critical_points{
  172     my $ra_x = shift;
  173     my $ra_y =shift;
  174     my $ra_yp = shift;
  175     my $rf_hermite_fun =shift;
  176     my %critical_points = ();
  177     my $last_index = @$ra_x -2;
  178     foreach my $i (0 .. $last_index ) {
  179         internal_critical_points($ra_x->[$i],   $ra_y->[$i],   $ra_yp->[$i],
  180                            $ra_x->[$i+1], $ra_y->[$i+1], $ra_yp->[$i+1],
  181          \%critical_points,$rf_hermite_fun);
  182     }
  183     %critical_points;
  184 }
  185 
  186 sub inflection_points{
  187     my $ra_x = shift;
  188     my $ra_y =shift;
  189     my $ra_yp = shift;
  190     my $rf_hermite_fun =shift;
  191     my %inflection_points = ();
  192     my $last_index = @$ra_x -2;
  193     foreach my $i (0 .. $last_index ) {
  194         internal_inflection_points($ra_x->[$i],   $ra_y->[$i],   $ra_yp->[$i],
  195                            $ra_x->[$i+1], $ra_y->[$i+1], $ra_yp->[$i+1],
  196          \%inflection_points,$rf_hermite_fun);
  197     }
  198     %inflection_points;
  199 }
  200 sub internal_critical_points{
  201     my ($x0,$l,$lp, $x1,$r,$rp, $rh_roots ,$rf_function) = @_;
  202      #data for one segment of the hermite spline
  203 
  204      # coefficients for the approximating polynomial
  205 
  206      my @a = (   $l,
  207                   $lp,
  208       -$lp/2 + (3*(-$lp - 2*($l - $r) - $rp))/2 + $rp/2,
  209       $lp + 2*($l - $r) + $rp
  210                );
  211 
  212     my ($root1, $root2, $z1,$z2);
  213     if ( $a[3] == 0 ) {
  214   if ( $a[2] == 0 ) {
  215   } else {
  216       $root1 = -$a[1]/( 2*$a[2] );
  217       if ( 0 <= $root1 and $root1 < 1) {
  218     $z1 = $root1*($x1-$x0) + $x0;
  219 
  220           $rh_roots -> {$z1} = &$rf_function($z1);
  221       }
  222   }
  223    } else {
  224   my $discriminent = (4*$a[2]**2 - 12*$a[1]*$a[3]);
  225 
  226 
  227   if ( $discriminent >= 0 ) {
  228       $discriminent = $discriminent**0.5;
  229       $root1 = (-2*$a[2] - $discriminent )/( 6*$a[3] );
  230       $root2 = (-2*$a[2] + $discriminent )/( 6*$a[3] );
  231       $z1 = $root1*($x1-$x0) + $x0;
  232       $z2 = $root2*($x1-$x0) + $x0;
  233       $rh_roots -> {$z1} = &$rf_function($z1) if  0 <= $root1 and $root1 < 1;
  234       $rh_roots -> {$z2} = &$rf_function($z1) if  0 <= $root2 and $root2 < 1;
  235   }
  236     }
  237 }
  238 
  239 sub internal_inflection_points {
  240     my ($x0,$l,$lp,,$x1,$r,$rp,$rh_roots,$rf_function) = @_;
  241      #data for one segment of the hermite spline
  242 
  243      # coefficients for the approximating polynomial
  244 
  245      my @a = (   $l,
  246                   $lp,
  247       -$lp/2 + (3*(-$lp - 2*($l - $r) - $rp))/2 + $rp/2,
  248       $lp + 2*($l - $r) + $rp
  249                );
  250     if ($a[3] == 0 ) {
  251     } else {
  252   my $root1 = -$a[2]/( 3*$a[3] );
  253   my $z1 = $root1*($x1-$x0) + $x0;
  254   $rh_roots -> {$z1} = &$rf_function($z1) if  0 <= $root1 and $root1 < 1;
  255     }
  256 
  257 }
  258 
  259 
  260 
  261 sub cubic_hermite {
  262     my ( $x0, $y0,$yp0, $x1, $y1, $yp1 ) = @_;
  263     my @a;
  264     my $width = $x1 - $x0;
  265     $yp0 = $yp0*$width;  # normalize to unit width
  266     $yp1 = $yp1*$width;
  267 
  268     $a[0] = $y0;
  269     $a[1] = $yp0;
  270     $a[2] = -3*$y0 - 2*$yp0 +3*$y1 -$yp1;
  271     $a[3] = 2*$y0 + $yp0 - 2*$y1 +$yp1;
  272 
  273     my $f = sub {
  274                         my $x = shift;
  275                         #normalize to unit width
  276       $x = ( $x - $x0 )/$width;
  277       ( ($a[3]*$x + $a[2]) * $x + $a[1] )*$x + $a[0];
  278 
  279                 };
  280     $f;
  281 }
  282 
  283 sub cubic_hermite_p {
  284     my ( $x0, $y0,$yp0, $x1, $y1, $yp1 )=@_;
  285     my @a;
  286     my $width = $x1 - $x0;
  287     $yp0 = $yp0*$width;  # normalize to unit width
  288     $yp1 = $yp1*$width;
  289 
  290     $a[0] = $y0;
  291     $a[1] = $yp0;
  292     $a[2] = -3*$y0 - 2*$yp0 +3*$y1 -$yp1;
  293     $a[3] = 2*$y0 + $yp0 - 2*$y1 +$yp1;
  294 
  295     my $fp = sub {
  296                         my $x = shift;
  297                         #normalize to unit width
  298       $x = ( $x - $x0 )/$width;
  299       ( (3*$a[3]*$x + 2*$a[2]) * $x + $a[1] )/$width ;
  300                 };
  301 
  302 
  303 
  304      $fp;
  305 
  306 }
  307 
  308 sub cubic_hermite_pp {
  309     my ( $x0, $y0,$yp0, $x1, $y1, $yp1 ) = @_;
  310     my @a;
  311     my $width = $x1 - $x0;
  312     $yp0 = $yp0*$width;  # normalize to unit width
  313     $yp1 = $yp1*$width;
  314 
  315     $a[0] = $y0;
  316     $a[1] = $yp0;
  317     $a[2] = -3*$y0 - 2*$yp0 +3*$y1 -$yp1;
  318     $a[3] = 2*$y0 + $yp0 - 2*$y1 +$yp1;
  319 
  320 
  321 
  322     my $fpp = sub {
  323                         my $x = shift;
  324                         #normalize to unit width
  325       $x = ( $x - $x0 )/$width;
  326        (6*$a[3]*$x + 2*$a[2])/$width**2 ;
  327                 };
  328 
  329     $fpp;
  330 }
  331 
  332 
  333 
  334 sub hermite_spline {
  335   my ($xref, $yref, $ypref) = @_;
  336   my @xvals  = @$xref;
  337   my @yvals  = @$yref;
  338   my @ypvals = @$ypref;
  339   my $x0 = shift @xvals;
  340   my $y0 = shift @yvals;
  341   my $yp0 = shift @ypvals;
  342   my ($x1,$y1,$yp1);
  343   my @polys;  #calculate a hermite polynomial evaluator for each region
  344 
  345   while (@xvals) {
  346     $x1 = shift @xvals;
  347     $y1 = shift @yvals;
  348     $yp1 = shift @ypvals;
  349 
  350     push @polys, cubic_hermite($x0, $y0, $yp0 , $x1, $y1 , $yp1);
  351     $x0  = $x1;
  352     $y0  = $y1;
  353     $yp0 = $yp1;
  354   }
  355 
  356 
  357   my $hermite_spline_function = sub {
  358     my $x = shift;
  359     my $y;
  360     my $fun;
  361     my @xvals = @$xref;
  362     my @fns = @polys;
  363     return &{$fns[0]} ($x) if $x == $xvals[0]; #handle left most endpoint
  364 
  365     while (@xvals && $x > $xvals[0]) {  # find the function for this range of x
  366       shift(@xvals);
  367       $fun = shift(@fns);
  368     }
  369 
  370     # now that we have the left hand of the input
  371     #check first that x isn't out of range to the left or right
  372     if (@xvals  && defined($fun) )  {
  373       $y =&$fun($x);
  374     }
  375     $y;
  376   };
  377   $hermite_spline_function;
  378 }
  379 
  380 sub hermite_spline_p {
  381   my ($xref, $yref, $ypref) = @_;
  382   my @xvals  = @$xref;
  383   my @yvals  = @$yref;
  384   my @ypvals = @$ypref;
  385   my $x0 = shift @xvals;
  386   my $y0 = shift @yvals;
  387   my $yp0 = shift @ypvals;
  388   my ($x1,$y1,$yp1);
  389   my @polys;  #calculate a hermite polynomial evaluator for each region
  390   while (@xvals) {
  391     $x1 = shift @xvals;
  392     $y1 = shift @yvals;
  393     $yp1 = shift @ypvals;
  394     push @polys, cubic_hermite_p($x0, $y0, $yp0 , $x1, $y1 , $yp1);
  395     $x0  = $x1;
  396     $y0  = $y1;
  397     $yp0 = $yp1;
  398   }
  399 
  400 
  401   my $hermite_spline_function_p = sub {
  402     my $x = shift;
  403     my $y;
  404     my $fun;
  405     my @xvals = @$xref;
  406     my @fns = @polys;
  407     return $y=&{$fns[0]} ($x) if $x == $xvals[0]; #handle left most endpoint
  408 
  409     while (@xvals && $x > $xvals[0]) {  # find the function for this range of x
  410       shift(@xvals);
  411       $fun = shift(@fns);
  412     }
  413 
  414     # now that we have the left hand of the input
  415     #check first that x isn't out of range to the left or right
  416     if (@xvals  && defined($fun) )  {
  417       $y =&$fun($x);
  418     }
  419     $y;
  420   };
  421   $hermite_spline_function_p;
  422 }
  423 sub hermite_spline_pp {
  424   my ($xref, $yref, $ypref) = @_;
  425   my @xvals  = @$xref;
  426   my @yvals  = @$yref;
  427   my @ypvals = @$ypref;
  428   my $x0 = shift @xvals;
  429   my $y0 = shift @yvals;
  430   my $yp0 = shift @ypvals;
  431   my ($x1,$y1,$yp1);
  432   my @polys;  #calculate a hermite polynomial evaluator for each region
  433   while (@xvals) {
  434     $x1 = shift @xvals;
  435     $y1 = shift @yvals;
  436     $yp1 = shift @ypvals;
  437     push @polys, cubic_hermite_pp($x0, $y0, $yp0 , $x1, $y1 , $yp1);
  438     $x0  = $x1;
  439     $y0  = $y1;
  440     $yp0 = $yp1;
  441   }
  442 
  443 
  444   my $hermite_spline_function_pp = sub {
  445     my $x = shift;
  446     my $y;
  447     my $fun;
  448     my @xvals = @$xref;
  449     my @fns = @polys;
  450     return $y=&{$fns[0]} ($x) if $x == $xvals[0]; #handle left most endpoint
  451 
  452     while (@xvals && $x > $xvals[0]) {  # find the function for this range of x
  453       shift(@xvals);
  454       $fun = shift(@fns);
  455     }
  456 
  457     # now that we have the left hand of the input
  458     #check first that x isn't out of range to the left or right
  459     if (@xvals  && defined($fun) )  {
  460       $y =&$fun($x);
  461     }
  462     $y;
  463   };
  464   $hermite_spline_function_pp;
  465 }
  466 
  467 
  468 1;

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9