/* Various operators as functions.
 */

logical_and a b = a && b;
logical_or a b = a || b;
bitwise_and a b = a & b;
bitwise_or a b = a | b;
eor a b = a ^ b;
left_shift a b = a << b;
right_shift a b = a >> b;
not a = !a;

less a b = a < b;
more a b = a > b;
less_equal a b = a <= b;
more_equal a b = a >= b;
equal a b = a == b;
not_equal a b = a != b;
pointer_equal a b = a === b;
not_pointer_equal a b = a !== b;

add a b = a + b;
subtract a b = a - b;
multiply a b = a * b;
divide a b = a / b;
power a b = a ** b;
square x = x * x;
remainder a b = a % b;

cons a b = a : b;
join a b = a ++ b;
subscript a b = a ? b;

generate s n f = [s, n .. f];
comma r i = (r, i);

compose f g = f @ g;

cast_unsigned_char x = (unsigned char) x;
cast_signed_char x = (signed char) x;
cast_unsigned_short x = (unsigned short) x;
cast_signed_short x = (signed short) x;
cast_unsigned_int x = (unsigned int) x;
cast_signed_int x = (signed int) x;
cast_float x = (float) x;
cast_double x = (double) x;
cast_complex x = (complex) x;
cast_double_complex x = (double complex) x;

unary_minus x = -x;
negate x = !x;
complement x = ~x;
unary_plus x = +x;

if_then_else a b c = if a then b else c;

// the vector ops ... im is an image, vec is a real_list
vec op_name im vec
	= im_lintra_vec ones im vec,
		op_name == "add" || op_name == "add'"
	= im_lintra_vec ones (-1 * im) vec,
		op_name == "subtract'" 
	= im_lintra_vec ones im inv,
		op_name == "subtract" 
	= im_lintra_vec vec im zeros,
		op_name == "multiply" || op_name == "multiply'"
	= im_lintra_vec vec (1 / im) zeros,
		op_name == "divide'" 
	= im_lintra_vec recip im zeros,
		op_name == "divide" 
	= im_expntra_vec im vec,
		op_name == "power'" 
	= im_powtra_vec im vec,
		op_name == "power" 
	= im_remainderconst_vec im vec,
		op_name == "remainder" 
	= im_andimage_vec im vec,
		op_name == "bitwise_and" || op_name == "bitwise_and'"
	= im_orimage_vec im vec,
		op_name == "bitwise_or" || op_name == "bitwise_or'"
	= im_eorimage_vec im vec,
		op_name == "eor" || op_name == "eor'"
	= im_equal_vec im vec,
		op_name == "equal" || op_name == "equal'"
	= im_notequal_vec im vec,
		op_name == "not_equal" || op_name == "not_equal'"
	= im_less_vec im vec,
		op_name == "less" 
	= im_moreeq_vec im vec,
		op_name == "less'" 
	= im_lesseq_vec im vec,
		op_name == "less_equal" 
	= im_more_vec im vec,
		op_name == "less_equal'" 
	= error "unimplemented vector operation"
{
	zeros = replicate (len vec) 0;
	ones = replicate (len vec) 1;
	recip = map (divide 1) vec;
	inv = map (multiply (-1)) vec;
}

/* Macbeth chart patch names.
 */
_macbeth_names = [
	"Dark skin",
	"Light skin",
	"Blue sky",
	"Foliage",
	"Blue flower",
	"Bluish green",
	"Orange",
	"Purplish blue",
	"Moderate red",
	"Purple",
	"Yellow green",
	"Orange yellow",
	"Blue",
	"Green",
	"Red",
	"Yellow",
	"Magenta",
	"Cyan",
	"White (density 0.05)",
	"Neutral 8 (density 0.23)",
	"Neutral 6.5 (density 0.44)",
	"Neutral 5 (density 0.70)",
	"Neutral 3.5 (density 1.05)",
	"Black (density 1.50)"
];

bandsplit x 
	= oo_unary_function bandsplit_op x, is_class x
	= map (subscript x) [0 .. bands - 1], is_image x
	= error (_ "bad arguments to " ++ "bandsplit")
{
	bands = im_header_int "Bands" x;
	bandsplit_op = Operator "bandsplit" (map Image @ bandsplit)
		Operator_type.COMPOUND false;
}

bandjoin l
	= wrapper joined,
		has_wrapper
	= joined, is_listof has_image l
	= error (_ "bad arguments to " ++ "bandjoin")
{
	has_wrapper = has_member_list (has_member "Image") l;
	wrapper = get_member_list (has_member "Image") (get_member "Image") l;
	joined = im_gbandjoin (map get_image l);
}

mean x
	= oo_unary_function mean_op x, is_class x
	= im_avg x, is_image x
	= mean_list x, is_real_list x || is_matrix x
	= error (_ "bad arguments to " ++ "mean")
{
	mean_op = Operator "mean" mean_object Operator_type.COMPOUND false;

	mean_object x
		= im_avg x, is_image x
		= mean_list x, is_real_list x || is_matrix x
		= error (_ "bad arguments to " ++ "mean");

	mean_list l 
		= s / n
	{
		totals = sum l;
		n = totals?0;
		s = totals?1;
	}

	// return [n, sum] for a list of numbers, or a list of list of num
	// etc.
	sum x
		= foldr accumulate [0, 0] x
	{
		accumulate x sofar
			= [n + 1, x + s], is_real x
			= [n + n', s + s'], is_list x
			= error "mean_list: not real or [real]"
		{
			n = sofar?0;
			s = sofar?1;

			sub_acc = sum x;

			n' = sub_acc?0;
			s' = sub_acc?1;
		}
	}
}

