#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3		// Use modern global access method and strict wave access.

function polyArea (polyx, polyy)		// differs from Igor function, because it depends on handedness
	wave polyx, polyy
	variable mysum=0, size, i
	
	duplicate /free polyx tempx
	duplicate /free polyy tempy
	size=numpnts (tempx)
	for (i = 0; i < size - 1; i += 1)
		mysum += (tempx[i] + tempx[i + 1]) * (tempy[i] - tempy[i + 1])
	endfor
	mysum += (tempx[i] + tempx[0]) * (tempy[i] - tempy[0])
	return mysum/2
end

static function getcrossings (w_inpoly, wcrossin, wcrossout)
	wave w_inpoly, wcrossin, wcrossout
	variable i, oldpts
	
	redimension /n=0 wcrossin, wcrossout
	for (i = 1; i < numpnts (w_inpoly); i += 1)
		if (w_inpoly[i - 1] == 1 && w_inpoly[i] == 0)
			oldpts = numpnts (wcrossout)
			redimension /n=(oldpts + 1) wcrossout
			wcrossout[oldpts] = i
		elseif (w_inpoly[i - 1] == 0 && w_inpoly[i] == 1)
			oldpts = numpnts (wcrossin)
			redimension /n=(oldpts + 1) wcrossin
			wcrossin[oldpts] = i - 1
		endif
	endfor
	if (w_inpoly[i - 1] == 1 && w_inpoly[0] == 0)
		oldpts = numpnts (wcrossout)
		redimension /n=(oldpts + 1) wcrossout
		wcrossout[oldpts] = 0
	elseif (w_inpoly[i - 1] == 0 && w_inpoly[0] == 1)
		oldpts = numpnts (wcrossin)
		redimension /n=(oldpts + 1) wcrossin
		wcrossin[oldpts] = i - 1
	endif
	
	if (wcrossout[0] > wcrossin[0] && numpnts (wcrossout) > 1)	// must have started in, so shift to delimit parts that are out
		i = wcrossout[0]
		wcrossout[0, oldpts - 1] = wcrossout[p + 1]
		wcrossout[oldpts] = i
	endif
end
static function copysegment (wfromx, wfromy, fromstart, fromend, wtox, wtoy, tostart)
	wave wfromx, wfromy, wtox, wtoy
	variable fromstart, fromend, tostart
	if (fromstart < fromend)	// i.e. segment well-behaved, and doesn't wrap around to beginning
		wtox[tostart, tostart + fromend - fromstart] = wfromx[p - tostart + fromstart]
		wtoy[tostart, tostart + fromend - fromstart] = wfromy[p - tostart + fromstart]
		tostart += fromend - fromstart + 1
	else	// segment wraps around to beginning of wave, so need to copy in two parts
		wtox[tostart, tostart + numpnts (wfromx) - fromstart - 1] =  wfromx[p - tostart + fromstart]
		wtoy[tostart, tostart + numpnts (wfromx) - fromstart - 1] =  wfromy[p - tostart + fromstart]
		tostart += numpnts (wfromx) - fromstart
		wtox[tostart, tostart + fromend] = wfromx[p - tostart]
		wtoy[tostart, tostart + fromend] = wfromy[p - tostart]
		tostart += fromend + 1
	endif
	return (tostart)
end

function polyUnion (wx1, wy1, wx2, wy2, destbase)
	// assumes polys are well-behaved, i.e. no twisting
	wave wx1, wy1, wx2, wy2
	string destbase
	variable index, i, i2
	
	//make clockwise
	duplicate /free wx1 tempx1
	duplicate /free wy1 tempy1
	if (polyarea (wx1, wy1) < 0)
		reverse tempx1; reverse tempy1
	endif
	duplicate /free wx2 tempx2
	duplicate /free wy2 tempy2
	if (polyarea (wx2, wy2) < 0)
		reverse tempx2; reverse tempy2
	endif
	if (tempx1[0] == tempx1[numpnts (tempx1) - 1] && tempy1[0] == tempy1[numpnts (tempy1) - 1])	// strip off extra point that closes loop?
		DeletePoints numpnts(tempx1) - 1, 1, tempx1, tempy1
	endif 
	if (tempx2[0] == tempx2[numpnts (tempx2) - 1] && tempy2[0] == tempy2[numpnts (tempy2) - 1])	// strip off extra point that closes loop?
		DeletePoints numpnts(tempx2) - 1, 1, tempx2, tempy2
	endif 
	
	// let Igor figure out overlap
	FindPointsInPoly tempx1, tempy1, tempx2, tempy2
	wave w_InPoly
	duplicate /free w_inPoly w1in2
	FindPointsInPoly tempx2, tempy2, tempx1, tempy1
	duplicate /free w_inPoly w2in1
	
	if (sum (w1in2) == 0 && sum (w2in1) == 0)
		Doalert /T="Warning" 0, "No overlap, so merge not possible"
		return 0
	elseif (sum (w1in2) == numpnts (w1in2))	// 1 fully contained in 2, so union is 2
		duplicate /o wx2 $(destbase + "x")
		duplicate /o wy2 $(destbase + "y")
		return 1
	elseif (sum (w2in1) == numpnts (w2in1))	// 2 fully contained in 1, so union is 1
		duplicate /o wx1 $(destbase + "x")
		duplicate /o wy1 $(destbase + "y")
		return 1
	endif
	
	// proceed with partial overlap
	make /o/n=(numpnts (tempx1) + numpnts (tempx2) - sum (w2in1) - sum (w1in2) + 1) $(destbase + "x"), $(destbase + "y")
	wave wdx = $(destbase + "x"), wdy = $(destbase + "y")
	// first find crossing points (from inside to outside) on poly1 and poly2
	make /n=0/free tempcrossin1, tempcrossin2, tempcrossout1, tempcrossout2
	getcrossings (w1in2, tempcrossin1, tempcrossout1)
	getcrossings (w2in1, tempcrossin2, tempcrossout2)
	
	// find segment on poly2 that picks up where first segment on poly1 goes inside poly2
	duplicate /free tempcrossin2 tempdist	
	tempdist = (tempx1[tempcrossin1[0]] - tempx2[tempcrossout2[p]])^2 + (tempy1[tempcrossin1[0]] - tempy2[tempcrossout2[p]])^2
	wavestats /M=1/Q tempdist
	i2 = v_minloc	// locate closest segment on poly2
	index = 0
	for (i = 0; i < numpnts (tempcrossout1); i += 1)		// go through all segments
		index = copysegment (tempx1, tempy1, tempcrossout1[i], tempcrossin1[i], wdx, wdy, index) // copy poly1 segment
		index = copysegment (tempx2, tempy2, tempcrossout2[i2], tempcrossin2[i2], wdx, wdy, index) // copy nearest poly2 segment
		i2 = (i2 > numpnts (tempcrossout2) - 2 ? 0 : i2 + 1)	// increment poly2 segment counter (and wrap to 0 if needed)
	endfor
	return 1
end