## -*- SM -*-
##
## Gabor Toth writes:
## 
## I wrote a couple of macros for plotting data defined on 2D irregular grids.
## The names are taken from the analogous IDL macro names, but the algorithm
## of propagating neighborhood and distance info was designed by me.
## 
## REGULARGRID is used to create a regular grid which covers the irregular points.
## 
## TRIANGULATE determines the 3 closest irregular points to each regular points, 
## and 3 weights proportional to the inverse of the distance to the irregular 
## points. The distances are measured in norm 1 now, but that or the weighting 
## scheme can easily be changed. 
## 
## Once the TRIANGULATE macro is done, one can interpolate any variable defined 
## on the irregular grid to the regular grid by a single command:
## 
## areg=TRIGRID(airr)
## 
## This is a single line macro, thus very fast.
## 
## You can try all this by running TRIGRIDTRY. To get a pretty picture, 
## run TRISHOW afterwards.
## 
## I hope you'll find these macros useful enough to put them into the 
## SM macro library.
## 
## Best regards,
## 
## Gabor
## 
## ============================================================================
trigridtry # Demonstrates the use of trigrid
		set xirreg={0   3   4   5   0}
		set yirreg={5   4   3   0   0}
		set zirreg={3.5 4.1 4.1 3.5 2}
		irregular_image xirreg yirreg zirreg 10
		define Min local define Max local
		minmax Min Max
		ERASE
		surface 3 $Min $Max
		#
irregular_image 44 # Given 3 vectors $1, $2, and $3, giving x, y, and z values,
		# construct an SM image of size $4x$4
		define v local
		foreach v (Aix Aiy Ax Ay nei1 nei2 nei3 x y wei1 wei2 wei3) {
		   SET $v LOCAL
		}

		triangulate $1 $2  $4 $4
		set zreg local
		set zreg=trigrid($3)
		# Visualize zreg
		image($(nx),$(ny)) $(x[0]) $(x[dimen(x)-1]) $(y[0]) $(y[dimen(y)-1])
		set image[Aix,Aiy]=zreg
		#
trishow	22	## Visualize neighbor info nei1,nei2,nei3
		# for the x,y points $1, $2
	# The irregular gridpoints are the biggest squares of various colors
	# For all regular points three squares are shown with colors 
	# corresponding to the 1st, 2nd and 3rd closest irregular gridpoints
		limits $1 $2
		define ctype | define expand | define ptype |
		ptype 4 0
		do i=0,dimen($1)-1{
		   expand 13
		   #ptype $($i+3) 0
		   ctype $($i+3)
		   rel $1[$i] $2[$i] dot
		   do n=1,3{
		      expand 12-3*$n
		      points Ax Ay if(nei$n==$i)
		   }
		}
		expand $expand ptype $ptype ctype $ctype
regulargrid 6 # Determine x and y coordinates for an nx*ny enclosing grid 
		# covering all the irregular points at positions xirr yirr.
		# Usage: regulargrid xirr yirr nx ny x y
		define Min local define Max local
		vecminmax $1 Min Max
		if($Min == $Max) {
		   USER ABORT "Error: min == max for vector $!1"
		}
		set $5=$Min,$Max,($Max-$Min)/($3-1)

		vecminmax $2 Min Max
		if($Min == $Max) {
		   USER ABORT "Error: min == max for vector $!2"
		}
		set $6=$Min,$Max,($Max-$Min)/($4-1)
		#
trigrid	1	# Interpolates an array on the irregular grid to the regular grid
		# Usage: Areg=trigrid(Airr)
		if(!is_vector(nei1)) {
		   USER ABORT "You must use triangulate before trigrid"
		}
		set $0=$1[nei1]*wei1+$1[nei2]*wei2+$1[nei3]*wei3
