## -*- 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+(Aix0)+nx*Aiy set ixjy=Aix+nx*(Aiy+(Aiy0)) # 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