deviation x
	= oo_unary_function deviation_op x, is_class x
	= im_deviate x, is_image x
	= deviation_list x, is_real_list x || is_matrix x
	= error (_ "bad arguments to " ++ "deviation")
{
	deviation_op = Operator "deviation" 
		deviation_object Operator_type.COMPOUND false;

	deviation_object x
		= im_deviate x, is_image x
		= deviation_list x, is_real_list x || is_matrix x
		= error (_ "bad arguments to " ++ "deviation");

	deviation_list l 
		= (abs (s2 - (s * s / n)) / (n - 1)) ** 0.5
	{
		totals = sum_sum2_list l;
		n = totals?0;
		s = totals?1;
		s2 = totals?2;
	}

	// return n, sum, sum of squares for a list of reals
	sum_sum2_list x
		= foldr accumulate [0, 0, 0] x
	{
		accumulate x sofar
			= [n + 1, x + s, x * x + s2], is_real x
			= [n + n', s + s', s2 + s2'], is_list x
			= error "sum_sum2_list: not real or [real]"
		{
			n = sofar?0;
			s = sofar?1;
			s2 = sofar?2;

			sub_acc = sum_sum2_list x;

			n' = sub_acc?0;
			s' = sub_acc?1;
			s2' = sub_acc?2;
		}
	}
}

abs x
	= oo_unary_function abs_op x, is_class x
	= im_abs x, is_image x
	= abs_cmplx x, is_complex x
	= abs_num x, is_real x
	= error (_ "bad arguments to " ++ "abs")
{
	abs_op = Operator "abs" abs_object Operator_type.COMPOUND false;

	abs_object x
		= im_abs x, is_image x
		= abs_cmplx x, is_complex x
		= abs_num x, is_real x
		= abs_list x, is_real_list x
		= abs_list (map abs_list x), is_matrix x
		= error (_ "bad arguments to " ++ "abs");

	abs_list l = (foldr1 add (map square l)) ** 0.5;

	abs_num n
		= n, n >= 0
		= -n;

	abs_cmplx c = ((re c)**2 + (im c)**2) ** 0.5;
}

copy x
	= oo_unary_function copy_op x, is_class x
	= im_copy x, is_image x
	= x
{
	copy_op = Operator "copy" copy Operator_type.COMPOUND_REWRAP false;
}

// like abs, but treat pixels as vectors ... ie. always get a 1-band image
// back ... also treat matricies as lists of vectors
// handy for dE from object difference
abs_vec x
	= oo_unary_function abs_vec_op x, is_class x
	= abs_vec_image x, is_image x
	= abs_vec_cmplx x, is_complex x
	= abs_vec_num x, is_real x
	= error (_ "bad arguments to " ++ "abs_vec")
{
	abs_vec_op = Operator "abs_vec" 
		abs_vec_object Operator_type.COMPOUND false;

	abs_vec_object x
		= abs_vec_image x, is_image x
		= abs_vec_cmplx x, is_complex x
		= abs_vec_num x, is_real x
		= abs_vec_list x, is_real_list x
		= mean (Vector (map abs_vec_list x)), is_matrix x
		= error (_ "bad arguments to " ++ "abs_vec");

	abs_vec_list l = (foldr1 add (map square l)) ** 0.5;

	abs_vec_num n
		= n, n >= 0
		= -n;

	abs_vec_cmplx c = ((re c)**2 + (im c)**2) ** 0.5;

	abs_vec_image im 
		= (foldr1 add (map square (bandsplit im))) ** 0.5;
}

transpose x
	= oo_unary_function transpose_op x, is_class x
	= transpose_image x, is_image x
	= transpose_matrix x, is_list x && is_list (hd x)
	= error (_ "bad arguments to " ++ "transpose")
{
	transpose_op = Operator "transpose" 
		transpose_object Operator_type.COMPOUND_REWRAP false;

	transpose_object x
		= transpose_matrix x, is_matrix x
		= transpose_image x, is_image x
		= error (_ "bad arguments to " ++ "transpose");

	transpose_matrix l
		= [], l' == []
		= (map hd l') : (transpose_matrix (map tl l'))
	{
		l' = takewhile (not_equal []) l;
	}

	transpose_image = im_flipver @ im_rot270;
}

rot45 x
	= oo_unary_function rot45_op x, is_class x
	= error "rot45 image: not implemented", is_image x 
	= error (_ "bad arguments to " ++ "rot45")
{
	rot45_op = Operator "rot45" 
		rot45_object Operator_type.COMPOUND_REWRAP false;

	rot45_object x
		= rot45_matrix x, is_odd_square_matrix x
		= error "rot45 image: not implemented", is_image x 
		= error (_ "bad arguments to " ++ "rot45");

	// slow, but what the heck
	rot45_matrix l = (im_rotate_dmask45 (Matrix l)).value;
}

// apply an image function to a [[real]] ... matrix is converted to a 1 band
// image for processing
apply_matrix_as_image fn m
	= (get_value @ im_vips2mask @ fn @ im_mask2vips @ Matrix) m;

rot90 x
	= oo_unary_function rot90_op x, is_class x
	= im_rot90 x, is_image x
	= error (_ "bad arguments to " ++ "rot90")
{
	rot90_op = Operator "rot90" 
		rot90_object Operator_type.COMPOUND_REWRAP false;

	rot90_object x
		= rot90_matrix x, is_matrix x
		= im_rot90 x, is_image x
		= error (_ "bad arguments to " ++ "rot90");

	// slow, but what the heck
	// avoid im_rotate_dmask90(), it can only do square odd-sided matricies
	rot90_matrix l = apply_matrix_as_image im_rot90 l;
}

rot180 x
	= oo_unary_function rot180_op x, is_class x
	= im_rot180 x, is_image x
	= error (_ "bad arguments to " ++ "rot180")
{
	rot180_op = Operator "rot180" 
		rot180_object Operator_type.COMPOUND_REWRAP false;

	rot180_object x
		= rot180_matrix x, is_matrix x
		= im_rot180 x, is_image x
		= error (_ "bad arguments to " ++ "rot180");

	// slow, but what the heck
	rot180_matrix l = apply_matrix_as_image im_rot180 l;
}

