#pragma rtGlobals=3		// Use modern global access method and strict wave access.
#pragma IgorVersion = 6.1
#pragma ModuleName = PearlAreaProfiles
#pragma version = 1.05

/// @file
/// @brief profile extraction for multi-dimensional datasets acquired from area detectors.
/// @ingroup ArpesPackage
///
///
///
/// @author matthias muntwiler, matthias.muntwiler@psi.ch
///
/// @copyright 2013-15 Paul Scherrer Institut @n
/// Licensed under the Apache License, Version 2.0 (the "License"); @n
///   you may not use this file except in compliance with the License. @n
///   You may obtain a copy of the License at
///   http://www.apache.org/licenses/LICENSE-2.0

/// @namespace PearlAreaProfiles
/// @brief profile extraction for multi-dimensional datasets acquired from area detectors.
///
/// PearlAnglescanTracker is declared in @ref pearl-area-profiles.ipf.
///

/// 1D cut through 3D dataset, integrate in normal dimensions
///
/// @param dataset
/// @param x1, x2, y1, y2, z1, z2
/// 		coordinates of integration region
///		by default, the coordinates use wave scaling
///		coordinates of rod dimensions (to be preserved) must be nan
///		coordinate pairs don't have to be ordered, i.e. both x1 <= x2 and x1 >= x2 are allowed.
/// @param destname
///		name of destination wave. to be created in current data folder.
///		if empty, the function returns a free wave
/// @param noavg
///		non-zero = calculate the sum, default = 0
///		as of version 1.05, this option should rather be called "calc_sum" or similar.
///		it is noavg for compatibility with older code.
/// @param sdev
///		non-zero = calculate the standard deviation, default = 0
///		by default, the function calculates the average of the integration region
///		set either the noavg or sdev option to select the sum or the standard deviation, respectively.
///		if both options are set, noavg (sum) takes precedence.
/// @param pscale
///		scaling of the slab coordinates x1, x2, ..., z2:
///		zero or default = wave scaling, non-zero = point scaling
///
/// @remark
/// * version 1.02: the specification of the destination coordinates has changed
/// * version 1.04: the function returns an empty wave reference if an error occurred
///
threadsafe function /wave ad_extract_rod(dataset, x1, x2, y1, y2, z1, z2, destname, [noavg, sdev, pscale])
	wave dataset
	variable x1, x2, y1, y2, z1, z2
	string destname
	variable noavg
	variable sdev
	variable pscale

	if (wavedims(dataset) != 3)
		return $""
	endif
	if (ParamIsDefault(noavg))
		noavg = 0
	endif
	if (ParamIsDefault(sdev))
		sdev = 0
	endif
	if (ParamIsDefault(pscale))
		pscale = 0
	endif
	
	variable p1, p2, q1, q2, r1, r2
	if (pscale)
		p1 = x1
		p2 = x2
		q1 = y1
		q2 = y2
		r1 = z1
		r2 = z2
	else
		p1 = round((x1 - DimOffset(dataset, 0)) / DimDelta(dataset, 0))
		p2 = round((x2 - DimOffset(dataset, 0)) / DimDelta(dataset, 0))
		q1 = round((y1 - DimOffset(dataset, 1)) / DimDelta(dataset, 1))
		q2 = round((y2 - DimOffset(dataset, 1)) / DimDelta(dataset, 1))
		r1 = round((z1 - DimOffset(dataset, 2)) / DimDelta(dataset, 2))
		r2 = round((z2 - DimOffset(dataset, 2)) / DimDelta(dataset, 2))
	endif
	
	if ((numtype(p1) == 2) || (numtype(p2) == 2))
		return ad_extract_rod_x(dataset, min(q1, q2), max(q1, q2), min(r1, r2), max(r1, r2), destname, noavg=noavg, sdev=sdev)
	elseif ((numtype(q1) == 2) || (numtype(q2) == 2))
		return ad_extract_rod_y(dataset, min(p1, p2), max(p1, p2), min(r1, r2), max(r1, r2), destname, noavg=noavg, sdev=sdev)
	elseif ((numtype(r1) == 2) || (numtype(r2) == 2))
		return ad_extract_rod_z(dataset, min(p1, p2), max(p1, p2), min(q1, q2), max(q1, q2), destname, noavg=noavg, sdev=sdev)
	else
		return $""
	endif