triangulate 24# Determine regular grid & interpolation weights (by G. Toth, 95)
        # Usage: triangulate xirreg yirreg [nx ny] 
        # Setup a regular grid enclosing the irregular and calculate
        # the neighbors nei1,nei2,nei3 and weights wei1,wei2,wei3 for each
        # regular grid point. The triangulation arrays are used in trigrid.
        # Vectors nx,ny,x,y,Aix,Aiy,Ax,Ay are also set for the regular grid.
		foreach 0 {xirr yirr indirr indreg ind0 iind0 jnd0 ind iind\
		   jxiy hxiy ixjy ixhy dst1 dst2 dst3 dd ww nn} {
		      SET $0 LOCAL
		}
		set xirr=$1 set yirr=$2
		if(!$?3){define 3 100}
		if(!$?4){define 4 100}
		set nx=$3 set ny=$4
		if(nx < 2 || ny < 2) {
		   echo "triangulate: the output image must be at least 2x2 (not $!(nx)x$!(ny))"
		   return
		}
		regulargrid xirr yirr nx ny x y
		set indreg=0,nx*ny-1
		set Aiy=int(indreg/nx+0.0001) set Aix=indreg-nx*Aiy
		set Ax=x[Aix] set Ay=y[Aiy]
		# Initialize neigbor and distance arrays for each regular grid point
		# Distance is measured in norm 1, ie. abs(dx)+abs(dy)
		set distmax local
		set distmax=2*(x[nx-1]-x[0]+y[ny-1]-y[0])+1
		set nei1=0*Aix set dst1=distmax+0*Aix
		set nei2=nei1  set dst2=dst1
		set nei3=nei1  set dst3=dst1
		# Assign an index to each irregular point. 
		set nn=1,dimen(xirr)
		# Indices of nn for indirect indexing
		set iind0=0,dimen(nn)-1 
		# Indices of regular points next to the irregular points
		set indirr=int((xirr-x[0])/(x[1]-x[0])-0.0001)+\
		    nx*int((yirr-y[0])/(y[1]-y[0])-0.0001)
		    # To avoid index out of range an extra irregular point is inserted
		    set xirr=0 concat xirr
		    set yirr=0 concat yirr
		    define verbose local define verbose | VERBOSE 0
		    # Put neighbor info into regular points around irregular points
		    triinic 0 0
		    triinic 0 1
		    triinic 1 0
		    triinic 1 1
		    # Setup shifted index arrays for propagation
		    set jxiy=Aix+(Aix<nx-1)+nx*Aiy
		    set hxiy=Aix-(Aix>0)+nx*Aiy
		    set ixjy=Aix+nx*(Aiy+(Aiy<ny-1))
		    set ixhy=Aix+nx*(Aiy-(Aiy>0))
		    # Propagate neighbor info into all regular points
		    set count=100
		    triiter
		    if(count<0){echo Did not iterate after 100 iteration}
		    VERBOSE $verbose
		    # Transform distances to weights: wi=1/di
		    set wei1=dst2*dst3 set wei2=dst1*dst3 set wei3=dst1*dst2
		    #  Normalize weights
		    set ww=wei1+wei2+wei3
		    set wei1=wei1/ww set wei2=wei2/ww set wei3=wei3/ww
		    # Change neighbor number into index number
		    set nei1=nei1-1 set nei2=nei2-1 set nei3=nei3-1
triinic 2 ## Determine neighbor and distance info for closest regular points
        # Assign nei1 and dst1 of regular grid points at indirr+$1+nx*$2
		set ind0=indirr+$1+nx*$2
		set dd=abs(xirr[nn]-Ax[ind0])+abs(yirr[nn]-Ay[ind0])
		# Propagate nn and dd into ind0. We allow for four irregular
		# points in a square, each gets a chance to propagate into a corner
		trinei1
		trinei1
		trinei1
		trinei1
triiter		## Determine neighborhood and distance info in the whole grid
		set count=count-1
		if(count<0){return}
		set ind0=indreg if(nei3==0 && nei2>0)
		trineis 3
		set ind0=indreg if(nei2==0 && nei1>0)
		trineis 2
		set ind0=indreg if(nei1==0)
		trineis 1
		set ind0=indreg if(nei3==0)
		if(dimen(ind0)>0){triiter}
		# Return when all points have third neighbors nei3 calculated
trineis	1	## Propagates $1th neighbor of ind0 from all directions
		if(dimen(ind0)>0){
		   set iind0=0,dimen(ind0)-1
		   trinei $1 jxiy
		   trinei $1 hxiy
		   trinei $1 ixjy
		   trinei $1 ixhy
		}
		# extra line for SM
trinei	2	## Propagates $1th neighbor of ind0 from shifted $2[ind0] locations 
		set jnd0=$2[ind0]
		do i=$1,1,-1{
		   set nn=nei$i[jnd0]
		   set dd=(nn>0)?abs(Ax[ind0]-xirr[nn])+abs(Ay[ind0]-yirr[nn]):distmax
		   trinei$1
		}
		# extra line for SM
trinei1		## Propagates 1st neighbor of ind0 from nn and dd
		set iind=iind0 if(dst1[ind0]>dd)
		triset 1
trinei2		## Propagates 2nd neighbor of ind0 from nn and dd
		set iind=iind0 if(dst2[ind0]>dd && nei1[ind0]!=nn)
		triset 2
trinei3		## Propagates 3rd neighbor of ind0 from nn and dd
		set iind=iind0 if(dst3[ind0]>dd && nei1[ind0]!=nn && nei2[ind0]!=nn)
		triset 3
triset	1	## Set $1th neighbor ans distance of ind0[iind] from nn and dd
		if(dimen(iind)>0){
		   set ind=ind0[iind]
		   set nei$1[ind]=nn[iind]
		   set dst$1[ind]=dd[iind]
		}
		# extra line for SM