rot270 x
	= oo_unary_function rot270_op x, is_class x
	= im_rot270 x, is_image x
	= error (_ "bad arguments to " ++ "rot270")
{
	rot270_op = Operator "rot270" 
		rot270_object Operator_type.COMPOUND_REWRAP false;

	rot270_object x
		= rot270_matrix x, is_matrix x
		= im_rot270 x, is_image x
		= error (_ "bad arguments to " ++ "rot270");

	// slow, but what the heck
	rot270_matrix l = apply_matrix_as_image im_rot270 l;
}

image_set_type type x 
	= oo_unary_function image_set_type_op x, is_class x
	= im_copy_set x (to_real type) 
		(im_header_double "Xres" x) (im_header_double "Yres" x)
		(im_header_int "Xoffset" x) (im_header_int "Yoffset" x),
		is_image x
	= error (_ "bad arguments to " ++ "image_set_type:" ++ 
		print type ++ " " ++ print x)
{
	image_set_type_op = Operator "image_set_type" 
		(image_set_type type) Operator_type.COMPOUND_REWRAP false;
}

image_set_origin xoff yoff x 
	= oo_unary_function image_set_origin_op x, is_class x
	= im_copy_set x 
		(im_header_int "Type" x) 
		(im_header_double "Xres" x) (im_header_double "Yres" x)
		(to_real xoff) (to_real yoff),
		is_image x
	= error (_ "bad arguments to " ++ "image_set_origin")
{
	image_set_origin_op = Operator "image_set_origin" 
		(image_set_origin xoff yoff) 
		Operator_type.COMPOUND_REWRAP false;
}

rotquad x
	= oo_unary_function rotquad_op x, is_class x
	= im_rotquad x, is_image x
	= error (_ "bad arguments to " ++ "rotquad")
{
	rotquad_op = Operator "rotquad" 
		rotquad_object Operator_type.COMPOUND_REWRAP false;

	rotquad_object x
		= rotquad_matrix x, is_matrix x
		= im_rotquad x, is_image x
		= error (_ "bad arguments to " ++ "rotquad");

	rotquad_matrix l = apply_matrix_as_image im_rotquad l;
}

cache tile_width tile_height max_tiles x
	= oo_unary_function cache_op x, is_class x
	= im_cache x (to_real tile_width) (to_real tile_height) 
		(to_real max_tiles), is_image x
	= error (_ "bad arguments to " ++ "cache")
{
	cache_op = Operator "cache" 
		(cache tile_width tile_height max_tiles) 
		Operator_type.COMPOUND_REWRAP false;
}

tile across down x
	= oo_unary_function tile_op x, is_class x
	= im_replicate x (to_real across) (to_real down), is_image x
	= error (_ "bad arguments to " ++ "tile")
{
	tile_op = Operator "tile" 
		(tile across down) Operator_type.COMPOUND_REWRAP false;
}

fliptb x
	= oo_unary_function fliptb_op x, is_class x
	= im_flipver x, is_image x
	= error (_ "bad arguments to " ++ "fliptb")
{
	fliptb_op = Operator "fliptb" 
		fliptb_object Operator_type.COMPOUND_REWRAP false;

	fliptb_object x
		= fliptb_matrix x, is_matrix x
		= im_flipver x, is_image x
		= error (_ "bad arguments to " ++ "fliptb");

	fliptb_matrix l = reverse l;
}

fliplr x
	= oo_unary_function fliplr_op x, is_class x
	= im_fliphor x, is_image x
	= error (_ "bad arguments to " ++ "fliplr")
{
	fliplr_op = Operator "fliplr" 
		fliplr_object Operator_type.COMPOUND_REWRAP false;

	fliplr_object x
		= fliplr_matrix x, is_matrix x
		= im_fliphor x, is_image x
		= error (_ "bad arguments to " ++ "fliplr");

	fliplr_matrix l = map reverse l;
}

max_pair a b
	= a, a > b
	= b;

min_pair a b
      =	a, a < b
      =	b;

range min value max = min_pair max (max_pair min value);

max x
	= oo_unary_function max_op x, is_class x
	= im_max x, is_image x
	= max_list x, is_real_list x || is_matrix x
	= x, is_number x
	= error (_ "bad arguments to " ++ "max")
{
	max_op = Operator "max" max Operator_type.COMPOUND false;

	max_list x
		= foldr1 max_pair x, is_real_list x
		= foldr1 max_pair (map max_list x), is_matrix x
		= max x;
}

min x
	= oo_unary_function min_op x, is_class x
	= im_min x, is_image x
	= min_list x, is_real_list x || is_matrix x
	= x, is_number x
	= error (_ "bad arguments to " ++ "min")
{
	min_op = Operator "min" min Operator_type.COMPOUND false;

	min_list x
		= foldr1 min_pair x, is_real_list x
		= foldr1 min_pair (map min_list x), is_matrix x
		= min x;
}

maxpos x
	= oo_unary_function maxpos_op x, is_class x
	= im_maxpos x, is_image x
	= maxpos_matrix x, is_matrix x
	= error (_ "bad arguments to " ++ "maxpos")
{
	maxpos_op = Operator "maxpos" maxpos Operator_type.COMPOUND false;

	maxpos_matrix m 
		= (indexes?row, row)
	{
		max_value = max (Matrix m);
		indexes = map (index (equal max_value)) m;
		row = index (not_equal (-1)) indexes;
	}
}

minpos x
	= oo_unary_function minpos_op x, is_class x
	= im_minpos x, is_image x
	= minpos_matrix x, is_matrix x
	= error (_ "bad arguments to " ++ "minpos")
{
	minpos_op = Operator "minpos" minpos Operator_type.COMPOUND false;

	minpos_matrix m 
		= (indexes?row, row)
	{
		min_value = min (Matrix m);
		indexes = map (index (equal min_value)) m;
		row = index (not_equal (-1)) indexes;
	}
}

stats x
	= oo_unary_function stats_op x, is_class x
	= im_stats x, is_image x
	= im_stats (to_image x).value, is_matrix x
	= error (_ "bad arguments to " ++ "stats")
{
	stats_op = Operator "stats" 
		stats Operator_type.COMPOUND false;
}