end

/// 1D cut through 3D dataset along X dimension.
///
/// see ad_extract_rod() for descriptions of common parameters.
threadsafe function /wave ad_extract_rod_x(dataset, q1, q2, r1, r2, destname, [noavg, sdev])
	wave dataset
	variable q1, q2, r1, r2
		// -inf < q1 < q2 < +inf, -inf < r1 < r2 < +inf
	string destname
	variable noavg
	variable sdev
	
	if (ParamIsDefault(noavg))
		noavg = 0
	endif
	if (ParamIsDefault(sdev))
		sdev = 0
	endif
	variable avg = !noavg && !sdev

	q1 = max(q1, 0)
	q2 = min(q2, dimsize(dataset, 1) - 1)
	r1 = max(r1, 0)
	r2 = min(r2, dimsize(dataset, 2) - 1)
	
	if (strlen(destname) > 0)
		duplicate /r=[][q1,q1][r1,r1]/o dataset, $destname
		wave w_dest = $destname
	else
		duplicate /r=[][q1,q1][r1,r1] /free dataset, w_dest
	endif
	redimension /n=(dimsize(w_dest, 0)) w_dest
	setscale /p x dimoffset(dataset, 0), dimdelta(dataset, 0), waveunits(dataset, 0), w_dest
	
	w_dest = 0
	variable qq, rr
	variable nn = 0
	for (qq = q1; qq <= q2; qq += 1)
		for (rr = r1; rr <= r2; rr += 1)
			w_dest += dataset[p][qq][rr]
			nn += 1
		endfor
	endfor
	
	if (sdev)
		duplicate /free w_dest, w_squares
		w_squares = 0
		for (qq = q1; qq <= q2; qq += 1)
			for (rr = r1; rr <= r2; rr += 1)
				w_squares += dataset[p][qq][rr]^2
			endfor
		endfor
	endif
	
	if (avg)
		w_dest /= nn
	elseif (sdev)
		w_dest = sqrt( (w_squares  - w_dest^2 / nn) / (nn - 1) )
	endif
	
	return w_dest
end

/// 1D cut through 3D dataset along Y dimension.
///
/// see ad_extract_rod() for descriptions of common parameters.
threadsafe function /wave ad_extract_rod_y(dataset, p1, p2, r1, r2, destname, [noavg, sdev])
	wave dataset
	variable p1, p2, r1, r2
		// 0 <= p1 < p2 < dimsize(0), 0 <= r1 < r2 < dimsize(2)
	string destname
	variable noavg
	variable sdev
	
	if (ParamIsDefault(noavg))
		noavg = 0
	endif
	if (ParamIsDefault(sdev))
		sdev = 0
	endif
	variable avg = !noavg && !sdev
	
	p1 = max(p1, 0)
	p2 = min(p2, dimsize(dataset, 0) - 1)
	r1 = max(r1, 0)
	r2 = min(r2, dimsize(dataset, 2) - 1)

	if (strlen(destname) > 0)
		duplicate /r=[p1,p1][][r1,r1]/o dataset, $destname
		wave w_dest = $destname
	else
		duplicate /r=[p1,p1][][r1,r1] /free dataset, w_dest
	endif
	redimension /n=(dimsize(w_dest, 1)) w_dest
	setscale /p x dimoffset(dataset, 1), dimdelta(dataset, 1), waveunits(dataset, 1), w_dest
	
	w_dest = 0
	variable pp, rr
	variable nn = 0
	for (pp = p1; pp <= p2; pp += 1)
		for (rr = r1; rr <= r2; rr += 1)
			w_dest += dataset[pp][p][rr]
			nn += 1
		endfor
	endfor

	if (sdev)
		duplicate /free w_dest, w_squares
		w_squares = 0
		for (pp = p1; pp <= p2; pp += 1)
			for (rr = r1; rr <= r2; rr += 1)
				w_squares += dataset[pp][p][rr]^2
				nn += 1
			endfor
		endfor
	endif
	
	if (avg)
		w_dest /= nn
	elseif (sdev)
		w_dest = sqrt( (w_squares  - w_dest^2 / nn) / (nn - 1) )
	endif
	
	return w_dest
