;+
;  NAME: spherical_image_dilate
; 
;  PURPOSE:
;    This function dilates binary data mapped on a spherical surface using a
;    circular structural element (kernel).
; 
;  CALLING SEQUENCE:
;    result=spherical_image_dilate(sph_image,radius) 
; 
;  INPUTS:
;    sph_image = a structure of type spherical_image_data (see
;                spherical_image_data__define.pro) with all fields defined.
;                If sph_image.image is not binary, then generate a binary
;                by doing a abs(*sph_image.image) ge 1.
;    radius = the radius of the circular dilation kernel in degrees
; 
;  OUTPUTS:
;    result = a structure of type spherical_image_data (see
;             spherical_image_data__define.pro) containing the dilated image.
;             This image is a binary image.
;
;  NOTES:
;    1.  Future improvement: allow arrays of structures and/or radii to be
;        passed into this routine.
;    2.  Another future improvement: allow routine to work on grayscale images
;        (i.e. assign to each point the maximum value of the surrounding
;        pixels defined by the kernel)
; 
;  MODIFICATION HISTORY:
;    M.DeRosa - 15 Dec 2005 - created
; 
;-

function spherical_image_dilate,sph_image,radius

;  usage message
if n_elements(sph_image) eq 0 then begin
  print,'  ERROR in spherical_image_dilate: no input data provided' 
  print,'  Calling sequence: result=spherical_image_dilate(sph_image,radius)'
  return,{spherical_image_data}
endif

;  get radius in radians
if n_elements(radius) eq 0 then begin
  print,'  ERROR in spherical_image_dilate: no dilation radius provided'
  return,{spherical_image_data}
endif
rad=radius(0)*!dpi/180
if rad le 0 then begin
  print,'  ERROR in spherical_image_dilate: negative dilation radius' 
  return,{spherical_image_data}
endif

;  get a binary input image
imin=abs(*sph_image.image) ge 1

;  create grids holding the theta and phi values of each point
nlat=sph_image.nlat
nlon=sph_image.nlon
phgr=(*sph_image.phi)#replicate(1,nlat)
thgr=replicate(1,nlon)#(*sph_image.theta)

;  loop through each row (latitude) in input image
result=bytarr(nlon,nlat)
for i=0,sph_image.nlat-1 do begin

  ;  compute spherical distances and determine latitude-dependent kernel
  anga=thgr(0,i)  ;  angle subtended from pole to current latitude
  angb=thgr  ;  angle subtended from poie to each point in the map
  capc=phgr(0,i)-phgr  ;  angle on surface between meridian lines
  sdist=acos( 0> (cos(anga)*cos(angb)+sin(anga)*sin(angb)*cos(capc)) <1 )
  kernel=sdist le rad

  ;  do convolution for current latitude
  wh=where(imin(*,i),nwh)
  if nwh gt 0 then for j=0,nwh-1 do result=result or shift(kernel,[wh(j),0])

endfor

;  format result
sph_dilated=sph_image
sph_dilated.image=ptr_new(result)
return,sph_dilated

end