e = 2.7182818284590452354;

pi = 3.14159265358979323846;

rad d = 2 * pi * (d / 360);

deg r = 360 * r / (2 * pi);

sign x
	= oo_unary_function sign_op x, is_class x
	= im_sign x, is_image x
	= sign_cmplx x, is_complex x
	= sign_num x, is_real x
	= error (_ "bad arguments to " ++ "sign")
{
	sign_op = Operator "sign" sign Operator_type.COMPOUND_REWRAP false;

	sign_num n
		= 0, n == 0
		= 1, n > 0
		= -1;

	sign_cmplx c 
		= (0, 0), mod == 0
		= (re c / mod, im c / mod)
	{
		mod = abs c;
	}
}

rint x
	= oo_unary_function rint_op x, is_class x
	= im_rint x, is_image x 
	= rint_value x, is_number x 
	= error (_ "bad arguments to " ++ "rint")
{
	rint_op = Operator "rint" rint Operator_type.ARITHMETIC false;

	rint_value x
		= (int) (x + 0.5), x > 0
		= (int) (x - 0.5);
}

scale x
	= oo_unary_function scale_op x, is_class x
	= (unsigned char) x, is_number x 
	= im_scale x, is_image x 
	= scale_list x, is_real_list x || is_matrix x
	= error (_ "bad arguments to " ++ "scale")
{
	scale_op = Operator "scale" scale Operator_type.COMPOUND_REWRAP false;

	scale_list l
		= apply_scale s o l
	{
		mn = find_limit min_pair l;
		mx = find_limit max_pair l;
		s = 255.0 / (mx - mn);
		o = -(mn * s);
	}

	find_limit fn l
		= find_limit fn (map (find_limit fn) l), is_listof is_list l
		= foldr1 fn l;

	apply_scale s o x
		= x * s + o, is_number x
		= map (apply_scale s o) x;
}

scaleps x
	= oo_unary_function scale_op x, is_class x
	= im_scaleps x, is_image x 
	= error (_ "bad arguments to " ++ "scale")
{
	scale_op = Operator "scaleps" 
		scaleps Operator_type.COMPOUND_REWRAP false;
}

fwfft x
	= oo_unary_function fwfft_op x, is_class x
	= im_fwfft x, is_image x 
	= error (_ "bad arguments to " ++ "fwfft")
{
	fwfft_op = Operator "fwfft" 
		fwfft Operator_type.COMPOUND_REWRAP false;
}

invfft x
	= oo_unary_function invfft_op x, is_class x
	= im_invfftr x, is_image x 
	= error (_ "bad arguments to " ++ "invfft")
{
	invfft_op = Operator "invfft" 
		invfft Operator_type.COMPOUND_REWRAP false;
}

falsecolour x
	= oo_unary_function falsecolour_op x, is_class x
	= image_set_type Image_type.sRGB (im_falsecolour x), is_image x 
	= error (_ "bad arguments to " ++ "falsecolour")
{
	falsecolour_op = Operator "falsecolour" 
		falsecolour Operator_type.COMPOUND_REWRAP false;
}

polar x
	= oo_unary_function polar_op x, is_class x
	= im_c2amph x, is_image x
	= polar_cmplx x, is_complex x
	= error (_ "bad arguments to " ++ "polar")
{
	polar_op = Operator "polar" polar Operator_type.COMPOUND false;

	polar_cmplx r
		= (l, a)
	{
		a
			= 270, x == 0 && y < 0
			= 90, x == 0 && y >= 0
			= 360 + atan (y / x), x > 0 && y < 0
			= atan (y / x), x > 0 && y >= 0
			= 180 + atan (y / x);

		l = (x ** 2 + y ** 2) ** 0.5;

		x = re r;
		y = im r;
	}
}

rectangular x
	= oo_unary_function rectangular_op x, is_class x
	= im_c2rect x, is_image x
	= rectangular_cmplx x, is_complex x
	= error (_ "bad arguments to " ++ "rectangular")
{
	rectangular_op = Operator "rectangular" 
		rectangular Operator_type.COMPOUND false;

	rectangular_cmplx p
		= (x, y)
	{
		l = re p;
		a = im p;

		x = l * cos a;
		y = l * sin a;
	}
}

recomb matrix image
	= colour_unary recomb_op image
{
	recomb_op x
		= im_recomb x (to_matrix matrix), is_image x
		= error (_ "bad arguments to " ++ "recomb");
}

extract_area x y w h obj
	= oo_unary_function extract_area_op obj, is_class obj
	= im_extract_area obj x' y' w' h', is_image obj
	= map (extract_range x' w') (extract_range y' h' obj), is_matrix obj
	= error (_ "bad arguments to " ++ "extract_area")
{
	x' = to_real x;
	y' = to_real y;
	w' = to_real w;
	h' = to_real h;

	extract_area_op = Operator "extract_area" (extract_area x y w h)
		Operator_type.COMPOUND_REWRAP false;

	extract_range from length list
		= (take length @ drop from) list;
}

extract_band b obj = subscript obj b;