end

/// 1D cut through 3D dataset along Z dimension.
///
/// see ad_extract_rod() for descriptions of common parameters.
threadsafe function /wave ad_extract_rod_z(dataset, p1, p2, q1, q2, destname, [noavg, sdev])
	wave dataset
	variable p1, p2, q1, q2
		// 0 <= p1 < p2 < dimsize(0), 0 <= q1 < q2 < dimsize(1)
	string destname
	variable noavg
	variable sdev
	
	if (ParamIsDefault(noavg))
		noavg = 0
	endif
	if (ParamIsDefault(sdev))
		sdev = 0
	endif
	variable avg = !noavg && !sdev

	p1 = max(p1, 0)
	p2 = min(p2, dimsize(dataset, 0) - 1)
	q1 = max(q1, 0)
	q2 = min(q2, dimsize(dataset, 1) - 1)

	if (strlen(destname) > 0)
		duplicate /r=[p1,p1][q1,q1][]/o dataset, $destname
		wave w_dest = $destname
	else
		duplicate /r=[p1,p1][q1,q1][] /free dataset, w_dest
	endif
	redimension /n=(dimsize(w_dest, 2)) w_dest
	setscale /p x dimoffset(dataset, 2), dimdelta(dataset, 2), waveunits(dataset, 2), w_dest
	
	w_dest = 0
	variable pp, qq
	variable nn = 0
	for (pp = p1; pp <= p2; pp += 1)
		for (qq = q1; qq <= q2; qq += 1)
			w_dest += dataset[pp][qq][p]
			nn += 1
		endfor
	endfor

	if (sdev)
		duplicate /free w_dest, w_squares
		w_squares = 0
		for (pp = p1; pp <= p2; pp += 1)
			for (qq = q1; qq <= q2; qq += 1)
				w_squares += dataset[pp][qq][p]^2
				nn += 1
			endfor
		endfor
	endif
	
	if (avg)
		w_dest /= nn
	elseif (sdev)
		w_dest = sqrt( (w_squares  - w_dest^2 / nn) / (nn - 1) )
	endif
	
	return w_dest
end

