package ICC::Support::bern;
use strict;
use Carp;
our $VERSION = 0.41;
# revised 2018-09-04
#
# Copyright © 2004-2020 by William B. Birkett
# inherit from Shared
use parent qw(ICC::Shared);
# enable static variables
use feature 'state';
# create new bern object
# hash keys are: 'input', 'output', 'fit', 'fix_hl', 'fix_sh'
# default values for 'input' and 'output' are []
# parameters: ([ref_to_attribute_hash])
# returns: (ref_to_object)
sub new {
# get object class
my $class = shift();
# create empty bern object
my $self = [
{}, # object header
[], # parameter array
[], # min/max t-values
[] # min/max i/o-values
];
# 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 $self = Storable::dclone(shift());
# swap 'input' and 'output' arrays
($self->[1][0], $self->[1][1]) = ($self->[1][1], $self->[1][0]);
($self->[2][0], $self->[2][1]) = ($self->[2][1], $self->[2][0]);
($self->[3][0], $self->[3][1]) = ($self->[3][1], $self->[3][0]);
# return inverted object
return($self);
}
# write bern 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 bern data to profile
goto &_writeICCbern;
}
# 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 4370)
my $n = $self->[0]{'curv_points'} // 4370;
# return size
return(12 + $n * 2);
}
# compute curve function
# parameters: (input_value)
# returns: (output_value)
sub transform {
# get parameters
my ($self, $in) = @_;
# verify input
(Scalar::Util::looks_like_number($in)) or croak('invalid transform input');
# compute intermediate value(s)
my @t = _rev($self->[1][0], $self->[2][0], $self->[3][0], $in);
# warn if multiple values
(@t == 1) or carp('multiple transform values');
# return first output value
return(_fwd($self->[1][1], $t[0]));
}
# compute inverse curve function
# parameters: (input_value)
# returns: (output_value)
sub inverse {
# get parameters
my ($self, $in) = @_;
# verify input
(Scalar::Util::looks_like_number($in)) or croak('invalid transform input');
# compute intermediate value(s)
my @t = _rev($self->[1][1], $self->[2][1], $self->[3][1], $in);
# warn if multiple values
(@t == 1) or carp('multiple transform values');
# return first output value
return(_fwd($self->[1][0], $t[0]));
}
# compute curve derivative
# parameters: (input_value)
# returns: (derivative_value)
sub derivative {
# get parameters
my ($self, $in) = @_;
# return forward derivative
return(_derivative($self, 0, $in));
}
# compute parametric partial derivatives
# parameters: (input_value)
# returns: (partial_derivative_array)
sub parametric {
# get parameters
my ($self, $in) = @_;
# return forward parametric partial derivatives
return(_parametric($self, 0, $in));
}
# 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 (between 1 and 21 numeric values)
(ICC::Shared::is_num_vector($_[0]) && 1 <= @{$_[0]} && 21 >= @{$_[0]}) or croak('invalid input values');
# store output values
$self->[1][0] = [@{$_[0]}];
# check for min/max values
_minmax($self, 0);
} elsif (@_) {
# error
carp('too many parameters');
}
# return array reference
return($self->[1][0]);
}
# 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 (between 1 and 21 numeric values)
(ICC::Shared::is_num_vector($_[0]) && 1 <= @{$_[0]} && 21 >= @{$_[0]}) or croak('invalid input values');
# store output values
$self->[1][1] = [@{$_[0]}];
# check for min/max values
_minmax($self, 1);
} elsif (@_) {
# error
carp('too many parameters');
}
# return array reference
return($self->[1][1]);
}
# 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);
# 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->[1][$i]} > 1 && ($x0 = $self->[1][$i][0]) != ($x1 = $self->[1][$i][-1])) {
# 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->[1][$i]} > 1 && $self->[1][$i][0] != $self->[1][$i][-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->[1][$i][0];
$x1 = $self->[1][$i][-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->[1][$i][0];
$x1 = $self->[1][$i][-1];
# if source is 'minmax'
} elsif (! ref($src) && $src eq 'minmax') {
# get all possible min/max values
@minmax = ($self->[1][$i][0], $self->[1][$i][-1], @{$self->[3][$i]});
# get min/max values
$x0 = List::Util::min(@minmax);
$x1 = List::Util::max(@minmax);
}
# 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->[1][$i][0] == $x0;
$f1 = $self->[1][$i][-1] == $x1;
# adjust Bernstein coefficients (y = mx + b)
@{$self->[1][$i]} = map {$m * $_ + $b} @{$self->[1][$i]};
# set endpoints (to ensure exact values)
$self->[1][$i][0] = $y0 if ($f0);
$self->[1][$i][-1] = $y1 if ($f1);
}
}
# re-calculate min/max values
_minmax($self);
# return
return($self);
}
# check if monotonic
# returns min/max values
# returns t-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 (@s);
# compute min/max values
_minmax($self);
# if format undefined
if (! defined($fmt)) {
# return sorted t-values ('input' and 'output')
return(sort {$a <=> $b} (@{$self->[2][0]}, @{$self->[2][1]}));
# if format 'input'
} elsif ($fmt eq 'input') {
# return input values
return(map {_fwd($self->[1][0], $_)} sort {$a <=> $b} (@{$self->[2][0]}, @{$self->[2][1]}));
# if format 'output'
} elsif ($fmt eq 'output') {
# return output values
return(map {_fwd($self->[1][1], $_)} sort {$a <=> $b} (@{$self->[2][0]}, @{$self->[2][1]}));
} 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
sub update {
# recompute min/max values
goto &_minmax;
}
# make table for profile 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][0])) {
# set format
$fmt2 = join(', ', ('%.3f') x @{$self->[1][0]});
# append parameter string
$s .= sprintf(" input parameters [$fmt2]\n", @{$self->[1][0]});
}
# if output_parameter_array defined
if (defined($self->[1][1])) {
# set format
$fmt2 = join(', ', ('%.3f') x @{$self->[1][1]});
# append parameter string
$s .= sprintf(" output parameters [$fmt2]\n", @{$self->[1][1]});
}
# return
return($s);
}
# find min/max values within the interval (0 - 1)
# used by '_rev' to determine solution regions
# called by the 'update' method
# i/o_index: 0 - input, 1 - output, undef - both
# parameters: (ref_to_object, [i/o_index])
sub _minmax {
# get parameters
my ($self, $io) = @_;
# local variables
my (@s, $t, $d);
my ($par, $fwd, $fwd2, $drv);
# for input and/or output curves
for my $i (defined($io) ? ($io ? 1 : 0) : (0, 1)) {
# initialize s-values array
@s = ();
# get ref to Bernstein parameters
$par = $self->[1][$i];
# if number of parameters > 2
if (@{$par} > 2) {
# compute forward differences
$fwd = [map {$par->[$_] - $par->[$_ - 1]} (1 .. $#{$par})];
# if forward differences have different signs (+/-)
if (grep {($fwd->[0] < 0) ^ ($fwd->[$_] < 0)} (1 .. $#{$fwd})) {
# compute second forward differences
$fwd2 = [map {$fwd->[$_] - $fwd->[$_ - 1]} (1 .. $#{$fwd})];
# compute derivative values over range (0 - 1)
$drv = [map {$#{$par} * _fwd($fwd, $_/(4 * $#{$par}))} (0 .. 4 * $#{$par})];
# for each pair of derivatives
for my $j (1 .. $#{$drv}) {
# if derivatives have different signs (+/-)
if (($drv->[$j - 1] < 0) ^ ($drv->[$j] < 0)) {
# estimate root from the derivative values
$t = ($j * $drv->[$j - 1] - ($j - 1) * $drv->[$j])/(($drv->[$j - 1] - $drv->[$j]) * 4 * $#{$par});
# loop until derivative approaches 0
while (abs($d = $#{$par} * _fwd($fwd, $t)) > 1E-9) {
# adjust root using Newton's method
$t -= $d/($#{$fwd} * $#{$par} * _fwd($fwd2, $t));
}
# push root on t-value array
push(@s, $t) if ($t > 0 && $t < 1);
}
}
}
}
# warn if root values were found
print "curve $i is not monotonic\n" if (@s && ! defined($self->[2][$i]));
# save root t-values
$self->[2][$i] = [@s];
# save root i/o-values
$self->[3][$i] = [map {_fwd($par, $_)} @s];
}
}
# combined transform
# nominal domain (0 - 1)
# parameters: (ref_to_object, direction, input_value)
# returns: (output_value)
sub _transform {
# get parameters
my ($self, $dir, $in) = @_;
# local variables
my ($p0, $p1, $p2, $p3, @t);
# verify input
(Scalar::Util::looks_like_number($in)) or croak('invalid transform input');
# get parameter arrays
$p0 = $self->[1][$dir];
$p1 = $self->[2][$dir];
$p2 = $self->[3][$dir];
$p3 = $self->[1][1 - $dir];
# compute intermediate value(s)
@t = _rev($p0, $p1, $p2, $in);
# warn if multiple values
(@t == 1) or carp('multiple transform values');
# return first output value
return(_fwd($p3, $t[0]));
}
# combined derivative
# nominal domain (0 - 1)
# parameters: (ref_to_object, direction, input_value)
# returns: (derivative_value)
sub _derivative {
# get parameters
my ($self, $dir, $in) = @_;
# local variables
my ($p0, $p1, $p2, $p3, @t);
my ($di, $do);
# verify input
(Scalar::Util::looks_like_number($in)) or croak('invalid transform input');
# get parameter arrays
$p0 = $self->[1][$dir];
$p1 = $self->[2][$dir];
$p2 = $self->[3][$dir];
$p3 = $self->[1][1 - $dir];
# compute intermediate value(s)
@t = _rev($p0, $p1, $p2, $in);
# warn if multiple values
(@t == 1) || print "multiple transform values\n";
# compute input derivative (using first value)
($di = _drv($p0, $t[0])) or croak('infinite derivative value');
# compute output derivative (using first value)
$do = _drv($p3, $t[0]);
# return combined derivative
return($do/$di);
}
# compute parametric partial derivatives
# direction: 0 - normal, 1 - inverse
# parameters: (ref_to_object, direction, input_value)
# returns: (partial_derivative_array)
sub _parametric {
# get parameters
my ($self, $dir, $in) = @_;
# local variables
my ($p0, $p1, $p2, $p3);
my (@t, $t, $di, $do, $bfi, $bfo, @pj);
# verify input
(Scalar::Util::looks_like_number($in)) or croak('invalid transform input');
# get parameter arrays
$p0 = $self->[1][$dir];
$p1 = $self->[2][$dir];
$p2 = $self->[3][$dir];
$p3 = $self->[1][1 - $dir];
# compute intermediate value(s)
@t = _rev($p0, $p1, $p2, $in);
# warn if multiple values
(@t == 1) || print "multiple transform values\n";
# use first value
$t = $t[0];
# compute input derivative
$di = _drv($p0, $t);
# compute input Bernstein polynomials
$bfi = @{$p0} ? _bernstein($#{$p0}, $t) : undef();
# compute output derivative
$do = _drv($p3, $t);
# compute output Bernstein polynomials
$bfo = @{$p3} ? _bernstein($#{$p3}, $t) : undef();
# initialize array
@pj = ();
# if there are input parameters
if (defined($bfi)) {
# for each parameter
for my $i (0 .. $#{$bfi}) {
# verify non-zero input derivative
($di) or croak('infinite Jacobian element');
# compute partial derivative
$pj[$i] = -$bfi->[$i] * $do/$di;
}
}
# add output parameters (if any)
push(@pj, @{$bfo}) if (defined($bfo));
# return array of partial derivatives
return(@pj);
}
# direct forward transform
# nominal domain (0 - 1)
# parameters: (parameter_array_reference, input_value)
# returns: (output_value)
sub _fwd {
# get parameters
my ($par, $in) = @_;
# if no parameters
if (! @{$par}) {
# return input
return($in);
# one parameter
} elsif (@{$par} == 1) {
# return constant
return($par->[0]);
# two parameters
} elsif (@{$par} == 2) {
# return linear interpolated value
return((1 - $in) * $par->[0] + $in * $par->[1]);
# three or more parameters
} else {
# check if ICC::Support::Lapack module is loaded
state $lapack = defined($INC{'ICC/Support/Lapack.pm'});
# if input < 0
if ($in < 0) {
# return extrapolated value
return($par->[0] + $in * _drv($par, 0));
# if input > 1
} elsif ($in > 1) {
# return extrapolated value
return($par->[-1] + ($in - 1) * _drv($par, 1));
} else {
# if Lapack module loaded
if ($lapack) {
# return Bernstein interpolated value
return(ICC::Support::Lapack::bernstein($in, $par));
} else {
# return Bernstein interpolated value
return(ICC::Shared::dotProduct(_bernstein($#{$par}, $in), $par));
}
}
}
}
# direct reverse transform
# nominal range (0 - 1)
# parameters: (parameter_array_reference, x-min/max_array_reference, y-min/max_array_reference, input_value)
# returns: (output_value -or- array_of_output_values)
sub _rev {
# get parameters
my ($par, $xminmax, $yminmax, $in) = @_;
# local variables
my (@xs, @ys, $slope, $ext, @sol);
my ($p, $loop, $xlo, $xhi, $xval, $yval);
# if no parameters
if (! @{$par}) {
# return input
return($in);
# one parameter
} elsif (@{$par} == 1) {
# warn if input not equal to constant
($in == $par->[0]) or carp('curve input not equal to constant parameter');
# return input
return($in);
# two parameters
} elsif (@{$par} == 2) {
# return linear interpolated value
return(($in - $par->[0])/($par->[1] - $par->[0]));
# three or more parameters
} else {
# setup segment arrays
@xs = (0, @{$xminmax}, 1);
@ys = ($par->[0], @{$yminmax}, $par->[-1]);
# if slope(0) not 0
if ($slope = _drv($par, 0)) {
# compute extrapolation from 0
$ext = ($in - $par->[0])/$slope;
# save if less than 0
push(@sol, $ext) if ($ext < 0);
}
# test first y-value
push(@sol, $xs[0]) if ($ys[0] == $in);
# for each array segment
for my $i (1 .. $#xs) {
# if input value falls within segment
if (($in > $ys[$i - 1] && $in < $ys[$i]) || ($in < $ys[$i - 1] && $in > $ys[$i])) {
# get segment polarity
$p = $ys[$i] > $ys[$i - 1];
# initialize limits
$xlo = $xs[$i - 1];
$xhi = $xs[$i];
# initialize loop counter
$loop = 0;
# bisection method loop
while ($xhi - $xlo > 1E-3 && $loop++ < 15) {
# compute midpoint
$xval = ($xlo + $xhi)/2;
# compute y-value
$yval = _fwd($par, $xval);
# if y-value > input -or- y-value < input
if ($p ? ($yval > $in) : ($yval < $in)) {
# set high limit
$xhi = $xval;
} else {
# set low limit
$xlo = $xval;
}
}
# initialize loop counter
$loop = 0;
# Newton's method loop
while (abs($in - $yval) > 1E-12 && $loop++ < 5 && ($slope = _drv($par, $xval))) {
# adjust x-value
$xval += ($in - $yval)/$slope;
# compute y-value
$yval = _fwd($par, $xval);
}
# print warning if failed to converge
printf("'_rev' function of 'bern' object failed to converge with input of %.2f\n", $in) if (abs($in - $yval) > 1E-12);
# save result
push(@sol, $xval);
}
# add segment y-value, if a solution
push(@sol, $xs[$i]) if ($ys[$i] == $in);
}
# if slope(1) not 0
if ($slope = _drv($par, 1)) {
# compute extrapolation from 1
$ext = ($in - $par->[-1])/$slope + 1;
# save if greater than 1
push(@sol, $ext) if ($ext > 1);
}
}
# return result (array or first solution)
return(wantarray ? @sol : $sol[0]);
}
# forward derivative
# nominal domain is (0 - 1)
# parameters: (parameter_array_reference, input_value)
# returns: (output_value)
sub _drv {
# get parameters
my ($par, $in) = @_;
# if no parameters
if (! @{$par}) {
# return identity derivative
return(1);
# one parameter
} elsif (@{$par} == 1) {
# return constant derivative
return(0);
# two parameters
} elsif (@{$par} == 2) {
# return linear derivative
return($par->[1] - $par->[0]);
# three or more parameters
} else {
# check if ICC::Support::Lapack module is loaded
state $lapack = defined($INC{'ICC/Support/Lapack.pm'});
# limit input value (0 - 1) to force linear extrapolation
$in = $in < 0 ? 0 : ($in > 1 ? 1 : $in);
# compute parameter differences
my $diff = [map {$par->[$_ + 1] - $par->[$_]} (0 .. $#{$par} - 1)];
# if Lapack module loaded
if ($lapack) {
# return Bernstein interpolated derivative
return($#{$par} * ICC::Support::Lapack::bernstein($in, $diff));
} else {
# return Bernstein interpolated derivative
return($#{$par} * ICC::Shared::dotProduct(_bernstein($#{$diff}, $in), $diff));
}
}
}
# compute Bernstein basis values
# parameters: (bernstein_degree, t)
# returns: (ref_to_bernstein_array)
sub _bernstein {
# get parameters
my ($degree, $t0) = @_;
# local variables
my ($bern, $t1, $m0, $m1, $n);
# binomial coefficients
state $pascal = [
[1],
[1, 1],
[1, 2, 1],
[1, 3, 3, 1],
[1, 4, 6, 4, 1],
[1, 5, 10, 10, 5, 1],
[1, 6, 15, 20, 15, 6, 1],
[1, 7, 21, 35, 35, 21, 7, 1],
[1, 8, 28, 56, 70, 56, 28, 8, 1],
[1, 9, 36, 84, 126, 126, 84, 36, 9, 1],
[1, 10, 45, 120, 210, 252, 210, 120, 45, 10, 1],
[1, 11, 55, 165, 330, 462, 462, 330, 165, 55, 11, 1],
[1, 12, 66, 220, 495, 792, 924, 792, 495, 220, 66, 12, 1],
[1, 13, 78, 286, 715, 1287, 1716, 1716, 1287, 715, 286, 78, 13, 1],
[1, 14, 91, 364, 1001, 2002, 3003, 3432, 3003, 2002, 1001, 364, 91, 14, 1],
[1, 15, 105, 455, 1365, 3003, 5005, 6435, 6435, 5005, 3003, 1365, 455, 105, 15, 1],
[1, 16, 120, 560, 1820, 4368, 8008, 11440, 12870, 11440, 8008, 4368, 1820, 560, 120, 16, 1],
[1, 17, 136, 680, 2380, 6188, 12376, 19448, 24310, 24310, 19448, 12376, 6188, 2380, 680, 136, 17, 1],
[1, 18, 153, 816, 3060, 8568, 18564, 31824, 43758, 48620, 43758, 31824, 18564, 8568, 3060, 816, 153, 18, 1],
[1, 19, 171, 969, 3876, 11628, 27132, 50388, 75582, 92378, 92378, 75582, 50388, 27132, 11628, 3876, 969, 171, 19, 1],
[1, 20, 190, 1140, 4845, 15504, 38760, 77520, 125970, 167960, 184756, 167960, 125970, 77520, 38760, 15504, 4845, 1140, 190, 20, 1]
];
# copy binomial coefficients for this degree
$bern = [@{$pascal->[$degree]}];
# compute input complement
$t1 = 1 - $t0;
# initialize multipliers
$m0 = $m1 = 1;
# get polynomial degree
$n = $#{$bern};
# for each polynomial
for (1 .. $n) {
# multiply
$bern->[$_] *= ($m0 *= $t0);
$bern->[$n - $_] *= ($m1 *= $t1);
}
# return array reference
return($bern);
}
# fit bern object to data
# uses LAPACK dgels function to perform a least-squares fit
# the parameters to be fitted have the value 'undef' and must all be within the
# output side of the parameter array, e.g. [[0, 100], [undef, undef, undef, undef]]
# fit data is an array containing input and output vectors
# parameters: (ref_to_object, fit_data_reference, ref_to_parameter_hash)
# returns: (dgels_info_value)
sub _fit {
# get parameters
my ($self, $fit, $hash) = @_;
# local variables
my ($m, $n, $d, $t, $p, $span, @sol);
my (@so, $outz, $bern, $nrhs, $info);
my ($outz2, $bern2, $fix_hl, $fix_sh, @c2);
# check if ICC::Support::Lapack module is loaded
state $lapack = defined($INC{'ICC/Support/Lapack.pm'});
# verify Lapack module loaded (need to add non-Lapack routine someday)
($lapack) or croak("'fit' option requires ICC::Support::Lapack module");
# get degree
$d = $#{$self->[1][1]};
# for each output parameter
for my $i (0 .. $d) {
# if parameter is undefined
if (! defined($self->[1][1][$i])) {
# push index
push(@so, $i);
# set value to 0
$self->[1][1][$i] = 0;
}
}
# if no 'undef' parameters
if (@so == 0) {
# warning
print "no 'undef' parameters to fit\n";
# return
return(0)
}
# check for min/max values
_minmax($self);
# for each sample
for my $i (0 .. $#{$fit->[0]}) {
# compute output difference
$outz->[$i] = [$fit->[1][$i] - transform($self, $fit->[0][$i])];
# get intermediate value
$t = _rev($self->[1][0], $self->[2][0], $self->[3][0], $fit->[0][$i]);
# compute Bernstein values for undefined parameters
$bern->[$i] = [@{_bernstein($d, $t)}[@so]];
}
# make copies
$outz2 = Storable::dclone($outz);
$bern2 = Storable::dclone($bern);
# get array sizes
$m = @{$bern};
$n = @{$bern->[0]};
$nrhs = @{$outz->[0]};
# fit the data, return if failed
($info = ICC::Support::Lapack::dgels('N', $m, $n, $nrhs, $bern, $m, $outz, $m)) && return($info);
# get curve polarity (0 - increasing, 1 - decreasing)
$p = ($span = ($so[-1] == $d ? $outz->[$n - 1][0] : $self->[1][1][-1]) - ($so[0] == 0 ? $outz->[0][0] : $self->[1][1][0])) < 0;
# get flags (indicating HL or SH contrast needs to be fixed)
$fix_hl = $hash->{'fix_hl'} && $so[0] == 0 && $so[1] == 1 && (($outz->[0][0] > $outz->[1][0]) ^ $p);
$fix_sh = $hash->{'fix_sh'} && $so[-1] == $d && $so[-2] == $d - 1 && (($outz->[$n - 1][0] < $outz->[$n - 2][0]) ^ $p);
# if degree >= 4 and HL or SH needs fixing
if ($d >= 4 && ($fix_hl || $fix_sh)) {
# for each sample
for my $i (0 .. $#{$bern2}) {
# if highlight fix
if ($fix_hl) {
# get first 2 Bernstein values
@c2 = splice(@{$bern2->[$i]}, 0, 2);
# combine to single value
unshift(@{$bern2->[$i]}, ($c2[0] + $c2[1]));
}
# if shadow fix
if ($fix_sh) {
# get last 2 Bernstein values
@c2 = splice(@{$bern2->[$i]}, -2);
# combine to single value
push(@{$bern2->[$i]}, ($c2[0] + $c2[1]));
}
}
# get array sizes
$m = @{$bern2};
$n = @{$bern2->[0]};
$nrhs = @{$outz2->[0]};
# fit the data, return if failed
($info = ICC::Support::Lapack::dgels('N', $m, $n, $nrhs, $bern2, $m, $outz2, $m)) && return($info);
# extract solution
@sol = map {$outz2->[$_][0]} (0 .. $n - 1);
# replace combined coefficient(s)
splice(@sol, 1, 0, $sol[0] + $span * 0.001) if ($fix_hl);
splice(@sol, -1, 0, $sol[-1] - $span * 0.001) if ($fix_sh);
} else {
# extract solution
@sol = map {$outz->[$_][0]} (0 .. $n - 1);
}
# copy solution to output array
@{$self->[1][1]}[@so] = @sol;
# return
return(0);
}
# 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, $fit);
# if 'input' is defined
if (defined($in = $hash->{'input'})) {
# if 'input' is a vector (21 values or less, some may be 'undef')
if (ICC::Shared::is_vector($in) && 21 >= @{$in}) {
# set object 'input' array
$self->[1][0] = [@{$in}];
# if 'input' is an integer between 1 and 20
} elsif (Scalar::Util::looks_like_number($in) && $in == int($in) && $in > 0 && $in < 21) {
# make identity curve
$self->[1][0] = [map {$_/$in} (0 .. $in)];
} else {
# error
croak("invalid 'input' values");
}
} else {
# disable input
$self->[1][0] = [];
}
# if 'output' is defined
if (defined($out = $hash->{'output'})) {
# if 'output' is a vector (21 values or less, some may be 'undef')
if (ICC::Shared::is_vector($out) && 21 >= @{$out}) {
# set object 'output' array
$self->[1][1] = [@{$out}];
# if 'output' is an integer between 1 and 20
} elsif (Scalar::Util::looks_like_number($out) && $out == int($out) && $out > 0 && $out < 21) {
# make identity curve
$self->[1][1] = [map {$_/$out} (0 .. $out)];
} else {
# error
croak("invalid 'output' values");
}
} else {
# disable output
$self->[1][1] = [];
}
# if 'fit' is defined and not empty
if (defined($fit = $hash->{'fit'}) && @{$fit}) {
# verify 'fit' vectors
(ICC::Shared::is_num_vector($fit->[0]) && 1 <= @{$fit->[0]}) or croak("invalid 'fit' input values");
(ICC::Shared::is_num_vector($fit->[1]) && 1 <= @{$fit->[1]}) or croak("invalid 'fit' output values");
(@{$fit->[0]} == @{$fit->[1]}) or croak("'fit' input and output vectors are different size");
# fit object to data
(_fit($self, $fit, $hash)) && croak("'fit' operation failed");
} else {
# verify object parameters are all defined
(@{$self->[1][0]} == grep {defined($_)} @{$self->[1][0]}) or croak("some 'input' values are 'undef'");
(@{$self->[1][1]} == grep {defined($_)} @{$self->[1][1]}) or croak("some 'output' values are 'undef'");
}
# add segment min/max y-values
_minmax($self);
}
# write bern 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 _writeICCbern {
# 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;