ICC-Profile
view release on metacpan or search on metacpan
lib/ICC/Support/spline.pm view on Meta::CPAN
push(@sol, $self->[1][0] + $sf * ($in - $self->[2][0])/$self->[3][0]);
}
# for each segment (0 <= s < 1)
for my $i (0 .. $n - 1) {
# if input is lower knot
if ($in == $self->[2][$i]) {
# add solution
push(@sol, $self->[1][0] + $i * $sf);
}
# if input lies within the y-value range
if ($in >= $self->[4][$i][0] && $in <= $self->[4][$i][1]) {
# compute inverse values (t-values)
@t = _rev($self, $in, $i);
# add solution(s)
push(@sol, map {$self->[1][0] + ($i + $_) * $sf} @t);
}
}
# if input is last knot (s == 1)
if ($in == $self->[2][-1]) {
# add solution
push(@sol, $self->[1][1]);
}
# if an extrapolated solution (s > 1)
if (($in > $self->[2][-1] && $self->[3][-1] > 0) || ($in < $self->[2][-1] && $self->[3][-1] < 0)) {
# add solution
push(@sol, $self->[1][1] + $sf * ($in - $self->[2][-1])/$self->[3][-1]);
}
# return result (array or first solution)
return(wantarray ? @sol : $sol[0]);
}
# compute curve derivative
# parameters: (input_value)
# returns: (derivative_value)
sub derivative {
# get parameters
my ($self, $in) = @_;
# if an extrapolated solution (s < 0)
if (($self->[1][1] > $self->[1][0] && $in < $self->[1][0]) || ($self->[1][1] < $self->[1][0] && $in > $self->[1][0])) {
# return endpoint derivative
return($self->[3][0] * $#{$self->[2]}/($self->[1][1] - $self->[1][0]));
# if an extrapolated solution (s >= 1)
} elsif (($self->[1][1] > $self->[1][0] && $in >= $self->[1][1]) || ($self->[1][1] < $self->[1][0] && $in <= $self->[1][1])) {
# return endpoint derivative
return($self->[3][-1] * $#{$self->[2]}/($self->[1][1] - $self->[1][0]));
} else {
# return interpolated derivative value
return(_derv($self, POSIX::modf($#{$self->[2]} * ($in - $self->[1][0])/($self->[1][1] - $self->[1][0]))) * $#{$self->[2]}/($self->[1][1] - $self->[1][0]));
}
}
# compute parametric partial derivatives
# parameters: (input_value)
# returns: (partial_derivative_array)
sub parametric {
# get parameters
my ($self, $in) = @_;
# local variables
my ($low, $t, $tc, $h00, $h01, $h10, $h11, @pj);
# update parametric derivative matrix, if necessary
_objpara($self) if (! defined($self->[5]) || $#{$self->[5]} != $#{$self->[2]});
# if an extrapolated solution (s < 0)
if (($self->[1][1] > $self->[1][0] && $in < $self->[1][0]) || ($self->[1][1] < $self->[1][0] && $in > $self->[1][0])) {
# for each output value (parameter)
for my $i (0 .. $#{$self->[2]}) {
# compute partial derivative
$pj[$i] = ($i == 0) + $self->[5][0][$i] * $#{$self->[2]} * ($in - $self->[1][0])/($self->[1][1] - $self->[1][0]);
}
# if an extrapolated solution (s >= 1)
} elsif (($self->[1][1] > $self->[1][0] && $in >= $self->[1][1]) || ($self->[1][1] < $self->[1][0] && $in <= $self->[1][1])) {
# for each output value (parameter)
for my $i (0 .. $#{$self->[2]}) {
# compute partial derivative
$pj[$i] = ($i == $#{$self->[2]}) + $self->[5][-1][$i] * $#{$self->[2]} * ($in - $self->[1][1])/($self->[1][1] - $self->[1][0]);
}
} else {
# compute position (0 - 1) and lower knot index
($t, $low) = POSIX::modf($#{$self->[2]} * ($in - $self->[1][0])/($self->[1][1] - $self->[1][0]));
# compute Hermite coefficients
$tc = 1 - $t;
$h00 = (1 + 2 * $t) * $tc * $tc;
$h01 = 1 - $h00;
$h10 = $t * $tc * $tc;
$h11 = -$t * $t * $tc;
# for each output value (parameter)
lib/ICC/Support/spline.pm view on Meta::CPAN
# return y-values
return(map {_fwd($self, POSIX::modf($_))} _minmax($self));
} else {
# error
croak("unsupported format for min/max values");
}
}
# make table for 'curv' objects
# assumes curve domain/range is (0 - 1)
# parameters: (number_of_table_entries, [direction])
# returns: (ref_to_table_array)
sub table {
# get parameters
my ($self, $n, $dir) = @_;
# local variables
my ($up, $table);
# validate number of table entries
($n == int($n) && $n >= 2) || carp('invalid number of table entries');
# purify direction flag
$dir = ($dir) ? 1 : 0;
# array upper index
$up = $n - 1;
# for each table entry
for my $i (0 .. $up) {
# compute table value
$table->[$i] = _transform($self, $dir, $i/$up);
}
# return table reference
return($table);
}
# make 'curv' object
# assumes curve domain/range is (0 - 1)
# parameters: (number_of_table_entries, [direction])
# returns: (ref_to_curv_object)
sub curv {
# return 'curv' object reference
return(ICC::Profile::curv->new(table(@_)));
}
# normalize transform
# adjusts object values linearly
# so that endpoint values are 0 or 1
# adjusts input-range and output-values by default
# adjusts input-values if format is 'x'
# adjusts output-values if format is 'y'
# parameters: ([format])
sub normalize {
# get parameters
my ($self, $fmt) = @_;
# local variables
my ($off, $range);
# verify flag parameter
(! defined($fmt) || $fmt eq 'x' || $fmt eq 'y') || croak('invalid normalize format parameter');
# if range selected
if (! defined($fmt) || $fmt eq 'x') {
# if range is increasing
if ($self->[1][1] > $self->[1][0]) {
# set range (0 - 1)
$self->[1][0] = 0;
$self->[1][1] = 1;
} else {
# set range (1 - 0)
$self->[1][0] = 1;
$self->[1][1] = 0;
}
}
# if output values selected
if (! defined($fmt) || $fmt eq 'y') {
# if output values are increasing
if ($self->[2][-1] > $self->[2][0]) {
# set offset
$off = $self->[2][0];
# compute range
$range = ($self->[2][-1] - $self->[2][0]);
} else {
# set offset
$off = $self->[2][-1];
# compute range
$range = ($self->[2][0] - $self->[2][-1]);
}
# verify range
($range != 0) || croak('cannot normalize output values');
lib/ICC/Support/spline.pm view on Meta::CPAN
} else {
# return derivative
return(derivative($self, $in));
}
}
# combined transform
# direction: 0 - normal, 1 - inverse
# parameters: (object_reference, direction, input_value)
# returns: (output_value)
sub _transform {
# get parameters
my ($self, $dir, $in) = @_;
# if inverse direction
if ($dir) {
# return inverse
return(scalar(inverse($self, $in)));
} else {
# return transform
return(transform($self, $in));
}
}
# compute knot derivatives for natural spline
# parameters: (ref_to_object)
sub _objderv {
# get parameter
my $self = shift();
# local variables
my ($ix, $rhs, $info, $derv);
# verify object has two or more knots
(($ix = $#{$self->[2]}) > 0) || croak('object must have two or more knots');
# check if ICC::Support::Lapack module is loaded
state $lapack = defined($INC{'ICC/Support/Lapack.pm'});
# make empty Math::Matrix object
$rhs = bless([], 'Math::Matrix');
# for each input element
for my $i (1 .. $ix - 1) {
# compute rhs (3 * (y[i + 1] - y[i - 1]))
$rhs->[$i][0] = 3 * ($self->[2][$i + 1] - $self->[2][$i - 1]);
}
# set rhs endpoint values
$rhs->[0][0] = 3 * ($self->[2][1] - $self->[2][0]);
$rhs->[$ix][0] = 3 * ($self->[2][$ix] - $self->[2][$ix - 1]);
# if ICC::Support::Lapack module is loaded
if ($lapack) {
# solve for derivative matrix
($info, $derv) = ICC::Support::Lapack::trisolve([(1) x $ix], [2, (4) x ($ix - 1), 2], [(1) x $ix], $rhs);
# otherwise, use Math::Matrix module
} else {
# solve for derivative matrix
$derv = Math::Matrix->tridiagonal([2, (4) x ($ix - 1), 2])->concat($rhs)->solve();
}
# for each knot
for my $i (0 .. $ix) {
# set derivative value
$self->[3][$i] = $derv->[$i][0];
}
}
# compute partial derivative matrix for natural spline
# parameters: (ref_to_object)
sub _objpara {
# get parameter
my $self = shift();
# local variables
my ($ix, $rhs, $info, $derv);
# verify object has two or more knots
(($ix = $#{$self->[2]}) > 0) || croak('object must have two or more knots');
# check if ICC::Support::Lapack module is loaded
state $lapack = defined($INC{'ICC/Support/Lapack.pm'});
# make rhs matrix (filled with zeros)
$rhs = $lapack ? ICC::Support::Lapack::zeros($ix + 1) : bless([map {[(0) x ($ix + 1)]} (0 .. $ix)], 'Math::Matrix');
# for each row
for my $i (1 .. $ix - 1) {
# set rhs diagonal values
$rhs->[$i - 1][$i] = 3;
$rhs->[$i + 1][$i] = -3;
}
# set rhs endpoint values
$rhs->[0][0] = -3;
$rhs->[1][0] = -3;
$rhs->[$ix][$ix] = 3;
$rhs->[$ix - 1][$ix] = 3;
# if ICC::Support::Lapack module is loaded
if ($lapack) {
# solve for derivative matrix
($info, $derv) = ICC::Support::Lapack::trisolve([(1) x $ix], [2, (4) x ($ix - 1), 2], [(1) x $ix], $rhs);
# set object
$self->[5] = bless($derv, 'Math::Matrix');
# otherwise, use Math::Matrix module
} else {
# solve for derivative matrix
$self->[5] = Math::Matrix->tridiagonal([2, (4) x ($ix - 1), 2])->concat($rhs)->solve();
}
}
# add local min/max values
# parameters: (ref_to_object)
# returns: (s-value_array)
sub _minmax {
# get object reference
my $self = shift();
# local variables
my (@t, @y, @s);
# for each interval
for my $i (0 .. $#{$self->[2]} - 1) {
# get min/max location(s)
@t = _local($self, $i);
# compute local min/max y-values
@y = map {_fwd($self, $_, $i)} @t;
# add the knot values
push(@y, ($self->[2][$i], $self->[2][$i + 1]));
# sort the values
@y = sort {$a <=> $b} @y;
# save y-value min/max span
$self->[4][$i] = [$y[0], $y[-1]];
# add min/max s-values
push(@s, map {$_ + $i} @t);
}
# return min/max s-values
return(sort {$a <=> $b} @s);
}
# compute local min/max value(s)
# values are within the segment (t)
# parameters: (ref_to_object, segment)
# returns: (t-value_array)
sub _local {
# get parameters
my ($self, $low) = @_;
# local variables
my ($y1, $y2, $m1, $m2);
my ($a, $b, $c, $dscr, @t);
# get endpoint values
$y1 = $self->[2][$low];
$y2 = $self->[2][$low + 1];
$m1 = $self->[3][$low];
$m2 = $self->[3][$low + 1];
# compute coefficients of quadratic equation (at^2 + bt + c = 0)
$a = 6 * ($y1 - $y2) + 3 * ($m1 + $m2);
$b = -6 * ($y1 - $y2) - 2 * $m2 - 4 * $m1;
$c = $m1;
# return if constant
return() if (abs($a) < 1E-15 && abs($b) < 1E-15);
# if linear equation (a is zero)
if (abs($a) < 1E-15) {
# add solution
push(@t, -$c/$b);
# if quadratic equation
} else {
# compute discriminant
$dscr = $b**2 - 4 * $a * $c;
# if discriminant > zero
if ($dscr > 0) {
# add solutions (critical points)
push(@t, -($b + sqrt($dscr))/(2 * $a));
push(@t, -($b - sqrt($dscr))/(2 * $a));
}
}
# return solution(s) within interval (0 < t < 0)
return (grep {$_ > 0 && $_ < 1} @t);
}
# compute second derivative for unit spline segment
# parameters: (ref_to_object, t-value, segment)
# returns: (second_derivative)
sub _derv2 {
# get parameters
my ($self, $t, $low) = @_;
# local variables
my ($h00, $h01, $h10, $h11);
# compute Hermite derivative coefficients
$h00 = 12 * $t - 6;
$h01 = - $h00;
$h10 = 6 * $t - 4;
$h11 = 6 * $t - 2;
# interpolate value
return($h00 * $self->[2][$low] + $h01 * $self->[2][$low + 1] + $h10 * $self->[3][$low] + $h11 * $self->[3][$low + 1]);
lib/ICC/Support/spline.pm view on Meta::CPAN
# compute transform value
# parameters: (ref_to_object, t-value, segment)
# returns: (y-value)
sub _fwd {
# get parameters
my ($self, $t, $low) = @_;
# local variables
my ($tc, $h00, $h01, $h10, $h11);
# check if ICC::Support::Lapack module is loaded
state $lapack = defined($INC{'ICC/Support/Lapack.pm'});
# if position is non-zero
if ($t) {
# if ICC::Support::Lapack module is loaded
if ($lapack) {
# interpolate value (Lapack XS)
return(ICC::Support::Lapack::hermite($t, $self->[2][$low], $self->[2][$low + 1], $self->[3][$low], $self->[3][$low + 1]));
} else {
# compute Hermite coefficients
$tc = 1 - $t;
$h00 = (1 + 2 * $t) * $tc * $tc;
$h01 = 1 - $h00;
$h10 = $t * $tc * $tc;
$h11 = -$t * $t * $tc;
# interpolate value (no Lapack XS)
return($h00 * $self->[2][$low] + $h01 * $self->[2][$low + 1] + $h10 * $self->[3][$low] + $h11 * $self->[3][$low + 1]);
}
} else {
# use lower knot output value
return($self->[2][$low]);
}
}
# compute inverse value
# parameters: (ref_to_object, y-value, segment)
# there are 0 to 3 solutions in the range 0 < t < 1
# returns: (t-value_array)
sub _rev {
# get parameters
my ($self, $in, $low) = @_;
# local variables
my ($y1, $y2, $m1, $m2);
my ($a, $b, $c, $d, $dscr, @t);
my ($d0, $d1, $cs, $cc, @r, $ccr, $sol, $lim0, $lim1);
# get endpoint values
$y1 = $self->[2][$low];
$y2 = $self->[2][$low + 1];
$m1 = $self->[3][$low];
$m2 = $self->[3][$low + 1];
# compute coefficients of cubic equation (at^3 + bt^2 + ct + d = 0)
$a = 2 * ($y1 - $y2) + $m1 + $m2;
$b = -3 * ($y1 - $y2) -2 * $m1 - $m2;
$c = $m1;
$d = $y1 - $in;
# constant equation (a, b, c are zero)
if (abs($a) < 5E-15 && abs($b) < 5E-15 && abs($c) < 5E-15) {
# add solution
push(@t, 0.5) if ($d == 0);
# linear equation (a, b are zero)
} elsif (abs($a) < 5E-15 && abs($b) < 5E-15) {
# add solution
push(@t, -$d/$c);
# quadratic equation (a is zero)
} elsif (abs($a) < 5E-15) {
# compute discriminant: > 0 two real, == 0 one real, < 0 two complex
$dscr = $c**2 - 4 * $b * $d;
# if discriminant is zero
if ($dscr == 0) {
# add solution (double root)
push(@t, -$c/(2 * $b));
# if discriminant > zero
} elsif ($dscr > 0) {
# add solutions
push(@t, -($c + sqrt($dscr))/(2 * $b));
push(@t, -($c - sqrt($dscr))/(2 * $b));
}
# cubic equation
} else {
# compute discriminant: > 0 three real, == 0 one or two real, < 0 one real and two complex
$dscr = 18 * $a * $b * $c * $d - 4 * $b**3 * $d + $b**2 * $c**2 - 4 * $a * $c**3 - 27 * $a**2 * $d**2;
# compute ∆0
$d0 = $b**2 - 3 * $a * $c;
# if discriminant is zero
if ($dscr == 0) {
# if ∆0 is zero
if ($d0 == 0) {
# add solution (triple root)
( run in 0.667 second using v1.01-cache-2.11-cpan-524268b4103 )