/// 2D cut through 3D dataset, integrate in normal dimension
///
/// @param	dataset
/// @param	x1, x2, y1, y2, z1, z2
///		coordinates of integration region.
///		by default, the coordinates use wave scaling.
///		coordinates of slab dimensions (to be preserved) must be nan.
///		coordinate pairs don't have to be ordered, i.e. both x1 <= x2 and x1 >= x2 are allowed.
///		coordinates can be out of range (-inf and +inf allowed) to select the whole range.
/// @param	destname
///		name of destination wave. to be created in current data folder.
///		if empty, the function returns a free wave.
/// @param	noavg
///		zero or default = average, non-zero = sum.
/// @param	pscale
///		scaling of the slab coordinates x1, x2, ..., z2:
///		zero or default = wave scaling, non-zero = point scaling.
///
/// @remark
/// * version 1.02: the specification of the destination coordinates has changed
/// * version 1.04: the function returns an empty wave reference if an error occurred
///
threadsafe function /wave ad_extract_slab(dataset, x1, x2, y1, y2, z1, z2, destname, [noavg, pscale])
	wave dataset
	variable x1, x2, y1, y2, z1, z2
	string destname
	variable noavg
	variable pscale
	
	if (wavedims(dataset) != 3)
		return $""
	endif
	if (ParamIsDefault(noavg))
		noavg = 0
	endif
	if (ParamIsDefault(pscale))
		pscale = 0
	endif
	
	variable p1, p2, q1, q2, r1, r2
	if (pscale)
		p1 = x1
		p2 = x2
		q1 = y1
		q2 = y2
		r1 = z1
		r2 = z2
	else
		p1 = round((x1 - DimOffset(dataset, 0)) / DimDelta(dataset, 0))
		p2 = round((x2 - DimOffset(dataset, 0)) / DimDelta(dataset, 0))
		q1 = round((y1 - DimOffset(dataset, 1)) / DimDelta(dataset, 1))
		q2 = round((y2 - DimOffset(dataset, 1)) / DimDelta(dataset, 1))
		r1 = round((z1 - DimOffset(dataset, 2)) / DimDelta(dataset, 2))
		r2 = round((z2 - DimOffset(dataset, 2)) / DimDelta(dataset, 2))
	endif
	
	if ((numtype(p1) < 2) && (numtype(p2) < 2))
		return ad_extract_slab_x(dataset, min(p1, p2), max(p1, p2), destname, noavg=noavg)
	elseif ((numtype(q1) < 2) && (numtype(q2) < 2))
		return ad_extract_slab_y(dataset, min(q1, q2), max(q1, q2), destname, noavg=noavg)
	elseif ((numtype(r1) < 2) && (numtype(r2) < 2))
		return ad_extract_slab_z(dataset, min(r1, r2), max(r1, r2), destname, noavg=noavg)
	else
		return $""
	endif
end

threadsafe function /wave ad_extract_slab_x(dataset, p1, p2, destname, [noavg])
	wave dataset
	variable p1, p2
		// x coordinate range (point scaling) to be integrated
		// -inf <= p1 < p2 <= +inf
	string destname // name of destination wave. to be created in current data folder. overrides existing.
		// if empty, the function returns a free wave
	variable noavg // zero or default = average, non-zero = sum
	
	if (ParamIsDefault(noavg))
		noavg = 0
	endif
	p1 = max(p1, 0)
	p2 = min(p2, dimsize(dataset, 0) - 1)

	if (strlen(destname) > 0)
		duplicate /r=[p1,p1][][]/o dataset, $destname
		wave w_dest = $destname
	else
		duplicate /r=[p1,p1][][] /free dataset, w_dest
	endif
	redimension /n=(dimsize(w_dest, 1), dimsize(w_dest, 2)) w_dest
	setscale /p x dimoffset(dataset, 1), dimdelta(dataset, 1), waveunits(dataset, 1), w_dest
	setscale /p y dimoffset(dataset, 2), dimdelta(dataset, 2), waveunits(dataset, 2), w_dest
	
	w_dest = 0
	variable pp
	variable nn = 0
	for (pp = p1; pp <= p2; pp += 1)
		w_dest += dataset[pp][p][q]
		nn += 1
	endfor

	if (noavg == 0)
		w_dest /= nn
	endif
	
	return w_dest
end

threadsafe function /wave ad_extract_slab_y(dataset, q1, q2, destname, [noavg])
	wave dataset
	variable q1, q2
		// y coordinate range (point scaling) to be integrated
		// -inf <= q1 < q2 <= +inf
	string destname // name of destination wave. to be created in current data folder. overrides existing.
		// if empty, the function returns a free wave
	variable noavg // zero or default = average, non-zero = sum
	
	if (ParamIsDefault(noavg))
		noavg = 0
	endif
	q1 = max(q1, 0)
	q2 = min(q2, dimsize(dataset, 1) - 1)

	if (strlen(destname) > 0)
		duplicate /r=[][q1,q1][]/o dataset, $destname
		wave w_dest = $destname
	else
		duplicate /r=[][q1,q1][] /free dataset, w_dest
	endif
	redimension /n=(dimsize(w_dest, 0), dimsize(w_dest, 2)) w_dest
	setscale /p x dimoffset(dataset, 0), dimdelta(dataset, 0), waveunits(dataset, 0), w_dest
	setscale /p y dimoffset(dataset, 2), dimdelta(dataset, 2), waveunits(dataset, 2), w_dest
	
	w_dest = 0
	variable qq
	variable nn = 0
	for (qq = q1; qq <= q2; qq += 1)
		w_dest += dataset[p][qq][q]
		nn += 1
	endfor

	if (noavg == 0)
		w_dest /= nn
	endif
	
	return w_dest
