File Coverage

blib/lib/ICC/Support/spline.pm

Criterion	Covered	Total	%
statement	19	687	2.7
branch	2	348	0.5
condition	1	175	0.5
subroutine	6	44	13.6
pod	15	19	78.9
total	43	1273	3.3

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;