extract_row y obj
	= oo_unary_function extract_row_op obj, is_class obj
	= extract_area 0 y' (get_width obj) 1 obj, is_image obj
	= [obj?y'], is_matrix obj
	= error (_ "bad arguments to " ++ "extract_row")
{
	y' = to_real y;

	extract_row_op = Operator "extract_row" (extract_row y)
		Operator_type.COMPOUND_REWRAP false;
}

extract_column x obj
	= oo_unary_function extract_column_op obj, is_class obj
	= extract_area x' 0 1 height obj, is_image obj
	= map (converse cons [] @ converse subscript x') obj, is_matrix obj
	= error (_ "bad arguments to " ++ "extract_column")
{
	x' = to_real x;
	height = im_header_int "Ysize" obj;

	extract_column_op = Operator "extract_column" (extract_column x)
		Operator_type.COMPOUND_REWRAP false;
}

blend cond in1 in2
	= oo_binary_function blend_op cond [in1,in2], is_class cond
	= im_blend (get_image cond) (get_image in1) (get_image in2),
		has_image cond && has_image in1 && has_image in2
	= error (_ "bad arguments to " ++ "blend")
{
	blend_op = Operator "blend" 
		blend_obj Operator_type.COMPOUND_REWRAP false;

	blend_obj cond x
		= blend_result_image
	{
		then_part = x?0;
		else_part = x?1;

		// get things about our output from inputs in this order
		objects = [then_part, else_part, cond];

		// properties of our output image
		target_width = get_member_list has_width get_width objects;
		target_height = get_member_list has_height get_height objects;
		target_bands = get_member_list has_bands get_bands objects;
		target_format = get_member_list has_format get_format objects;
		target_type = get_member_list has_type get_type objects;

		to_image x
			= x, is_image x
			= x.value, is_Image x
			= black + x
		{
			black = im_black target_width target_height target_bands;
		}

		then_image = to_image then_part;
		else_image = to_image else_part;

		then_image' = clip2fmt target_format then_image;
		else_image' = clip2fmt target_format else_image;

		resized = size_alike [cond, then_image', else_image'];

		blend_result_image = image_set_type target_type
			(im_blend resized?0 resized?1 resized?2);
	}
}

insert x y small big
	= oo_binary_function insert_op small big, is_class small
	= oo_binary'_function insert_op small big, is_class big
	= im_insert big small (to_real x) (to_real y), 
		is_image small && is_image big
	= error (_ "bad arguments to " ++ "insert")
{
	insert_op = Operator "insert" 
		(insert x y) Operator_type.COMPOUND_REWRAP false;
}

insert_noexpand x y small big
	= oo_binary_function insert_noexpand_op small big, is_class small
	= oo_binary'_function insert_noexpand_op small big, is_class big
	= im_insert_noexpand big small (to_real x) (to_real y), 
		is_image small && is_image big
	= error (_ "bad arguments to " ++ "insert_noexpand")
{
	insert_noexpand_op = Operator "insert_noexpand" 
		(insert_noexpand x y) Operator_type.COMPOUND_REWRAP false;
}

measure x y w h u v image
	= oo_unary_function measure_op image, is_class image
	= im_measure image 
		(to_real x) (to_real y) (to_real w) (to_real h)
		(to_real u) (to_real v), 
			is_image image
	= error (_ "bad arguments to " ++ "measure")
{
	measure_op = Operator "measure" 
		(measure x y w h u v) Operator_type.COMPOUND_REWRAP false;
}

extract_bands b n obj 
	= oo_unary_function extract_bands_op obj, is_class obj
	= im_extract_bands obj (to_real b) (to_real n), is_image obj
	= error (_ "bad arguments to " ++ "extract_bands")
{
	extract_bands_op = Operator "extract_bands" 
		(extract_bands b n) Operator_type.COMPOUND_REWRAP false;
}

transform ipol wrap params image
	= oo_unary_function transform_op image, is_class image
	= im_transform image 
		(to_matrix params) (to_real ipol) (to_real wrap), is_image image
	= error (_ "bad arguments to " ++ "transform")
{
	transform_op = Operator "transform" 
		(transform ipol wrap params) 
		Operator_type.COMPOUND_REWRAP false;
}

transform_search max_error max_iterations order ipol wrap sample reference
	= oo_binary_function transform_search_op sample reference, is_class sample
	= oo_binary'_function transform_search_op sample reference, 
		is_class reference
	= im_transform_search sample reference
		(to_real max_error) (to_real max_iterations) (to_real order) 
		(to_real ipol) (to_real wrap), 
			is_image sample && is_image reference
	= error (_ "bad arguments to " ++ "transform_search")
{
	transform_search_op = Operator "transform_search" 
		(transform_search max_error max_iterations order ipol wrap) 
		Operator_type.COMPOUND false;
}

rotate angle image
	= oo_binary_function rotate_op angle image, is_class angle
	= oo_binary'_function rotate_op angle image, is_class image
	= im_similarity image (cos angle) (sin angle) 0 0, 
		is_real angle && is_image image 
	= error (_ "bad arguments to " ++ "rotate")
{
	rotate_op = Operator "rotate" 
		rotate Operator_type.COMPOUND_REWRAP false;
}

conj x
	= oo_unary_function conj_op x, is_class x
	= (re x, -im x), 
		is_complex x || 
		(is_image x && format == Image_format.COMPLEX) ||
		(is_image x && format == Image_format.DPCOMPLEX)
	// assume it's some sort of real 
	= x
{
	format = im_header_int "BandFmt" x;
	conj_op = Operator "conj" conj Operator_type.COMPOUND false;
}

clip2fmt format image
	= oo_unary_function clip2fmt_op image, is_class image
	= im_clip2fmt image (to_real format), is_image image
	= error (_ "bad arguments to " ++ "clip2fmt")
{
	clip2fmt_op = Operator "clip2fmt" 
		(clip2fmt format) Operator_type.COMPOUND_REWRAP false;
}

embed type x y w h im
	= oo_unary_function embed_op im, is_class im
	= im_embed im (to_real type) 
		(to_real x) (to_real y) (to_real w) (to_real h), is_image im
	= error (_ "bad arguments to " ++ "embed")
{
	embed_op = Operator "embed" 
		(embed type x y w h) Operator_type.COMPOUND_REWRAP false;
}

/* Morph a mask with a [[real]] matrix ... turn m2 into an image, morph it 
 * with m1, turn it back to a matrix again.
 */
_morph_2_masks fn m1 m2
	= m''
{
	image = (unsigned char) im_mask2vips (Matrix m2);
	m2_width = get_width image;
	m2_height = get_height image;

	// need to embed m2 in an image large enough for us to be able to 
	// position m1 all around the edges, with a 1 pixel overlap
	image' = embed 0 
		(m1.width / 2) (m1.height / 2) 
		(m2_width + (m1.width - 1)) (m2_height + (m1.height - 1)) 
		image;

	// morph!
	image'' = fn m1 image';

	// back to mask
	m' = im_vips2mask ((double) image'');

	// Turn 0 in output to 128 (don't care).
	m'' 
		= map (map fn) m'.value
	{
		fn a
			= 128, a == 0;
			= a;
	}
}

dilate mask image
	= oo_unary_function dilate_op image, is_class image
	= im_dilate image (to_matrix mask), is_image image
	= error (_ "bad arguments to " ++ "dilate")
{
	dilate_op = Operator "dilate" 
		dilate_object Operator_type.COMPOUND_REWRAP false;
	
	dilate_object x
		= _morph_2_masks dilate mask x, is_matrix x
		= dilate mask x;
}

erode mask image
	= oo_unary_function erode_op image, is_class image
	= im_erode image (to_matrix mask), is_image image
	= error (_ "bad arguments to " ++ "erode")
{
	erode_op = Operator "erode" 
		erode_object Operator_type.COMPOUND_REWRAP false;

	erode_object x
		= _morph_2_masks erode mask x, is_matrix x
		= erode mask x;
}

conv mask image
	= oo_unary_function conv_op image, is_class image
	= im_conv image (to_matrix mask), is_image image
	= error (_ "bad arguments to " ++ "conv")
{
	conv_op = Operator "conv" 
		(conv mask) Operator_type.COMPOUND_REWRAP false;
}

convsep mask image
	= oo_unary_function convsep_op image, is_class image
	= im_convsep image (to_matrix mask), is_image image
	= error (_ "bad arguments to " ++ "convsep")
{
	convsep_op = Operator "convsep" 
		(convsep mask) Operator_type.COMPOUND_REWRAP false;
}

rank w h n image
	= oo_unary_function rank_op image, is_class image
	= im_rank image (to_real w) (to_real h) (to_real n), is_image image
	= error (_ "bad arguments to " ++ "rank")
{
	rank_op = Operator "rank" 
		(rank w h n) Operator_type.COMPOUND_REWRAP false;
}

rank_image n x
	// work for groups too (convenient)
	= rlist x.value, is_Group x
	= rlist x, is_list x
	= error (_ "bad arguments to " ++ "rank_image")
{
	rlist l
		= wrapper ranked, has_wrapper
		= ranked
	{
		has_wrapper = has_member_list (has_member "Image") l;
		wrapper = get_member_list (has_member "Image") (get_member "Image") l;
		ranked = im_rank_image (map get_image l) (to_real n);
	}
}

hist_find image
	= oo_unary_function hist_find_op image, is_class image
	= im_histgr image (-1), is_image image
	= error (_ "bad arguments to " ++ "hist_find")
{
	hist_find_op = Operator "hist_find" 
		hist_find Operator_type.COMPOUND_REWRAP false;
}

hist_find_nD bins image
	= oo_unary_function hist_find_nD_op image, is_class image
	= im_histnD image (to_real bins), is_image image
	= error (_ "bad arguments to " ++ "hist_find_nD")
{
	hist_find_nD_op = Operator "hist_find_nD" 
		(hist_find_nD bins) Operator_type.COMPOUND_REWRAP false;
}

hist_map hist image
	= oo_binary_function hist_map_op hist image, is_class hist
	= oo_binary'_function hist_map_op hist image, is_class image
	= im_maplut image hist, is_image hist && is_image image
	= error (_ "bad arguments to " ++ "hist_map")
{
	hist_map_op = Operator "hist_map" 
		hist_map Operator_type.COMPOUND_REWRAP false;
}

hist_cum hist 
	= oo_unary_function hist_cum_op hist, is_class hist
	= im_histcum hist, is_image hist
	= error (_ "bad arguments to " ++ "hist_cum")
{
	hist_cum_op = Operator "hist_cum" 
		hist_cum Operator_type.COMPOUND_REWRAP false;
}

hist_norm hist 
	= oo_unary_function hist_norm_op hist, is_class hist
	= im_histnorm hist, is_image hist
	= error (_ "bad arguments to " ++ "hist_norm")
{
	hist_norm_op = Operator "hist_norm" 
		hist_norm Operator_type.COMPOUND_REWRAP false;
}

hist_match in ref
	= oo_binary_function hist_match_op in ref, is_class in
	= oo_binary'_function hist_match_op in ref, is_class ref
	= im_histspec in ref, is_image in && is_image ref
	= error (_ "bad arguments to " ++ "hist_match")
{
	hist_match_op = Operator "hist_match" 
		hist_match Operator_type.COMPOUND_REWRAP false;
}

hist_equalize x = hist_map ((hist_norm @ hist_cum @ hist_find) x) x;

hist_equalize_local w h image
	= oo_unary_function hist_equalize_local_op image, is_class image
	= lhisteq image, is_image image
	= error (_ "bad arguments to " ++ "hist_equalize_local")
{
	hist_equalize_local_op = Operator "hist_equalize_local" 
		(hist_equalize_local w h) Operator_type.COMPOUND_REWRAP false;

	// loop over bands, if necessary
	lhisteq im
		= im_lhisteq im (to_real w) (to_real h), get_bands im == 1
		= (foldl1 join @ map lhisteq @ bandsplit) im;
}

// find the threshold below which are percent of the image (percent in [0,1])
// eg. hist_thresh 0.1 x == 12, then x < 12 will light up 10% of the pixels
hist_thresh percent image
	= x
{
	// our own normaliser ... we don't want to norm channels separately
	// norm to [0,1]
	my_hist_norm h = h / max h;
	sum = foldr1 add;

	// normalised cumulative hist
	// we sum the channels before we normalise, because we want to treat them
	// all the same
	h = (my_hist_norm @ sum @ bandsplit @ hist_cum @ hist_find) 
		image.value;

	// threshold that, then use im_profile to search for the x position in the
	// histogram
	x = mean (im_profile (h > percent) 1);
}

resize xfac yfac interp image
	= oo_unary_function resize_op image, is_class image
	= resize_im image, is_image image
	= error (_ "bad arguments to " ++ "resize")
{
	resize_op = Operator "resize" 
		resize_im Operator_type.COMPOUND_REWRAP false;

	xfac' = to_real xfac;
	yfac' = to_real yfac;

	rxfac' = 1 / xfac';
	ryfac' = 1 / yfac';

	resize_im im
		// upscale by integer factor, nearest neighbour
		= im_zoom im xfac' yfac', 
			is_int xfac' && is_int yfac' && 
			xfac' >= 1 && yfac' >= 1 &&
			interp == Interpolate.NEAREST_NEIGHBOUR

		// downscale by integer factor, nearest neighbour
		= im_subsample im rxfac' ryfac', 
			is_int rxfac' && is_int ryfac' && 
			rxfac' >= 1 && ryfac' >= 1 &&
			interp == Interpolate.NEAREST_NEIGHBOUR

		// upscale by any factor, nearest neighbour 
		// can't really do this right ... upscale by integer part, then
		// bilinear to exact size
		= scale (break xfac')?1 (break yfac')?1
			(im_zoom im (break xfac')?0 (break yfac')?0),
			xfac' >= 1 && yfac' >= 1 &&
			interp == Interpolate.NEAREST_NEIGHBOUR

		// downscale by any factor, nearest neighbour 
		// can't really do this right ... downscale by integer part, 
		// then bilinear to exact size
		= scale (1 / (break rxfac')?1) (1 / (break ryfac')?1)
			(im_subsample im (break rxfac')?0 (break ryfac')?0),
			rxfac' >= 1 && ryfac' >= 1 &&
			interp == Interpolate.NEAREST_NEIGHBOUR

		// upscale by any factor, bilinear
		= scale xfac' yfac' im,
			xfac' >= 1 && yfac' >= 1 &&
			interp == Interpolate.BILINEAR

		// downscale by any factor, bilinear
		// block shrink by integer factor, then bilinear resample to 
		// exact
		= scale (1 / (break rxfac')?1) (1 / (break ryfac')?1)
			(im_shrink im (break rxfac')?0 (break ryfac')?0),
			rxfac' >= 1 && ryfac' >= 1 &&
			interp == Interpolate.BILINEAR

		= error ("resize: unimplemented argument combination:\n" ++ 
			"  xfac = " ++ print xfac' ++ "\n" ++
			"  yfac = " ++ print yfac' ++ "\n" ++
			"  interp = " ++ print interp ++ " (" ++
				Interpolate.names.lookup 1 0 interp ++ ")")
	{
		// convert a float scale to integer plus fraction
		// eg. scale by 2.5 becomes [2, 1.25] (x * 2.5 == x * 2 * 1.25)
		break f = [floor f, f / floor f];

		// binlinear resize
		scale xfac yfac im
			= im_affine im 
				xfac 0 0 yfac 
				0 0 
				0 0 
				(rint (get_width im * xfac)) 
				(rint (get_height im * yfac));
	}
}

sharpen radius x1 y2 y3 m1 m2 in
	= oo_unary_function sharpen_op in, is_class in
	= im_sharpen in (to_real radius)
		(to_real x1) (to_real y2) (to_real y3) 
		(to_real m1) (to_real m2), is_image in 
	= error (_ "bad arguments to " ++ "sharpen")
{
	sharpen_op = Operator "sharpen" 
		(sharpen radius x1 y2 y3 m1 m2) 
		Operator_type.COMPOUND_REWRAP false;
}

tone_analyse s m h sa ma ha in
	= oo_unary_function tone_analyse_op in, is_class in
	= im_tone_analyse in
		(to_real s) (to_real m) (to_real h)
		(to_real sa) (to_real ma) (to_real ha), is_image in 
	= error (_ "bad arguments to " ++ "tone_analyse")
{
	tone_analyse_op = Operator "tone_analyse" 
		(tone_analyse s m h sa ma ha)
		Operator_type.COMPOUND_REWRAP false;
}

tone_map hist image
	= oo_binary_function tone_map_op hist image, is_class hist
	= oo_binary'_function tone_map_op hist image, is_class image
	= im_tone_map image hist, is_image hist && is_image image
	= error (_ "bad arguments to " ++ "tone_map")
{
	tone_map_op = Operator "tone_map" 
		tone_map Operator_type.COMPOUND_REWRAP false;
}

tone_build fmt b w s m h sa ma ha 
	= (Image @ clip2fmt fmt) 
		(im_tone_build_range mx mx 
			(to_real b) (to_real w) 
			(to_real s) (to_real m) (to_real h)
			(to_real sa) (to_real ma) (to_real ha))
{
	mx = Image_format.maxval fmt;
}

icc_export depth profile intent in
	= oo_unary_function icc_export_op in, is_class in
	= im_icc_export_depth in
		(to_real depth) (expand profile) (to_real intent), is_image in 
	= error (_ "bad arguments to " ++ "icc_export")
{
	icc_export_op = Operator "icc_export" 
		(icc_export depth profile intent)
		Operator_type.COMPOUND_REWRAP false;
}

icc_import profile intent in
	= oo_unary_function icc_import_op in, is_class in
	= im_icc_import in
		(expand profile) (to_real intent), is_image in 
	= error (_ "bad arguments to " ++ "icc_import")
{
	icc_import_op = Operator "icc_import" 
		(icc_import profile intent)
		Operator_type.COMPOUND_REWRAP false;
}

icc_transform in_profile out_profile intent in
	= oo_unary_function icc_transform_op in, is_class in
	= im_icc_transform in
		(expand in_profile) (expand out_profile) 
		(to_real intent), is_image in 
	= error (_ "bad arguments to " ++ "icc_transform")
{
	icc_transform_op = Operator "icc_transform" 
		(icc_transform in_profile out_profile intent)
		Operator_type.COMPOUND_REWRAP false;
}

icc_ac2rc profile in
	= oo_unary_function icc_ac2rc_op in, is_class in
	= im_icc_ac2rc in (expand profile), is_image in 
	= error (_ "bad arguments to " ++ "icc_ac2rc")
{
	icc_ac2rc_op = Operator "icc_ac2rc" 
		(icc_ac2rc profile)
		Operator_type.COMPOUND_REWRAP false;
}

print_base base in
	= oo_unary_function print_base_op in, is_class in
	= map (print_base base) in, is_list in 
	= print_base_real, is_real in 
	= error (_ "bad arguments to " ++ "print_base")
{
	print_base_op 
		= Operator "print_base" (print_base base) Operator_type.COMPOUND false;

	print_base_real 
		= error "print_base: bad base", base < 2 || base > 16
		= "0", in < 0 || chars == [] 
		= reverse chars
	{
		digits = map (converse remainder base) 
			(takewhile (not_equal 0) 
				(iterate (converse idiv base) in));
		chars = map tohd digits;

		tohd x
			= (char) ((int) '0' + x), x < 10
			= (char) ((int) 'A' + (x - 10));

		idiv a b = (int) (a / b);
	}
}

/* id x: the identity function
 *
 * id :: * -> *
 */
id x = x;

/* const x y: junk y, return x
 * 
 * (const 3) is the function that always returns 3.
 * const :: * -> ** -> *
 */
const x y = x;

/* converse fn a b: swap order of args to fn
 * 
 * converse fn a b == fn b a
 * converse :: (* -> ** -> ***) -> ** -> * -> ***
 */
converse fn a b = fn b a;

/* fix fn x: find the fixed point of a function
 */
fix fn x = limit (iterate fn x);

/* until pred fn n: apply fn to n until pred succeeds; return that value
 *
 * until (more 1000) (multiply 2) 1 = 1024
 * until :: (* -> bool) -> (* -> *) -> * -> *
 */
until pred fn n
	= n, pred n
	= until pred fn (fn n);

/* Infinite list of primes.
 */
primes 
	= 1 : (sieve [2..])
{
	sieve l = hd l : sieve (filter (nmultiple (hd l)) (tl l));
	nmultiple n x = x / n != (int) (x / n);
}

/* Map a 3-ary function over three objects.
 */
map_trinary fn a b c
	= wrap (map3 (map_trinary fn) a' b' c'), 
		is_list a' && is_list b' && is_list c'

	= wrap (map2 (map_trinary fn a') b' c'), 
		is_list b' && is_list c'
	= wrap (map2 (map_trinary (converse31 fn) b') a' c'), 
		is_list a' && is_list c'
	= wrap (map2 (map_trinary (converse32 fn) c') a' b'), 
		is_list a' && is_list b'

	= wrap (map (map_trinary fn a' b') c'), 
		is_list c'
	= wrap (map (map_trinary (converse32 fn) a' c') b'), 
		is_list b'
	= wrap (map (map_trinary (converse34 fn) b' c') a'), 
		is_list a'

	= fn a b c
{
	converse31 fn a b c = fn b a c;
	converse32 fn a b c = fn c a b;
	converse33 fn a b c = fn a c b;
	converse34 fn a b c = fn b c a;

	a'
		= a.value, is_Group a
		= a;
	b'
		= b.value, is_Group b
		= b;
	c'
		= c.value, is_Group c
		= c;
	wrap
		= Group, is_Group a || is_Group b || is_Group c
		= id;
}

/* Map a 2-ary function over a pair of objects.
 */
map_binary fn a b
	= wrap (map2 (map_binary fn) a' b'), is_list a' && is_list b'
	= wrap (map (map_binary fn a') b'), is_list b'
	= wrap (map (map_binary (converse fn) b') a'), is_list a'
	= fn a b
{
	a'
		= a.value, is_Group a
		= a;
	b'
		= b.value, is_Group b
		= b;
	wrap
		= Group, is_Group a || is_Group b 
		= id;
}

/* Map a 1-ary function over an object.
 */
map_unary fn a 
	= wrap (map (map_unary fn) a'), is_list a'
	= fn a
{
	a'
		= a.value, is_Group a
		= a;
	wrap
		= Group, is_Group a
		= id;
}

/* Remove features smaller than x pixels across from an image. This used to be
 * rather complex ... convsep is now good enough to use.
 */
smooth x image = convsep (matrix_gaussian_blur (to_real x * 2)) image;

/* Chop up an image into a list of lists of smaller images. Pad edges with 
 * black.
 */
imagearray_chop tile_width tile_height hoverlap voverlap i
	= map chop' [0, vstep .. height]
{
	width = get_width i;
	height = get_height i;
	bands = get_bands i;
	format = get_format i;
	type = get_type i;

	tile_width' = to_real tile_width;
	tile_height' = to_real tile_height;
	hoverlap' = to_real hoverlap;
	voverlap' = to_real voverlap;

	/* Unique pixels per tile.
	 */
	hstep = tile_width' - hoverlap';
	vstep = tile_height' - voverlap';

	/* Calculate padding ... pad up to tile_size pixel boundary.
	 */
	sx = tile_width' + (width - width % hstep);
	sy = tile_height' + (height - height % vstep);

	/* Expand image with black to pad size.
	 */
	pad = embed 0 0 0 sx sy i;

	/* Chop up a row.
	 */
	chop' y 
		= map chop'' [0, hstep .. width]
	{
		chop'' x = extract_area x y tile_width' tile_height' pad;
	}
}

/* Reassemble image.
 */
imagearray_assemble hoverlap voverlap il
	= (image_set_origin 0 0 @ foldl1 tbj @ map (foldl1 lrj)) il
{
	lrj l r = insert (get_width l + hoverlap) 0 r l;
	tbj t b = insert 0 (get_height t + voverlap) b t;
}

/* Generate an nxn identity matrix.
 */
identity_matrix n
	= error "identity_matrix: n > 0", n < 1
	= map line [0 .. n - 1]
{
	line p = take p [0, 0 ..] ++ [1] ++ take (n - p - 1) [0, 0 ..];
}