end

threadsafe function /wave ad_extract_slab_z(dataset, r1, r2, destname, [noavg])
	wave dataset
	variable r1, r2
		// z coordinate range (point scaling) to be integrated
		// -inf <= r1 < r2 <= +inf
	string destname // name of destination wave. to be created in current data folder. overrides existing.
		// if empty, the function returns a free wave
	variable noavg // zero or default = average, non-zero = sum
	
	if (ParamIsDefault(noavg))
		noavg = 0
	endif
	r1 = max(r1, 0)
	r2 = min(r2, dimsize(dataset, 2) - 1)

	if (strlen(destname) > 0)
		duplicate /r=[][][r1,r1]/o dataset, $destname
		wave w_dest = $destname
	else
		duplicate /r=[][][r1,r1] /free dataset, w_dest
	endif
	redimension /n=(dimsize(w_dest, 0), dimsize(w_dest, 1)) w_dest
	setscale /p x dimoffset(dataset, 0), dimdelta(dataset, 0), waveunits(dataset, 0), w_dest
	setscale /p y dimoffset(dataset, 1), dimdelta(dataset, 1), waveunits(dataset, 1), w_dest
	
	w_dest = 0
	variable rr
	variable nn = 0
	for (rr = r1; rr <= r2; rr += 1)
		w_dest += dataset[p][q][rr]
		nn += 1
	endfor

	if (noavg == 0)
		w_dest /= nn
	endif
	
	return w_dest
end

/// 1D cut through 2D dataset along X dimension, new destination wave.
///
threadsafe function /wave ad_profile_x(dataset, q1, q2, destname, [noavg])
	wave dataset
	variable q1, q2 // -inf <= q1 < q2 <= +inf
		// deprecated: q2 = -1 stands for dimsize(0) - 1
	string destname // name of destination wave. to be created in current data folder. overrides existing.
		// if empty, the function returns a free wave
	variable noavg // zero or default = average, non-zero = sum
	
	if (ParamIsDefault(noavg))
		noavg = 0
	endif

	if (strlen(destname) > 0)
		duplicate /r=[0,0][] /o dataset, $destname
		wave w_dest = $destname
	else
		duplicate /r=[0,0][] /free dataset, w_dest
	endif
	
	return ad_profile_x_w(dataset, q1, q2, w_dest, noavg=noavg)
end

/// 1D cut through 2D dataset along X dimension, existing destination wave.
///
threadsafe function /wave ad_profile_x_w(dataset, q1, q2, destwave, [noavg])
	wave dataset
	variable q1, q2 // -inf <= q1 < q2 <= +inf
		// deprecated: q2 = -1 stands for dimsize(0) - 1
	wave destwave // existing destination wave
	variable noavg // zero or default = average, non-zero = sum
	
	if (ParamIsDefault(noavg))
		noavg = 0
	endif

	redimension /n=(dimsize(dataset, 0)) destwave
	setscale /p x dimoffset(dataset, 0), dimdelta(dataset, 0), waveunits(dataset, 0), destwave
	setscale d 0, 0, waveunits(dataset, -1), destwave
	
	q1 = max(q1, 0)
	if (q2 < 0)
		q2 = inf
	endif
	q2 = min(q2, dimsize(dataset, 1) - 1)
	
	destwave = 0
	variable qq
	variable nn = 0
	for (qq = q1; qq <= q2; qq += 1)
		destwave += dataset[p][qq]
		nn += 1
	endfor

	if (noavg == 0)
		destwave /= nn
	endif
	
	return destwave
