File Coverage

lib/ICC/Support/bern.pm

Criterion	Covered	Total	%
statement	19	458	4.1
branch	2	230	0.8
condition	1	115	0.8
subroutine	6	33	18.1
pod	15	17	88.2
total	43	853	5.0

package ICC::Support::bern;

use strict;
use Carp;

our $VERSION = 0.41;

# revised 2018-09-04
#
# Copyright © 2004-2019 by William B. Birkett

# add development directory
use lib 'lib';

# 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;