*! clippolyline v1.1 Asjad Naqvi 31Jul2022
* Clip polyline vectors based on a custom bounding box


cap program drop clipline


***********************
// 	  clippolyline   //
***********************

program define clippolyline, eclass
	version 15
	
	syntax namelist(max=1), Box(numlist min=4 max=4)
	
	*** checks

	// check gtools
	capture findfile gtools.ado
	if _rc != 0 {
		display as error "gtools package is missing. Click here to install: {stata ssc install clipline, replace}"
		exit
	}		
	
	
	
	capture confirm file "`namelist'.dta"
	if _rc {
	   di as err "{p}File {bf:`namelist'} not found{p_end}"
	   exit 601
	}	
		

	// push the bounds to locals	
	tokenize `box'		
		local xmin = `1'	
		local xmax = `2'
		local ymin = `3'
		local ymax = `4'
	
		if `xmax' < `xmin' {
			di as error "xmax < xmin. Please make sure the coordinates are specified in the correct order."
			exit
		}
		
		if `ymax' < `ymin' {
			di as error "ymax < ymin. Please make sure the coordinates are specified in the correct order."
			exit
		}
				
		
		

// main routine	
	
qui {
	preserve	
		
		use "`namelist'.dta", clear
		
		count if shape_order > 3
		
		
		if r(N) > 0 {
			noisily di in red "`namelist' is not a valid polyline file. It contains polygons. Try {bf:clipolygon} instead."
			exit
		}

	
		// check if the box is valid (v1.1)
		
		cap drop box
		gen  box = 1 if inrange(_X, `xmin', `xmax') & inrange(_Y, `ymin', `ymax')


		bysort _ID: egen double boxcheck = mean(box)
		drop if boxcheck==.
		drop boxcheck	
		drop box
		sort _ID shape_order		
		
		count 
		
		if `r(N)' == 0 {
			di as error "No shapes are contained within the box. Please check the box bounds."
			exit
		}	
	
	
	
		keep _ID _X _Y shape_order		
		drop if shape_order==1

		// convert to line coordinates
		greshape wide _X _Y, i(_ID) j(shape_order)


		// call clipline program
		local box2: subinstr local box " " ", ", all
		clipline _X2 _Y2 _X3 _Y3, box(`box2')


		cap drop clip*

		gen clip_X2 = .
		gen clip_Y2 = .

		gen clip_X3 = .
		gen clip_Y3 = .

		mata: clipdata = .
		mata: st_view(clipdata, .,"clip_X2 clip_Y2 clip_X3 clip_Y3")
		mata: cliplist2 = ((1::rows(cliplist)), cliplist)
		mata: clipdata[cliplist2[.,1], .] = cliplist2[.,2::5]

		drop if clip_X2 ==.
		drop _X2 _Y2 _X3 _Y3

		gen clip_X1 = .
		gen clip_Y1 = .


		greshape long clip_X clip_Y, i(_ID) j(shape_order)

		ren clip* *
		order _ID _X _Y shape_order

		sort _ID shape_order

		compress
		save "`namelist'_clipped.dta", replace

	restore
}	

	di in green "Done! File exported as `namelist'_clipped.dta"

	
end	





cap program drop clipline

*********************
// 	   clipline    //
*********************

program define clipline, sortpreserve
	version 15
	syntax varlist(min=4 max=4 numeric), 	///
		[Box(numlist min=4 max=4)] 			///	  // manually define the bounding box 	
		[OFFset(real 0)] 					///   // manually define the offset of the box
		[lines] 							///   // export lines inside the box
		[addbox]								  // export bounding box coordinates	
	
	di in green "clipline: Initializing"
	
	// get the variables
	tokenize "`varlist'"
	local x1 `1'
	local y1 `2'
	local x2 `3'
	local y2 `4'
	
	
	mata: points = st_data(.,("`x1'", "`y1'", "`x2'", "`y2'"))

	
	mata: xmin = .
	mata: xmax = .
	mata: ymin = .
	mata: ymax = .	
	
	
	// bounding box
	//di in green "clipline: Bounding box"
	
	if "`box'" != "" {
        local box2: subinstr local box " " ", ", all
		mata: getmanbounds(points,`box2',xmin,xmax,ymin,ymax,`offset')	
	}
	else {
		mata: getautobounds(points,xmin,xmax,ymin,ymax,`offset') 
	}
	
	
	//mata: xmin,xmax,ymin,ymax

	
	mata: cliplist = clipline_core(points, xmin, xmax, ymin, ymax)
	di in green "clipline: Finished"
		
	//cliplist	
	
	if "`lines'" != "" {
		mata: st_local("newobs", strofreal(rows(cliplist)))

		if `newobs' > _N {
			set obs `newobs'
		}
		
		cap drop clip*	
		mata: st_matrix("clip",cliplist)
		mat colnames clip = "clip_x1" "clip_y1" "clip_x2" "clip_y2"

		svmat clip, n(col)
		
		di in green "clipline: Clipped lines exported back to Stata"
	}
	
	if "`addbox'" != "" {
		cap drop box*

		qui {
			gen box_x = .
			gen box_y = .
			
			replace box_x = xmin in 1
			replace box_y = ymin in 1

			replace box_x = xmin in 2
			replace box_y = ymax in 2
			
			replace box_x = xmax in 3
			replace box_y = ymax in 3
			
			replace box_x = xmax in 4
			replace box_y = ymin in 4
			
			replace box_x = xmin in 5
			replace box_y = ymin in 5		
		}

		di in green "clipline: Box exported back to Stata"
	}
	
	
end


************************
// 	  clipline_core	  //		
************************


cap mata: mata drop clipline_core()

mata // clipline_core
function clipline_core(points, xmin, xmax, ymin, ymax)
{
	cliplist = J(rows(points),4,.)
	
	for (i=1; i <= rows(points); i++) {			
		cliplist[i,.] = cliproutine(points[i,1], points[i,2], points[i,3], points[i,4], xmin, xmax, ymin, ymax)	
		}
			
	return(cliplist)
}
end


********************
//  getmanbounds  //		
********************

cap mata: mata drop getmanbounds()

mata // getmanbounds
void getmanbounds(points, xlo, xhi, ylo, yhi, xmin, xmax, ymin, ymax, real scalar value)
{
	
	
	
	displacex = abs((xhi - xlo) * value)
	displacey = abs((yhi - ylo) * value)

	xmin 	  = xlo + displacex
	xmax 	  = xhi - displacex

	ymin 	  = ylo + displacey
	ymax 	  = yhi - displacey

	
	// pass them to stata scalars for drawing later
	st_numscalar("xmin", xmin)
	st_numscalar("xmax", xmax)

	st_numscalar("ymin", ymin)
	st_numscalar("ymax", ymax)

}

end



************************
//  getautobounds	  //		
************************

cap mata: mata drop getautobounds()

mata // getautobounds
void getautobounds(points,xmin,xmax,ymin,ymax,real scalar value)
{
	real matrix mymat
	mymat = colshape(points,2)
	
	displacex = abs((max(mymat[.,1]) - min(mymat[.,1])) * value)
	displacey = abs((max(mymat[.,2]) - min(mymat[.,2])) * value)

	xmin 	  = min(mymat[.,1]) + displacex
	xmax 	  = max(mymat[.,1]) - displacex

	ymin 	  = min(mymat[.,2]) + displacey
	ymax 	  = max(mymat[.,2]) - displacey

	
	// pass them to stata scalars for drawing later
	st_numscalar("xmin", xmin)
	st_numscalar("xmax", xmax)

	st_numscalar("ymin", ymin)
	st_numscalar("ymax", ymax)
}
end

********************
//  cliproutine	  //		
********************

cap mata: mata drop cliproutine()
mata:  // cliproutine
	function cliproutine(x1, y1, x2, y2, minX, maxX, minY, maxY)  
	{
		real scalar code1, code2, accept, x, y

		// Defining region codes 
		LEFT   = 1  // 0001 
		RIGHT  = 2  // 0010 
		BOTTOM = 4  // 0100 
		TOP    = 8  // 1000 
		
		
		code1 = computeCode(x1, y1, minX, maxX, minY, maxY)
		code2 = computeCode(x2, y2, minX, maxX, minY, maxY)
		
		accept = 0
		t = 0 // counter
		
		zz = 0
		while (zz == 0)  {
		
			// If both endpoints lie within rectangle 
			if (code1 == 0 & code2 == 0) {
				accept = 1
				break
			}
			// If both endpoints are outside rectangle 
			else if (code1==code2 & code1!=0 & code2!=0) {
				accept = 0
				
				x1 = .
				y1 = .
				x2 = .
				y2 = .	
				break
			}
			
			// Some segment lies within the rectangle 			
			else {	
				x = 1
				y = 1
				
				if (code1 != 0) {
					codeout = code1
				}
				else {
					codeout = code2
				}
				
				// point is above the clip rectangle 
				if      (codeout >= TOP) {   		
					x = x1 + ((x2 - x1) * (maxY - y1) / (y2 - y1))
					y = maxY
				}
				
				// point is below the clip rectangle
				else if (codeout >= BOTTOM) {
					x = x1 + ((x2 - x1) * (minY - y1) / (y2 - y1))
					y = minY
				}
				
				// point is to the right of the clip rectangle 
				else if (codeout >= RIGHT) {
					y = y1 + ((y2 - y1) * (maxX - x1) / (x2 - x1))
					x = maxX
				}
				
				// point is to the left of the clip rectangle 
				else if (codeout >= LEFT) {
					y = y1 + ((y2 - y1) * (minX - x1) / (x2 - x1))
					x = minX
				}	
				
			if (codeout == code1) { 
                x1 = x 
                y1 = y 
                code1 = computeCode(x1, y1, minX, maxX, minY, maxY) 
			}
			else {
                x2 = x 
                y2 = y 
                code2 = computeCode(x2, y2, minX, maxX, minY, maxY) 
			}
		
		t = t + 1
		if (t > 100) {
		//	printf("bad egg\n")
				x1 = .
				y1 = .
				x2 = .
				y2 = .
			break
			}
		
		}

	}
			
	return(x1,y1,x2,y2)
}
end




*********************
// 	 computeCode   //
*********************

cap mata: mata drop computeCode()
mata // computeCode
	function computeCode(xx, yy, minX, maxX, minY, maxY)  
	{
		real scalar code
		code = 0   // 1 = left, 2 = right, 4 = bottom, 8 = top (binary)
			if (xx < minX) code = code + 1
			if (xx > maxX) code = code + 2		
			if (yy < minY) code = code + 4
			if (yy > maxY) code = code + 8	
		return(code)
	}
end



**** END OF CLIPLINE ****