end

/// 1D cut through 2D dataset along Y dimension, new destination wave.
///
threadsafe function /wave ad_profile_y(dataset, p1, p2, destname, [noavg])
	wave dataset
	variable p1, p2 // -inf <= p1 < p2 < inf
		// deprecated: p2 = -1 stands for dimsize(0) - 1
	string destname // name of destination wave. to be created in current data folder. overrides existing.
		// if empty, the function returns a free wave
	variable noavg // zero or default = average, non-zero = sum
	
	if (ParamIsDefault(noavg))
		noavg = 0
	endif

	if (strlen(destname) > 0)
		duplicate /r=[][0,0] /o dataset, $destname
		wave w_dest = $destname
	else
		duplicate /r=[][0,0] /free dataset, w_dest
	endif
	MatrixTranspose w_dest
	
	return ad_profile_y_w(dataset, p1, p2, w_dest, noavg=noavg)
end

/// 1D cut through 2D dataset along X dimension, existing destination wave.
///
threadsafe function /wave ad_profile_y_w(dataset, p1, p2, destwave, [noavg])
	wave dataset
	variable p1, p2 // -inf <= p1 < p2 < inf
		// deprecated: p2 = -1 stands for dimsize(0) - 1
	wave destwave // existing destination wave
	variable noavg // zero or default = average, non-zero = sum
	
	if (ParamIsDefault(noavg))
		noavg = 0
	endif

	redimension /n=(dimsize(dataset, 1)) destwave
	setscale /p x dimoffset(dataset, 1), dimdelta(dataset, 1), waveunits(dataset, 1), destwave
	setscale d 0, 0, waveunits(dataset, -1), destwave
	
	p1 = max(p1, 0)
	if (p2 < 0)
		p2 = inf
	endif
	p2 = min(p2, dimsize(dataset, 0) - 1)
	
	destwave = 0
	variable pp
	variable nn = 0
	for (pp = p1; pp <= p2; pp += 1)
		destwave += dataset[pp][p]
		nn += 1
	endfor

	if (noavg == 0)
		destwave /= nn
	endif
	
	return destwave
end

threadsafe function calc_y_profile_mins(image)
	// experimental
	wave image
	
	wave yminlocs = ad_profile_x(image, 0, 0, "ymins", noavg=1)
	variable nx = dimsize(image, 0)
	variable ix
	for (ix = 0; ix < nx; ix += 1)
		wave profile = ad_profile_y(image, ix, ix, "", noavg=1)
		wavestats /q/m=1 profile
		yminlocs[ix] = v_minloc
	endfor
end

/// collect profiles from a multi-scan.
///
/// @warning experimental: name and interface of this function may change.
///
function ad_collect_multiscan_y(dataset, positions, destwave, [noavg])
	wave dataset
	wave positions
	wave destwave
	variable noavg
	
	variable tol = (wavemax(positions) - wavemin(positions)) / numpnts(positions) / 100
	
	duplicate /free positions, positions_sorted
	sort positions_sorted, positions_sorted
	duplicate /free positions_sorted, positions_diff
	differentiate /p /meth=2 positions_sorted /d=positions_diff
	positions_diff[0] = 1
	extract /free positions_sorted, positions_unique, positions_diff > tol
	variable n_unique = numpnts(positions_unique)
	redimension /n=(dimsize(dataset, 0), n_unique) destwave
	
	variable i
	variable nx, ny
	for (i = 0; i < n_unique; i += 1)
		extract /free dataset, data_extract, abs(positions[q] - positions_unique[i]) < tol
		nx = dimsize(dataset, 0)
		ny = dimsize(data_extract, 0) / nx
		redimension /n=(nx, ny) data_extract
		wave profile = ad_profile_x(data_extract, -inf, inf, "", noavg=noavg)
		destwave[][i] = profile[p]
	endfor
end