package ICC::Support::spline;
use strict;
use Carp;
our $VERSION = 0.12;
# revised 2018-09-04
#
# Copyright © 2004-2020 by William B. Birkett
# inherit from Shared
use parent qw(ICC::Shared);
# add complex math functions
use Math::Complex qw(root sqrt cbrt Im Re);
# enable static variables
use feature 'state';
# create new spline object
# hash keys are: 'input', 'output', 'type', 'damp', 'fix_hl', 'fix_sh'
# default values for 'input' and 'output' are [0, 1]
# parameters: ([ref_to_attribute_hash])
# returns: (ref_to_object)
sub new {
# get object class
my $class = shift();
# create empty spline object
my $self = [
{}, # object header
[], # input values
[], # output values
[], # derivative values
[], # min/max output values
[], # parametric derivative matrix
];
# if one parameter, a hash reference
if (@_ == 1 && ref($_[0]) eq 'HASH') {
# set object contents from hash
_new_from_hash($self, $_[0]);
} elsif (@_) {
# error
croak('invalid parameter(s)');
}
# return blessed object
return(bless($self, $class));
}
# create inverse object
# swaps the input and output arrays
# returns: (ref_to_object)
sub inv {
# clone new object
my $inv = Storable::dclone(shift());
# get min/max values
my @s = _minmax($inv);
# verify object is monotonic
(! @s) or croak("cannot invert a non-monotonic 'spline' object");
# if input is a range
if (@{$inv->[1]} == 2 && @{$inv->[2]} > 2) {
# interpolate input array
$inv->[1] = [map {(1 - $_/$#{$inv->[2]}) * $inv->[1][0] + $_/$#{$inv->[2]} * $inv->[1][-1]} (0 .. $#{$inv->[2]})];
}
# swap input and output arrays
($inv->[1], $inv->[2]) = ($inv->[2], $inv->[1]);
# clear parameteric derivative matrix
$inv->[5] = [];
# sort input/output values
_sort_values($inv);
# compute knot derivatives
_objderv($inv);
# add segment min/max y-values
_minmax($inv);
# return inverted object
return($inv);
}
# write spline object to ICC profile
# note: writes an equivalent 'curv' object
# parameters: (ref_to_parent_object, file_handle, ref_to_tag_table_entry)
sub write_fh {
# verify 4 parameters
(@_ == 4) or croak('wrong number of parameters');
# write spline data to profile
goto &_writeICCspline;
}
# get tag size (for writing to profile)
# note: writes an equivalent curv object
# returns: (tag_size)
sub size {
# get parameters
my ($self) = @_;
# get count value (defaults to 4096, the maximum allowed in an 'mft2' tag)
my $n = $self->[0]{'curv_points'} // 4096;
# return size
return(12 + $n * 2);
}
# compute curve function
# parameters: (input_value)
# returns: (output_value)
sub transform {
# get parameters
my ($self, $in) = @_;
# local variable
my ($low);
# if an extrapolated solution (s < 0)
if ($in < $self->[1][0]) {
# return extrapolated value
return($self->[2][0] + $self->[3][0] * ($in - $self->[1][0]));
# if an extrapolated solution (s >= 1)
} elsif ($in >= $self->[1][-1]) {
# return extrapolated value
return($self->[2][-1] + $self->[3][-1] * ($in - $self->[1][-1]));
# if input array contains x-values
} elsif ($#{$self->[1]} > 1) {
# initilize segment
$low = 0;
# while input value > upper knot value
while ($in > $self->[1][$low + 1]) {
# increment segment
$low++;
}
# return interpolated value
return(_fwd($self, ($self->[1][$low] - $in)/($self->[1][$low] - $self->[1][$low + 1]), $low));
# input array contains range
} else {
# return interpolated value
return(_fwd($self, POSIX::modf($#{$self->[2]} * ($in - $self->[1][0])/($self->[1][1] - $self->[1][0]))));
}
}
# compute inverse curve function
# note: there may be multiple solutions
# parameters: (input_value)
# returns: (output_value -or- array_of_output_values)
sub inverse {
# get parameters
my ($self, $in) = @_;
# local variables
my ($ix, @sol, @t);
# get output array upper index
$ix = $#{$self->[2]};
# if an extrapolated solution (s < 0)
if (($in < $self->[2][0] && $self->[3][0] > 0) || ($in > $self->[2][0] && $self->[3][0] < 0)) {
# add solution
push(@sol, $self->[1][0] + ($in - $self->[2][0])/$self->[3][0]);
}
# for each segment (0 <= s < 1)
for my $i (0 .. $ix - 1) {
# if input is lower knot
if ($in == $self->[2][$i]) {
# add solution (x-values : range)
push(@sol, ($#{$self->[1]} > 1) ? $self->[1][$i] : ($self->[1][0] * ($ix - $i) + $self->[1][1] * $i)/$ix);
}
# 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);
# if input array contains x-values
if ($#{$self->[1]} > 1) {
# add solution(s)
push(@sol, map {(1 - $_) * $self->[1][$i] + $_ * $self->[1][$i + 1]} @t);
# input array contains range
} else {
# add solution(s)
push(@sol, map {($self->[1][0] * ($ix - $i - $_) + $self->[1][1] * ($i + $_))/$ix} @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] + ($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) = @_;
# local variable
my ($low);
# if an extrapolated solution (s < 0)
if ($in < $self->[1][0]) {
# return endpoint derivative
return($self->[3][0]);
# if an extrapolated solution (s >= 1)
} elsif ($in >= $self->[1][-1]) {
# return endpoint derivative
return($self->[3][-1]);
# if input array contains x-values
} elsif ($#{$self->[1]} > 1) {
# initilize segment
$low = 0;
# while input value > upper knot value
while ($in > $self->[1][$low + 1]) {
# increment segment
$low++;
}
# return interpolated derivative value
return(_derv($self, ($self->[1][$low] - $in)/($self->[1][$low] - $self->[1][$low + 1]), $low));
# input array contains range
} else {
# return interpolated derivative value
return(_derv($self, POSIX::modf($#{$self->[2]} * ($in - $self->[1][0])/($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 ($dx, $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 ($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] * ($in - $self->[1][0]);
}
# if an extrapolated solution (s >= 1)
} elsif ($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] * ($in - $self->[1][-1]);
}
} else {
# if input array contains x-values
if ($#{$self->[1]} > 1) {
# initialize segment
$low = 0;
# while input value > upper knot value
while ($in > $self->[1][$low + 1]) {
# increment segment
$low++;
}
# compute delta x-value
$dx = $self->[1][$low + 1] - $self->[1][$low];
# compute t-value
$t = ($in - $self->[1][$low])/$dx;
# input array contains range
} else {
# compute delta x-value
$dx = ($self->[1][1] - $self->[1][0])/$#{$self->[2]};
# compute t-value, index
($t, $low) = POSIX::modf(($in - $self->[1][0])/$dx);
}
# 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)
for my $i (0 .. $#{$self->[2]}) {
# compute partial derivative
$pj[$i] = $h00 * ($low == $i) + $h01 * ($low + 1 == $i) + $h10 * $dx * $self->[5][$low][$i] + $h11 * $dx * $self->[5][$low + 1][$i];
}
}
# return
return(@pj);
}
# get/set header hash reference
# parameters: ([ref_to_new_hash])
# returns: (ref_to_hash)
sub header {
# get object reference
my $self = shift();
# if there are parameters
if (@_) {
# if one parameter, a hash reference
if (@_ == 1 && ref($_[0]) eq 'HASH') {
# set header to copy of hash
$self->[0] = Storable::dclone(shift());
} else {
# error
croak('parameter must be a hash reference');
}
}
# return reference
return($self->[0]);
}
# get/set input array reference
# an array of two values is a range
# otherwise, it contains input values
# parameters: ([ref_to_array])
# returns: (ref_to_array)
sub input {
# get object reference
my $self = shift();
# if one parameter supplied
if (@_ == 1) {
# verify 'input' vector (two or more numeric values)
(ICC::Shared::is_num_vector($_[0]) && 2 <= @{$_[0]}) or croak('invalid input values');
# verify 'input' vector size
(@{$_[0]} == 2 || @{$_[0]} == @{$self->[2]}) or carp('number of input values differs from object output values');
# store output values
$self->[1] = [@{$_[0]}];
# update object internal elements
update($self, 1);
} elsif (@_) {
# error
carp('too many parameters');
}
# return array reference
return($self->[1]);
}
# get/set output array reference
# parameters: ([ref_to_array])
# returns: (ref_to_array)
sub output {
# get object reference
my $self = shift();
# if one parameter supplied
if (@_ == 1) {
# verify 'output' vector (two or more numeric values)
(ICC::Shared::is_num_vector($_[0]) && 2 <= @{$_[0]}) or croak('invalid output values');
# verify 'output' vector size
(@{$self->[1]} == 2 || @{$_[0]} == @{$self->[1]}) or carp('number of output values differs from object input values');
# store output values
$self->[2] = [@{$_[0]}];
# update object internal elements
update($self);
} elsif (@_) {
# error
carp('too many parameters');
}
# return array reference
return($self->[2]);
}
# normalize transform
# adjusts Bernstein coefficients linearly
# adjusts 'input' and 'output' by default
# hash keys: 'input', 'output'
# parameter: ([hash_ref])
# returns: (ref_to_object)
sub normalize {
# get parameters
my ($self, $hash) = @_;
# local variables
my ($m, $b, $val, $src, $x0, $x1, $y0, $y1, $f0, $f1, @minmax, $p);
# set hash key array
my @key = qw(input output);
# for ('input', 'output')
for my $i (0, 1) {
# undefine slope
$m = undef;
# if no parameter (default)
if (! defined($hash)) {
# if array has 2 or more coefficients
if (@{$self->[$i + 1]} > 1 && ($x0 = $self->[$i + 1][0]) != ($x1 = $self->[$i + 1][-1])) {
# special case, if 'output' and 'fix_hl' was called
$x0 = 0 if ($i && defined($self->[0]{'hl_retention'}));
# compute adjustment values
$m = abs(1/($x0 - $x1));
$b = -$m * ($x0 < $x1 ? $x0 : $x1);
# set endpoints
($y0, $y1) = $x0 < $x1 ? (0, 1) : (1, 0);
}
} else {
# if hash contains key
if (defined($val = $hash->{$key[$i]})) {
# if array has 2 or more coefficients
if (@{$self->[$i + 1]} > 1 && $self->[$i + 1][0] != $self->[$i + 1][-1]) {
# if hash value is an array reference
if (ref($val) eq 'ARRAY') {
# if two parameters
if (@{$val} == 2) {
# get parameters
($y0, $y1) = @{$val};
# get endpoints
$x0 = $self->[$i + 1][0];
$x1 = $self->[$i + 1][-1];
# if three parameters
} elsif (@{$val} == 3) {
# get parameters
($src, $y0, $y1) = @{$val};
# if source is 'endpoints'
if (! ref($src) && $src eq 'endpoints') {
# get endpoints
$x0 = $self->[$i + 1][0];
$x1 = $self->[$i + 1][-1];
# if source is 'minmax'
} elsif (! ref($src) && $src eq 'minmax') {
# if 'output' curve
if ($i) {
# get all min/max values
@minmax = (map {@{$_}} @{$self->[4]});
# get min/max values
$x0 = List::Util::min(@minmax);
$x1 = List::Util::max(@minmax);
# if 'input' curve
} else {
# get min/max values (endpoints)
$x0 = $self->[1][0];
$x1 = $self->[1][-1];
}
}
# if four parameters
} elsif (@{$val} == 4) {
# get parameters
($x0, $x1, $y0, $y1) = @{$val};
}
# if x/y values are numeric
if ((4 == grep {Scalar::Util::looks_like_number($_)} ($x0, $x1, $y0, $y1)) && ($x0 != $x1)) {
# compute slope
$m = ($y0 - $y1)/($x0 - $x1);
# compute offset
$b = $y0 - $m * $x0;
} else {
# warning
carp("invalid '$key[$i]' parameter(s)");
}
}
} else {
# warning
carp("can't normalize '$key[$i]' coefficients");
}
}
}
# if slope is defined
if (defined($m)) {
# set endpoint flags
$f0 = $self->[$i + 1][0] == $x0;
$f1 = $self->[$i + 1][-1] == $x1;
# adjust Bernstein coefficients (y = mx + b)
@{$self->[$i + 1]} = map {$m * $_ + $b} @{$self->[$i + 1]};
# set endpoints (to ensure exact values)
$self->[$i + 1][0] = $y0 if ($f0);
$self->[$i + 1][-1] = $y1 if ($f1);
# save input polarity
$p = $m < 0 if ($i == 0);
# if 'input' and 'hl_retention' is defined
if (! $i && defined($val = $self->[0]{'hl_retention'})) {
# adjust value
$self->[0]{'hl_retention'} = $m * $val + $b;
}
# if 'output' and 'hl_original' is defined
if ($i && defined($val = $self->[0]{'hl_original'})) {
# adjust values
$self->[0]{'hl_original'} = [map {$m * $_ + $b} @{$val}];
}
# if 'output' and 'sh_original' is defined
if ($i && defined($val = $self->[0]{'sh_original'})) {
# adjust values
$self->[0]{'sh_original'} = [map {$m * $_ + $b} @{$val}];
}
}
}
# update object internals, sorting values if 'input' polarity changed
update($self, $p);
# return
return($self);
}
# check if monotonic
# returns min/max values
# returns s-values by default
# returns input-values if format is 'input'
# returns output-values if format is 'output'
# curve is monotonic if no min/max values
# parameters: ([format])
# returns: (array_of_values)
sub monotonic {
# get object reference
my ($self, $fmt) = @_;
# local variables
my ($t, $low);
# if format undefined
if (! defined($fmt)) {
# return s-values
return(_minmax($self));
# if format 'input'
} elsif ($fmt eq 'input') {
# if input array contains x-values
if ($#{$self->[1]} > 1) {
# return input values
return(map {($t, $low) = POSIX::modf($_); (1 - $t) * $self->[1][$low] + $t * $self->[1][$low + 1]} _minmax($self));
# input array contains range
} else {
# return input values
return(map {$t = $_/$#{$self->[2]}; (1 - $t) * $self->[1][0] + $t * $self->[1][1]} _minmax($self));
}
# if format 'output'
} elsif ($fmt eq 'output') {
# return output values
return(map {_fwd($self, POSIX::modf($_))} _minmax($self));
} else {
# error
croak('unsupported format for min/max values');
}
}
# update internal object elements
# this method used primarily when optimizing
# call after changing 'input' or 'output' values
# if the 'input' values were changed, set the flag
# parameter: ([flag])
sub update {
# get parameters
my ($self, $flag) = @_;
# sort values if flag set
_sort_values($self) if ($flag);
# recompute knot derivatives
_objderv($self);
# recompute segment min/max y-values
_minmax($self);
# recompute parametric derivative matrix if flag set, or 'akima' type spline
_objpara($self) if ($flag || (defined($self->[0]{'type'} && $self->[0]{'type'} eq 'akima')));
}
# 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) or 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 equivalent '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(@_)));
}
# print object contents to string
# format is an array structure
# parameter: ([format])
# returns: (string)
sub sdump {
# get parameters
my ($self, $p) = @_;
# local variables
my ($s, $fmt, $fmt2);
# resolve parameter to an array reference
$p = defined($p) ? ref($p) eq 'ARRAY' ? $p : [$p] : [];
# get format string
$fmt = defined($p->[0]) && ! ref($p->[0]) ? $p->[0] : 'undef';
# set string to object ID
$s = sprintf("'%s' object, (0x%x)\n", ref($self), $self);
# if input_parameter_array defined
if (defined($self->[1])) {
# set format
$fmt2 = join(', ', ('%.3f') x @{$self->[1]});
# append parameter string
$s .= sprintf(" input parameters [$fmt2]\n", @{$self->[1]});
}
# if output_parameter_array defined
if (defined($self->[2])) {
# set format
$fmt2 = join(', ', ('%.3f') x @{$self->[2]});
# append parameter string
$s .= sprintf(" output parameters [$fmt2]\n", @{$self->[2]});
}
# return
return($s);
}
# combined transform
# direction: 0 - normal, 1 - inverse
# parameters: (ref_to_object, 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));
}
}
# combined derivative
# direction: 0 - normal, 1 - inverse
# parameters: (ref_to_object, direction, input_value)
# returns: (derivative_value)
sub _derivative {
# get parameters
my ($self, $dir, $in) = @_;
# if inverse direction
if ($dir) {
# unimplemented function error
croak('unimplemented function');
} else {
# return derivative
return(derivative($self, $in));
}
}
# compute knot derivatives
# parameter: (ref_to_object)
sub _objderv {
# get parameter
my $self = shift();
# verify input values are unique
_unique($self->[1]) or croak('input values not unique');
# if spline type undefined or 'natural'
if (! defined($self->[0]{'type'}) || $self->[0]{'type'} eq 'natural') {
# compute knot derivatives for natural spline
_natural_derv($self)
# if spline type is 'akima'
} elsif ($self->[0]{'type'} eq 'akima') {
# compute knot derivatives for Akima spline
_akima_derv($self);
} else {
# error
croak('undefined spline type');
}
}
# compute knot derivatives for natural spline
# parameter: (ref_to_object)
sub _natural_derv {
# get parameter
my $self = shift();
# local variables
my ($ix, $rhs, $dx, $info, $derv);
# verify object has two or more knots
(($ix = $#{$self->[2]}) > 0) or 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');
# if input array contains input values
if ($#{$self->[1]} > 1) {
# for each inner element
for my $i (1 .. $ix - 1) {
# compute rhs (6 * (y[i + 1] - y[i - 1])/(x[i + 1] - x[i - 1]))
$rhs->[$i][0] = 6 * ($self->[2][$i + 1] - $self->[2][$i - 1])/($self->[1][$i + 1] - $self->[1][$i - 1]);
}
# set rhs endpoint values
$rhs->[0][0] = 3 * ($self->[2][1] - $self->[2][0])/($self->[1][1] - $self->[1][0]);
$rhs->[$ix][0] = 3 * ($self->[2][$ix] - $self->[2][$ix - 1])/($self->[1][$ix] - $self->[1][$ix - 1]);
# input array contains range
} else {
# compute x-segment length
$dx = ($self->[1][1] - $self->[1][0])/$ix;
# for each inner element
for my $i (1 .. $ix - 1) {
# compute rhs (6 * (y[i + 1] - y[i - 1]))/(x[i + 1] - x[i - 1])
$rhs->[$i][0] = 3 * ($self->[2][$i + 1] - $self->[2][$i - 1])/$dx;
}
# set rhs endpoint values
$rhs->[0][0] = 3 * ($self->[2][1] - $self->[2][0])/$dx;
$rhs->[$ix][0] = 3 * ($self->[2][$ix] - $self->[2][$ix - 1])/$dx;
}
# 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 knot derivatives for Akima spline
# parameter: (ref_to_object)
sub _akima_derv {
# get parameter
my $self = shift();
# local variables
my ($ix, @m, $s, $b, $d1, $d2);
# verify object has two or more knots
(($ix = $#{$self->[2]}) > 0) or croak('object must have two or more knots');
# compute uniform segment length (for 'input' range)
$s = ($self->[1][1] - $self->[1][0])/$ix if ($#{$self->[1]} == 1);
# for each segment
for my $i (0 .. ($ix - 1)) {
# compute segment slope (range // x-values)
$m[$i + 2] = ($self->[2][$i + 1] - $self->[2][$i])/($s // ($self->[1][$i + 1] - $self->[1][$i]));
}
# if 2 points
if ($ix == 1) {
# linear slopes
$m[0] = $m[1] = $m[3] = $m[4] = $m[2];
# if 3 or more points
} else {
# quadratic slopes
$m[1] = 2 * $m[2] - $m[3];
$m[0] = 2 * $m[1] - $m[2];
$m[$ix + 2] = 2 * $m[$ix + 1] - $m[$ix];
$m[$ix + 3] = 2 * $m[$ix + 2] - $m[$ix + 1];
}
# if damping value is undefined
if (! defined($self->[0]{'damp'})) {
# set default damping value
$self->[0]{'damp'} = List::Util::sum(map {$_ * $_} @m)/(@m * 50);
}
# Akima used abs() functions to calulate the slopes
# this can lead to abrupt changes in curve shape and parametric derivatives
# we replace abs(x) with sqrt(x^2 + b) to add stability
# get damping value and verify >= 0
(($b = $self->[0]{'damp'}) >= 0) or carp('akima damping value is negative');
# compute Akima derivative
for my $i (0 .. $ix) {
# compute slope differences (with damping)
$d1 = sqrt(($m[$i + 1] - $m[$i])**2 + $b);
$d2 = sqrt(($m[$i + 3] - $m[$i + 2])**2 + $b);
# if denominator not 0
if ($d1 + $d2) {
# use Akima slope (weighted average with damping)
$self->[3][$i] = ($d2 * $m[$i + 1] + $d1 * $m[$i + 2])/($d1 + $d2);
} else {
# otherwise, use average slope of adjoining segments
$self->[3][$i] = ($m[$i + 1] + $m[$i + 2])/2;
}
}
}
# compute partial derivative matrix
# parameter: (ref_to_object)
sub _objpara {
# get parameter
my $self = shift();
# verify input values are unique
_unique($self->[1]) or croak('input values not unique');
# if spline type undefined or 'natural'
if (! defined($self->[0]{'type'}) || $self->[0]{'type'} eq 'natural') {
# compute knot derivatives for natural spline
_natural_para($self)
# if spline type is 'akima'
} elsif ($self->[0]{'type'} eq 'akima') {
# compute knot derivatives for Akima spline
_akima_para($self);
} else {
# error
croak('undefined spline type');
}
}
# compute partial derivative matrix for natural spline
# parameter: (ref_to_object)
sub _natural_para {
# get parameter
my $self = shift();
# local variables
my ($ix, $n, $rhs, $x, $info, $derv);
# verify object has two or more knots
(($ix = $#{$self->[2]}) > 0) or 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');
# if input array contains input values
if ($#{$self->[1]} > 1) {
# for each inner element
for my $i (1 .. $ix - 1) {
# set rhs diagonal values
$rhs->[$i][$i - 1] = $x = 6/($self->[1][$i - 1] - $self->[1][$i + 1]);
$rhs->[$i][$i + 1] = -$x;
}
# set rhs endpoint values
$rhs->[0][0] = $x = 3/($self->[1][0] - $self->[1][1]);
$rhs->[0][1] = -$x;
$rhs->[$ix][$ix] = $x = 3/($self->[1][$ix] - $self->[1][$ix - 1]);
$rhs->[$ix][$ix -1] = -$x;
# input array contains range
} else {
# compute rhs value
$x = 3 * $ix/($self->[1][1] - $self->[1][0]);
# for each inner element
for my $i (1 .. $ix - 1) {
# set rhs diagonal values
$rhs->[$i - 1][$i] = $x;
$rhs->[$i + 1][$i] = -$x;
}
# set rhs endpoint values
$rhs->[0][0] = -$x;
$rhs->[1][0] = -$x;
$rhs->[$ix][$ix] = $x;
$rhs->[$ix - 1][$ix] = $x;
}
# 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();
}
}
# compute partial derivative matrix for Akima spline
# parameter: (ref_to_object)
sub _akima_para {
# get object
my $self = shift();
# local variables
my ($ix, $s, @m, @x, $b);
my ($d1, $d2, $dif, @dm, $dd1, $dd2, $derv);
# see 'akima_parametric.plx' for details on this function
# verify object has two or more knots
(($ix = $#{$self->[2]}) > 0) or croak('object must have two or more knots');
# compute uniform segment length (for 'input' range)
$s = ($self->[1][1] - $self->[1][0])/$ix if ($#{$self->[1]} == 1);
# for each segment
for my $i (0 .. ($ix - 1)) {
# compute segment x-length
$x[$i + 2] = $s // ($self->[1][$i + 1] - $self->[1][$i]);
# compute segment slope (range // x-values)
$m[$i + 2] = ($self->[2][$i + 1] - $self->[2][$i])/$x[$i + 2];
}
# if 2 points
if ($ix == 1) {
# linear slopes
$m[0] = $m[1] = $m[3] = $m[4] = $m[2];
# if 3 or more points
} else {
# quadratic slopes
$m[1] = 2 * $m[2] - $m[3];
$m[0] = 2 * $m[1] - $m[2];
$m[$ix + 2] = 2 * $m[$ix + 1] - $m[$ix];
$m[$ix + 3] = 2 * $m[$ix + 2] - $m[$ix + 1];
}
# get the damping value
$b = $self->[0]{'damp'};
# for each node
for my $i (0 .. $ix) {
# add row of zeros
$derv->[$i] = [(0) x ($ix + 1)];
}
# if 2 knots
if ($ix == 1) {
# set partial derivatives
$derv = [[-1/$x[2], 1/$x[2]], [-1/$x[2], 1/$x[2]]];
} else {
# for each row
for my $i (0 .. $ix) {
# compute slope differences (with damping)
$d1 = sqrt(($m[$i + 1] - $m[$i])**2 + $b);
$d2 = sqrt(($m[$i + 3] - $m[$i + 2])**2 + $b);
# for each column
for my $j (0 .. $ix) {
# compute difference
$dif = $j - $i;
# skip outer corner regions
next if (abs($dif) > 2);
# the @dm array contains (∂m0/∂n, ∂m1/∂n, ∂m2/∂n, ∂m3/∂n)
# where m0 - m3 are the chord slopes, and n is the node value
# if difference -2
if ($dif == -2) {
# if row 1
if ($i == $ix) {
# set partial derivatives
@dm = (-1/$x[$i], 0, 1/$x[$i], 2/$x[$i]);
} else {
# set partial derivatives
@dm = (-1/$x[$i], 0, 0, 0);
}
# if difference -1
} elsif ($dif == -1) {
# if row 1
if ($i == 1) {
# if more than 3 knots
if ($ix > 2) {
# set partial derivatives
@dm = (-2/$x[$i + 1], -1/$x[$i + 1], 0, 0);
} else {
# set partial derivatives
@dm = (-2/$x[$i + 1], -1/$x[$i + 1], 0, 1/$x[$i + 1]);
}
# if next to last row
} elsif ($i == $ix - 1) {
# set partial derivatives
@dm = (1/$x[$i], -1/$x[$i + 1], 0, 1/$x[$i + 1]);
# if last row
} elsif ($i == $ix) {
# set partial derivatives
@dm = (1/$x[$i], -1/$x[$i + 1], -2/$x[$i + 1] - 1/$x[$i], -3/$x[$i + 1] - 2/$x[$i]);
} else {
# set partial derivatives
@dm = (1/$x[$i], -1/$x[$i + 1], 0, 0);
}
# if difference 0
} elsif ($dif == 0) {
# if row 0
if ($i == 0) {
# set partial derivatives
@dm = (-3/$x[$i + 2], -2/$x[$i + 2], -1/$x[$i + 2], 0);
# if row 1
} elsif ($i == 1) {
# if more than three knots
if ($ix > 2) {
# set partial derivatives
@dm = (2/$x[$i + 1] + 1/$x[$i + 2], 1/$x[$i + 1], -1/$x[$i + 2], 0);
} else {
# set partial derivatives
@dm = (2/$x[$i + 1] + 1/$x[$i + 2], 1/$x[$i + 1], -1/$x[$i + 2], -2/$x[$i + 2] - 1/$x[$i + 1]);
}
# if next to last row
} elsif ($i == $ix - 1) {
# set partial derivatives
@dm = (0, 1/$x[$i + 1], -1/$x[$i + 2], -2/$x[$i + 2] - 1/$x[$i + 1]);
# if last row
} elsif ($i == $ix) {
# set partial derivatives
@dm = (0, 1/$x[$i + 1], 2/$x[$i + 1], 3/$x[$i + 1]);
} else {
# set partial derivatives
@dm = (0, 1/$x[$i + 1], -1/$x[$i + 2], 0);
}
# if difference 1
} elsif ($dif == 1) {
# if row 0
if ($i == 0) {
# set partial derivatives
@dm = (3/$x[$i + 2] + 2/$x[$i + 3], 2/$x[$i + 2] + 1/$x[$i + 3], 1/$x[$i + 2], -1/$x[$i + 3]);
# if row 1
} elsif ($i == 1) {
# if more than 3 knots
if ($ix > 2) {
# set partial derivatives
@dm = (-1/$x[$i + 2], 0, 1/$x[$i + 2], -1/$x[$i + 3]);
} else {
# set partial derivatives
@dm = (-1/$x[$i + 2], 0, 1/$x[$i + 2], 2/$x[$i + 2]);
}
# if next to last row
} elsif ($i == $ix - 1) {
# set partial derivatives
@dm = (0, 0, 1/$x[$i + 2], 2/$x[$i + 2]);
} else {
# set partial derivatives
@dm = (0, 0, 1/$x[$i + 2], -1/$x[$i + 3]);
}
# if difference 2
} elsif ($dif == 2) {
# if row 0
if ($i == 0) {
# set partial derivatives
@dm = (-2/$x[$i + 3], -1/$x[$i + 3], 0, 1/$x[$i + 3]);
} else {
# set partial derivatives
@dm = (0, 0, 0, 1/$x[$i + 3]);
}
}
# compute ∂d1/∂n
$dd1 = akima_dpd($m[$i], $m[$i + 1], $dm[0], $dm[1], $d1);
# compute ∂d2/∂n
$dd2 = akima_dpd($m[$i + 2], $m[$i + 3], $dm[2], $dm[3], $d2);
# compute ∂s/∂n
$derv->[$i][$j] = akima_spd($d1, $d2, $m[$i + 1], $m[$i + 2], $dd1, $dd2, $dm[1], $dm[2]);
}
}
}
# set object
$self->[5] = bless($derv, 'Math::Matrix');
}
# compute partial derivative function for _akima_para
# parameters: (d1, d2, m1, m2, ∂d1/∂n, ∂d2/∂n, ∂m1/∂n, ∂m2/∂n)
# returns; (∂s/∂n)
sub akima_spd {
=pod
https://www.wolframalpha.com/input/?i=d%2Fdx+(a(x)+*+c(x)+%2B+b(x)+*+d(x))%2F(a(x)+%2B+b(x))
d/dx((a(x) c(x) + b(x) d(x))/(a(x) + b(x))) =
(c(x) a'(x) + a(x) c'(x) + d(x) b'(x) + b(x) d'(x))/(a(x) + b(x)) - ((a'(x) + b'(x)) (a(x) c(x) + b(x) d(x)))/(a(x) + b(x))^2
=cut
# return partial derivative
return(($_[2] * $_[5] + $_[1] * $_[6] + $_[3] * $_[4] + $_[0] * $_[7])/($_[1] + $_[0]) - (($_[5] + $_[4]) * ($_[1] * $_[2] + $_[0] * $_[3]))/($_[1] + $_[0])**2);
}
# compute partial derivative function for _akima_para
# parameters: (m1, m2, ∂m1/∂n, ∂m2/∂n, d)
# returns; (∂d/∂n)
sub akima_dpd {
=pod
https://www.wolframalpha.com/input/?i=d%2Fdx+(sqrt((a(x)+-+b(x))%5E2+%2B+c))
d/dx(sqrt((a(x) - b(x))^2 + c)) =
((a(x) - b(x)) (a'(x) - b'(x)))/sqrt((a(x) - b(x))^2 + c)
=cut
# return partial derivative
return(($_[0] - $_[1]) * ($_[2] - $_[3])/$_[4]);
}
# add local min/max values
# parameter: (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);
# add min/max s-values
push(@s, map {$_ + $i} @t);
# add inner knot s-value, if derivative is 0
push(@s, $i) if ($i && $self->[3][$i] == 0);
# 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]];
}
# 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 ($dx, $y1, $y2, $m1, $m2);
my ($a, $b, $c, $dscr, @t);
# compute delta x-value (x-values : range)
$dx = ($#{$self->[1]} > 1) ? ($self->[1][$low + 1] - $self->[1][$low]) : ($self->[1][1] - $self->[1][0])/$#{$self->[2]};
# get endpoint values
$y1 = $self->[2][$low];
$y2 = $self->[2][$low + 1];
$m1 = $self->[3][$low] * $dx;
$m2 = $self->[3][$low + 1] * $dx;
# 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 ($dx, $h00, $h01, $h10, $h11);
# compute delta x-value (x-values : range)
$dx = ($#{$self->[1]} > 1) ? ($self->[1][$low + 1] - $self->[1][$low]) : ($self->[1][1] - $self->[1][0])/$#{$self->[2]};
# 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]/$dx + $h01 * $self->[2][$low + 1]/$dx + $h10 * $self->[3][$low] + $h11 * $self->[3][$low + 1])/$dx);
}
# compute derivative for unit spline segment
# parameters: (ref_to_object, t-value, segment)
# returns: (derivative)
sub _derv {
# get parameters
my ($self, $t, $low) = @_;
# local variables
my ($dx, $tc, $ttc, $h00, $h01, $h10, $h11);
# compute delta x-value (x-values : range)
$dx = ($#{$self->[1]} > 1) ? ($self->[1][$low + 1] - $self->[1][$low]) : ($self->[1][1] - $self->[1][0])/$#{$self->[2]};
# if position is non-zero
if ($t) {
# compute Hermite derivative coefficients
$tc = (1 - $t);
$ttc = -3 * $t * $tc;
$h00 = 2 * $ttc;
$h01 = -2 * $ttc;
$h10 = $ttc + $tc;
$h11 = $ttc + $t;
# interpolate value
return($h00 * $self->[2][$low]/$dx + $h01 * $self->[2][$low + 1]/$dx + $h10 * $self->[3][$low] + $h11 * $self->[3][$low + 1]);
} else {
# use lower knot derivative value
return($self->[3][$low]);
}
}
# compute transform value
# parameters: (ref_to_object, t-value, segment)
# returns: (y-value)
sub _fwd {
# get parameters
my ($self, $t, $low) = @_;
# local variables
my ($dx, $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) {
# compute delta x-value (x-values : range)
$dx = ($#{$self->[1]} > 1) ? ($self->[1][$low + 1] - $self->[1][$low]) : ($self->[1][1] - $self->[1][0])/$#{$self->[2]};
# interpolate value (Lapack XS)
# if ICC::Support::Lapack module is loaded
if ($lapack) {
return(ICC::Support::Lapack::hermite($t, $self->[2][$low], $self->[2][$low + 1], $dx * $self->[3][$low], $dx * $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 * $dx * $self->[3][$low] + $h11 * $dx * $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 ($dx, $y1, $y2, $m1, $m2);
my ($A, $B, $C, $D, $dscr, @t);
my ($d0, $d1, $cs, $cc, @r, $ccr, $sol, $lim0, $lim1);
# compute delta x-value (x-values : range)
$dx = ($#{$self->[1]} > 1) ? ($self->[1][$low + 1] - $self->[1][$low]) : ($self->[1][1] - $self->[1][0])/$#{$self->[2]};
# get endpoint values
$y1 = $self->[2][$low];
$y2 = $self->[2][$low + 1];
$m1 = $self->[3][$low] * $dx;
$m2 = $self->[3][$low + 1] * $dx;
# 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 (abs($dscr) < 5E-15) {
# 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 (abs($dscr) < 5E-15) {
# if ∆0 is zero
if (abs($d0) < 5E-15) {
# add solution (triple root)
push(@t, -$B/(3 * $A));
} else {
# add solutions (double root and single root)
push(@t, (9 * $A * $D - $B * $C)/(2 * $d0));
push(@t, (4 * $A * $B * $C - 9 * $A**2 * $D - $B**3)/($A * $d0));
}
} else {
# compute ∆1
$d1 = 2 * $B**3 - 9 * $A * $B * $C + 27 * $A**2 * $D;
# compute C (a complex number)
$cs = sqrt($d1**2 - 4 * $d0**3);
$cc = cbrt($d1 == $cs ? ($d1 + $cs)/2 : ($d1 - $cs)/2);
# compute cube roots of 1 (three complex numbers)
@r = root(1, 3);
# for each root
for my $i (0 .. 2) {
# multiply by cube root of 1
$ccr = $r[$i] * $cc;
# compute solution (a complex number)
$sol = ($B + $ccr + $d0/$ccr)/(-3 * $A);
# add solution, if real
push(@t, Re($sol)) if ($dscr > 0 || abs(Im($sol)) < 1E-15);
}
}
}
# set test limits to avoid duplicates at knots
$lim0 = ($in == $y1) ? 1E-14 : 0;
$lim1 = ($in == $y2) ? (1 - 1E-14) : 1;
# return valid solutions
return(sort {$a <=> $b} grep {$_ > $lim0 && $_ < $lim1} @t);
}
# test vector values
# returns true if values are unique
# parameter: (vector)
# returns: (flag)
sub _unique {
# get parameter
my $vec = shift();
# local variable
my (%x);
# make count hash
for (@{$vec}) {$x{$_}++}
# return
return(@{$vec} == keys(%x));
}
# sort object input/output values
# parameter: (ref_to_object)
sub _sort_values {
# get parameter
my $self = shift();
# local variable
my (@p);
# if input contains x-values
if ($#{$self->[1]} > 1) {
# warn if vectors not same length
($#{$self->[1]} == $#{$self->[2]}) or carp('input and output vectors not same length');
# for each point
for my $i (0 .. $#{$self->[1]}) {
# copy x-y values
$p[$i] = [$self->[1][$i], $self->[2][$i]];
}
# sort points by input value
@p = sort {$a->[0] <=> $b->[0]} @p;
# for each point
for my $i (0 .. $#{$self->[1]}) {
# copy x-y values
$self->[1][$i] = $p[$i]->[0];
$self->[2][$i] = $p[$i]->[1];
}
} else {
# if x-range is decreasing
if ($self->[1][0] > $self->[1][-1]) {
# reverse x and y vectors
@{$self->[1]} = reverse(@{$self->[1]});
@{$self->[2]} = reverse(@{$self->[2]});
}
}
}
# fix highlight region data
# corrects the print-specific
# problem of highlight dot retention
# samples below the limit (0 - 1) are fixed
# parameter: (ref_to_object, limit)
sub _fix_hl {
# get parameter
my ($self, $limit) = @_;
# local variable
my ($span, $lo, $hi, $n, @ix, @xr, @x, @y, $bern, @new);
# compute 'output' range
$span = ($self->[2][-1] - $self->[2][0]);
# set 'output' limits
$lo = $self->[2][0] + 0.002 * $span;
$hi = $self->[2][0] + 2 * $limit * $span;
# get upper index of 'output' array
$n = $#{$self->[2]};
# find samples within limits
@ix = grep {$self->[2][$_] >= $lo && $self->[2][$_] <= $hi} (1 .. $n);
# use at least 4 samples
@ix = ($ix[0] .. $ix[0] + 3) if (@ix < 4);
# if input is a range
if (@{$self->[1]} == 2) {
# interpolate x-values
@xr = map {(1 - $_/$n) * $self->[1][0] + ($_/$n) * $self->[1][1]} (0 .. $n);
# get 'fit' x-values
@x = @xr[@ix];
} else {
# get 'fit' x-values
@x = @{$self->[1]}[@ix];
}
# get 'fit' y-values
@y = @{$self->[2]}[@ix];
# make 'bern' object of degree 3, fitted to x-y values (requires Lapack module)
$bern = ICC::Support::bern->new({'input' => [$x[0], $x[-1]], 'output' => [(undef) x 4], 'fit' => [\@x, \@y]});
# set 'output' limit
$hi = $self->[2][0] + $limit * $span;
# for each sample
for my $i (reverse(0 .. $n)) {
# if 'output' value ≤ limit
if (@new || $self->[2][$i] <= $hi) {
# compute new 'output' value using 'bern' curve
push(@new, $bern->transform(@{$self->[1]} > 2 ? $self->[1][$i] : $xr[$i]));
}
}
# splice new values, saving originals in hash
$self->[0]{'hl_original'} = [splice(@{$self->[2]}, 0, @new, reverse(@new))];
# save x-axis intercept (highlight dot retention value)
$self->[0]{'hl_retention'} = $bern->inverse(0);
# compute knot derivatives
_objderv($self);
# add segment min/max y-values
_minmax($self);
}
# fix shadow region data
# corrects the print-specific problem of
# non-monotonic data is the shadow region
# samples above the limit (0 - 1) are fixed
# parameter: (ref_to_object, limit)
sub _fix_sh {
# get parameter
my ($self, $limit) = @_;
# local variable
my ($span, $lo, $n, @ix, @xr, @x, @y, $bern, @new);
# compute 'output' range
$span = ($self->[2][-1] - $self->[2][0]);
# set 'output' limit
$lo = $self->[2][-1] - 2 * (1 - $limit) * $span;
# get upper index of 'output' array
$n = $#{$self->[2]};
# find samples with 'output' > limit
@ix = grep {$self->[2][$_] >= $lo} (0 .. $n);
# use at least 4 samples
@ix = (($ix[-1] - 3) .. $ix[-1]) if (@ix < 4);
# if input is a range
if (@{$self->[1]} == 2) {
# interpolate x-values
@xr = map {(1 - $_/$n) * $self->[1][0] + ($_/$n) * $self->[1][1]} (0 .. $n);
# get 'fit' x-values
@x = @xr[@ix];
} else {
# get 'fit' x-values
@x = @{$self->[1]}[@ix];
}
# get 'fit' y-values
@y = @{$self->[2]}[@ix];
# make 'bern' object of degree 3, fitted to x-y values (requires Lapack module)
$bern = ICC::Support::bern->new({'input' => [$x[0], $x[-1]], 'output' => [(undef) x 4], 'fit' => [\@x, \@y]});
# set 'output' limit
$lo = $self->[2][0] + $limit * $span;
# for each sample
for my $i (0 .. $n) {
# if 'output' value ≥ limit
if (@new || $self->[2][$i] >= $lo) {
# compute new 'output' value using 'bern' curve
push(@new, $bern->transform(@{$self->[1]} > 2 ? $self->[1][$i] : $xr[$i]));
}
}
# splice new values, saving originals in hash
$self->[0]{'sh_original'} = [splice(@{$self->[2]}, -@new, @new, @new)];
# compute knot derivatives
_objderv($self);
# add segment min/max y-values
_minmax($self);
}
# set object contents from parameter hash
# parameters: (ref_to_object, ref_to_parameter_hash)
sub _new_from_hash {
# get parameters
my ($self, $hash) = @_;
# local variables
my ($in, $out, $type, $damp, $fix, $limit, $pok, $fok);
# if 'input' is defined
if (defined($in = $hash->{'input'})) {
# if 'input' is a vector (2 or more values)
if (ICC::Shared::is_num_vector($in) && 2 <= @{$in}) {
# set object 'input' array
$self->[1] = [@{$in}];
# if 'input' is an integer > 0 (2 or more values)
} elsif (Scalar::Util::looks_like_number($in) && $in == int($in) && $in > 0) {
# make identity curve
$self->[1] = [map {$_/$in} (0 .. $in)];
} else {
# error
croak("invalid 'input' values");
}
} else {
# use default
$self->[1] = [0, 1];
}
# if 'output' is defined
if (defined($out = $hash->{'output'})) {
# if 'output' is a vector (2 or more values)
if (ICC::Shared::is_num_vector($out) && 2 <= @{$out}) {
# set object 'output' array
$self->[2] = [@{$out}];
# if 'output' is an integer > 0 (2 or more values)
} elsif (Scalar::Util::looks_like_number($out) && $out == int($out) && $out > 0) {
# make identity curve
$self->[2] = [map {$_/$out} (0 .. $out)];
} else {
# error
croak("invalid 'output' values");
}
} else {
# use default
$self->[2] = [0, 1];
}
# if 'type' is defined
if (defined($type = $hash->{'type'})) {
# verify supported type
($type =~ m/^(natural|akima)$/) or croak('unsupported spline type');
# set object type in hash
$self->[0]{'type'} = $type;
}
# if 'damp' is defined
if (defined($damp = $hash->{'damp'})) {
# verify supported type
(Scalar::Util::looks_like_number($damp) && $damp >= 0) or croak('invalid akima damping value');
# set object type in hash
$self->[0]{'damp'} = $damp;
}
# verify number of input values
(@{$self->[1]} == 2 || @{$self->[1]} == @{$self->[2]}) or croak('wrong number of input values');
# sort input/output values
_sort_values($self);
# compute knot derivatives
_objderv($self);
# add segment min/max y-values
_minmax($self);
# if 'fix_hl' is defined
if (defined($fix = $hash->{'fix_hl'})) {
# verify polarity
($pok = $self->[2][-1] > $self->[2][0]) or carp("'fix_hl' requires increasing output values");
# verify fix limit (-1 to 1)
($fok = (Scalar::Util::looks_like_number($fix) && $fix >= -1 && $fix <= 1)) or carp("invalid 'fix_hl' limit");
# compute 'output' limit
$limit = $self->[2][0] * (1 - $fix) + $self->[2][-1] * $fix;
# fix the data if polarity and limit are OK, and spline is non-monotonic in highlights
_fix_hl($self, abs($fix)) if ($pok && $fok && ($fix < 0 || grep {$_ <= $limit} monotonic($self, 'output')));
}
# if 'fix_sh' is defined
if (defined($fix = $hash->{'fix_sh'})) {
# verify polarity
($pok = $self->[2][-1] > $self->[2][0]) or carp("'fix_sh' requires increasing output values");
# verify fix limit (-1 to 1)
($fok = (Scalar::Util::looks_like_number($fix) && $fix >= -1 && $fix <= 1)) or carp("invalid 'fix_sh' limit");
# compute 'output' limit
$limit = $self->[2][0] * (1 - $fix) + $self->[2][-1] * $fix;
# fix the data if polarity and limit are OK, and spline is non-monotonic in shadows
_fix_sh($self, abs($fix)) if ($pok && $fok && ($fix < 0 || grep {$_ >= $limit} monotonic($self, 'output')));
}
}
# write spline object to ICC profile
# note: writes an equivalent curv object
# parameters: (ref_to_object, ref_to_parent_object, file_handle, ref_to_tag_table_entry)
sub _writeICCspline {
# get parameters
my ($self, $parent, $fh, $tag) = @_;
# local variables
my ($n, $up, @table);
# get count value (defaults to 4096, the maximum allowed in an 'mft2' tag)
$n = $self->[0]{'curv_points'} // 4096;
# validate number of table entries
($n == int($n) && $n >= 2) or carp('invalid number of table entries');
# array upper index
$up = $n - 1;
# for each table entry
for my $i (0 .. $up) {
# transform value
$table[$i] = transform($self, $i/$up);
}
# seek start of tag
seek($fh, $tag->[1], 0);
# write tag type signature and count
print $fh pack('a4 x4 N', 'curv', $n);
# write array
print $fh pack('n*', map {$_ < 0 ? 0 : ($_ > 1 ? 65535 : $_ * 65535 + 0.5)} @table);
}
1;