;note - some non-obvious routines used are HISTR and POLATE
 ; h = HISTR(array)  computes the running sum of the histogram of array
 ; and then normalizes to 1.0 (i.e., divides by the # of points) and puts
 ; this result in the FP vector h
 ; yp = POLATE(x,y,xp) uses the relationship (x,y) to determine values
 ; of yp corresponding to xp using linear interpolation
 ; it is assumed that both x and xp are monotonic !
 ;===========================================================================
func histgiven,xin,hm
 ;returns an image which is the input image XIN remapped such that its
 ;histogram matches hm which must be a running sum histogram beginning at 0
 hmx=indgen(hm)
 xq=min(xin)
 xq=xin-xq
 hin=histr(xq)
 g=polate(hm,hmx,hin)
 return,g(xq)
 endfunc
  ;===========================================================================
func histmatch,xin,xm
 ;returns an image which is the input image XIN remapped such that its
 ;histogram matches that of the image in XM
 ;both inputs must be arrays but not necessarily the same size
 xmmin=min(xm)
 if symdtype(xin) ne 0 then hm=histr(xm - xmmin) else hm=histr(word(xm)-xmmin)
 hmx=indgen(hm) + xmmin
 xq=min(xin)
 xq=xin-xq
 if symdtype(xin) ne 0 then hin=histr(xq) else hin=histr(word(xq))
 hm(0) = hm(0) > hin(0)	;this avoids extrapolation problems at botton
 g=polate(hm,hmx,hin)
 yq=g(xq)
 xq=0		;save memory
 return,yq
 endfunc
 ;===========================================================================
func histmatchsub,xin,xinsub,xm
 ;similar to histmatch but does the histogram matching on subimages
 ;instead of full image
 xmmin=min(xm)
 if symdtype(xin) ne 0 then hm=histr(xm - xmmin) else hm=histr(word(xm)-xmmin)
 hmx=indgen(hm) + xmmin
 mq=min(xinsub)
 xq=xinsub-mq
 if symdtype(xin) ne 0 then hin=histr(xq) else hin=histr(word(xq))
 hm(0) = hm(0) > hin(0)	;this avoids extrapolation problems at botton
 g=polate(hm,hmx,hin)
 if symdtype(xin) ne 0 then {
 yq=g( ((xin-mq)>0)<(num_elem(g)-1)) + xmmin
 } else {
 yq=g( ((word(xin)-mq)>0)<(num_elem(g)-1)) + xmmin
 }
 xq=0		;save memory
 return,yq
 endfunc
 ;===========================================================================
func histeq,xin
 ;returns a byte image which is the input image XIN remapped such that it
 ;has an equally populated histogram (histogram equalization) over the range
 ;0 to !scalemax
 xq=min(xin)
 xq=xin-xq
 hin=histr(xq)
 hm=indgen(fltarr((!scalemax+1)))
 g=byte(word(polate(hm/!scalemax,hm,hin)
 return,g(xq)
 endfunc
 ;===========================================================================
func histgauss,xin,dw
 ;returns a byte image which is the input image XIN remapped such that it
 ;has a gaussian histogram with a gaussian width dw over the range
 ;0 to !scalemax
 ;if dw is not specified, a value of 3 is assumed
 if !narg lt 2 then dw=3.0
 xq=min(xin)
 xq=xin-xq
 hin=histr(xq)
 hx=indgen(fltarr((!scalemax+1)))
 hg=0.5*(erf( (hx-128.)*dw/128.)+1.0)
 ;we have to scale the min and max of hg to the 0,1 range
 hg=hg/(hg(!scalemax)-hg(0))	hg=hg-hg(0)	hg=hg<1.0
  $g=byte(word(polate(hg,hx,hin)
 return, $g(xq)
 endfunc
 ;===========================================================================
func pscale,x,bot,top,pow
 ;generalized version of scale, generates and uses lookup table
 ;10/6/96 modified to handle fp cases with narrow ranges
 ;for I*4 use the fp branch, otherwise we can get very slow performance
 ;if someone uses I*4's with very large values
 type = symdtype(x)
 if type lt 2 then { 
 ;for integer cases, just a finite set of cases
 n=abs(top-bot)+1
 ;indgen won't go above 2^24 for FP, so use lonarr here rather then fltarr
 g=indgen(lonarr(n))
 } else {
 ;for fp and dp, do 256 here, scale will reduce to screen range
 n = 256
 range = abs(top-bot)
 if range eq 0 then range = 1.0
 ;now invert range for a multiplier
 range = 256./range
 g=indgen(fltarr(n))
 }
 if pow ne 1.0 then g=g^pow
 g=scale(g)
 if top lt bot then g=reverse(g)
 offset=bot<top
 if type lt 2 then {
 if type eq 1 or type eq 0 then offset=word(rfix(offset))
 y=g( (x-offset)>0<(n-1) )
 } else {
 y = g( ((x-offset)*range) > 0 < 255)
 }
 return,y
 endfunc
 ;===========================================================================
func trim,xin,cutb,cutt,pq
 ;returns a byte image which is the input image XIN scaled over a range
 ;which cuts off the top and bottom fractions CUTB and CUTT
 ;if CUTT is not given, a value of 1.0-CUTB is used; if neither is given,
 ;then .001 (0.1%) is trimmed off both ends
 ncase !narg-1
 { cq=0.001	cq2=0.999	pow=1.0}
 { cq=cutb	cq2=1.0-cq	pow=1.0}
 { cq=cutb	cq2=cutt	pow=1.0}
 { cq=cutb	cq2=cutt	pow=pq }
 else { ty,'TRIM doesn''t accept',!narg,' arguments'	return,0 }
 endcase
 xq=min(xin)
 if symdtype(xin) ne 0 then {
   hin=histr(xin-xq)
   ;to get a good display when there is a "no data" value that is a
   ;large negative #, check the first entry
   if (xq le -10000) {
   	if (hin(0) lt 1.0) {
		hin = hin - hin(0)
		hin = hin/total(hin)	;renormalize
	}
   }	
 } else {
 ;the byte case
 xq=word(xq)
 hin=histr(xin)	hin=hin(xq:*)   }
 ;catch narrow range FP images here and just scale them
 if symdtype(xin) ge 3 and num_elem(hin) lt 20 then {
 range = [min(xin), max(xin)]
 } else {
 ;normal case
 hinx=indgen(hin)+xq
 if num_elem(hin) gt 1 then range=polate(hin,hinx,[cq,cq2]) else range=[hinx,hinx]
 ;note, if the bottom or top bins are populated more than the trim factor,
 ;we have a problem, therefore insure that the range doesn't exceed
 ;the data limits as given by hinx(0) and max(hinx)
 range(0)=range(0)>hinx(0)
 range(1)=range(1)<max(hinx)
 }
 ;type,'range',range
 $trim_range = range		;save for possible use
 if pow ne 1.0 then return,pscale(xin,range(0),range(1),pow)
 return,scale(xin,range(0),range(1))
 endfunc
 ;===========================================================================
func trimtell,xin,cutb,cutt,pq
 ;returns a byte image which is the input image XIN scaled over a range
 ;which cuts off the top and bottom fractions CUTB and CUTT
 ;if CUTT is not given, a value of 1.0-CUTB is used; if neither is given,
 ;then .001 (0.1%) is trimmed off both ends
 ncase !narg-1
 { cq=0.001	cq2=0.999	pow=1.0}
 { cq=cutb	cq2=1.0-cq	pow=1.0}
 { cq=cutb	cq2=cutt	pow=1.0}
 { cq=cutb	cq2=cutt	pow=pq }
 else { ty,'TRIM doesn''t accept',!narg,' arguments'	return,0 }
 endcase
 xq=min(xin)
 if symdtype(xin) ne 0 then {
   hin=histr(xin-xq)
   ;to get a good display when there is a "no data" value that is a
   ;large negative #, check the first entry
   if (xq le -10000) {
   	if (hin(0) lt 1.0) {
		hin = hin - hin(0)
		hin = hin/total(hin)	;renormalize
	}
   }	
 } else {
 ;the byte case
 xq=word(xq)
 hin=histr(xin)	hin=hin(xq:*)   }
 hinx=indgen(hin)+xq
 ;catch narrow range FP images here and just scale them
 if symdtype(xin) ge 3 and num_elem(hin) lt 20 then {
 range = [min(xin), max(xin)]
 } else {
 ;normal case
 if num_elem(hin) gt 1 then range=polate(hin,hinx,[cq,cq2]) else range=[hinx,hinx]
 ;note, if the bottom or top bins are populated more than the trim factor,
 ;we have a problem, therefore insure that the range doesn't exceed
 ;the data limits as given by hinx(0) and max(hinx)
 range(0)=range(0)>hinx(0)
 range(1)=range(1)<max(hinx)
 }
 type,'range',range
 $trim_range = range		;save for possible use
 if pow ne 1.0 then return,pscale(xin,range(0),range(1),pow)
 return,scale(xin,range(0),range(1))
 endfunc
 ;===========================================================================
func ft,xin,cutb,cutt,pq
 ;1/17/91 name change from fasttrim to ft (tired of typing the longer one)
 ;returns a byte image which is the input image XIN scaled over a range
 ;which cuts off the top and bottom fractions CUTB and CUTT
 ;if CUTT is not given, a value of 1.0-CUTB is used; if neither is given,
 ;then .001 (0.1%) is trimmed off both ends
 ncase !narg-1
 { cq=0.001	cq2=0.999	pow=1.0}
 { cq=cutb	cq2=1.0-cq	pow=1.0}
 { cq=cutb	cq2=cutt	pow=1.0}
 { cq=cutb	cq2=cutt	pow=pq }
 else { ty,'FT doesn''t accept',!narg,' arguments'	return,0 }
 endcase
 ;use a subset of the image for the histogram
 nx=dimen(xin,0)		ny=dimen(xin,1)
 if nx gt 10 then {
 np=(ny/10)>1	dy=ny/np	im=indgen(lonarr(1,np))*nx*dy+10*nx
 in=lonarr(nx,np)	in=indgen(in,0)
 in=in+im	subim=xin(in) } else	subim=xin
 xq=min(subim)  mxq = !lastmax
 sfac = 1.0
 if symdtype(subim) gt 1 then {
   ;consider range reduction if we have a very large spread here
   if (mxq - xq) gt 200000 then {
     sfac  = 200000.0/(mxq - xq)
     subim = subim * sfac
     xq = xq * sfac
     ;this will also make subim FP
     sfac = 1.0/sfac	;for use on the ranges later
     ty,'sfac =', sfac
   }
 }
 if symdtype(subim) ne 0 then {
   ;3/25/98 - if the range is too large, we cannot subtract a large negative min
   ;from a large positive value without rollover, but safe for any positive min
   if xq lt 0 then hin=histr(subim) else hin=histr(subim-xq)
   ;when submin has neg values, the histr routine accounts for them but if
   ;all positive, the first hin will be for 0
   ;to get a good display when there is a "no data" value that is a
   ;large negative #, check the first entry
   if (xq le -10000) {
   	if (hin(0) lt 1.0) {
		hin = hin - hin(0)
		hin = hin/max(hin)	;renormalize
	}
   }	
 } else {
   ;the byte case
   xq=word(xq)
   hin=histr(subim)	hin=hin(xq:*)
 }
 if isscalar(hin) then hin = [hin]
 hinx=indgen(hin)+xq
 ;catch narrow range FP images here and just scale them
 if symdtype(xin) ge 3 and num_elem(hin) lt 20 then {
   range = sfac * [min(xin), max(xin)]
 } else {
   ;normal case
   if num_elem(hin) gt 1 then range = sfac * polate(hin,hinx,[cq,cq2]) else range = sfac * [hinx,hinx]
   ;note, if the bottom or top bins are populated more than the trim factor,
   ;we have a problem, therefore insure that the range doesn't exceed
   ;the data limits as given by hinx(0) and max(hinx)
   range(0)=range(0)> (sfac * hinx(0))
   range(1)=range(1)< (sfac * max(hinx))
 }
 ;;type,'range',range
 $trim_range = range		;save for possible use
 r0 = range(0)		r1 = range(1)
 ;;ty,'r0, r1 =', r0, r1
 if pow ne 1.0 then return,pscale(xin, r0, r1, pow)
 ;conditional below to bypass bug with F*8 arrays, remove when distributions debugged
 ;;if symdtype(xin)  eq 4 then return, scale(float(xin), r0, r1)
 ;;  else return, scale(xin, r0, r1)
 return, scale(xin, r0, r1)
 endfunc
 ;===========================================================================
func ftsigned(xin, cutt, pq)
 ;a variation of ft for signed quantities (like QUV) where we want a +/- limit applied
 ;only cutt is passed
 cq = 0
 ncase !narg-1
   { cq2=0.999	pow=1.0}
   { cq2=cutt	pow=1.0}
   { cq2=cutt	pow=pq }
   else { ty,'FT doesn''t accept',!narg,' arguments'	return,0 }
 endcase
 ;use a subset of the image for the histogram
 nx=dimen(xin,0)		ny=dimen(xin,1)
 if nx gt 10 then {
   np=(ny/10)>1      dy=ny/np          im=indgen(lonarr(1,np))*nx*dy+10*nx
   in=lonarr(nx,np)  in=indgen(in,0)
   in=in+im          subim=xin(in)
 } else	subim=xin
 ;we have to do this eventually, do now before other tests
 subim = abs(subim)
 ;here we use 0 for the min
 xq=0
 mxq = max(subim)
 sfac = 1.0
 if symdtype(subim) gt 1 then {
   ;consider range reduction if we have a very large spread here
   if (mxq - xq) gt 200000 then {
     sfac  = 200000.0/(mxq - xq)
     subim = subim * sfac
     xq = xq * sfac
     ;this will also make subim FP
     sfac = 1.0/sfac	;for use on the ranges later
     ty,'sfac =', sfac
   }
 }
 if symdtype(subim) ne 0 then {
   ;3/25/98 - if the range is too large, we cannot subtract a large negative min
   ;from a large positive value without rollover, but safe for any positive min
   if xq lt 0 then hin=histr(subim) else hin=histr(subim-xq)
   ;when submin has neg values, the histr routine accounts for them but if
   ;all positive, the first hin will be for 0
   ;to get a good display when there is a "no data" value that is a
   ;large negative #, check the first entry
   if (xq le -10000) {
   	if (hin(0) lt 1.0) {
		hin = hin - hin(0)
		hin = hin/max(hin)	;renormalize
	}
   }	
 } else {
   ;the byte case
   xq=word(xq)
   hin=histr(subim)	hin=hin(xq:*)
 }
 if isscalar(hin) then hin = [hin]
 hinx = indgen(hin)+xq
 ;catch narrow range FP images here and just scale them
 if symdtype(xin) ge 3 and num_elem(hin) lt 20 then {
   range = sfac * [min(xin), max(xin)]
 } else {
   ;normal case
   if num_elem(hin) gt 1 then range = sfac * polate(hin,hinx,[cq,cq2]) else range = sfac * [hinx,hinx]
   ;note, if the bottom or top bins are populated more than the trim factor,
   ;we have a problem, therefore insure that the range doesn't exceed
   ;the data limits as given by hinx(0) and max(hinx)
   range(1) = range(1)< (sfac * max(hinx))
   ;and for this symmetric mode for signed data, we have
   range(0) = -range(1)
 }
 ;;type,'range',range
 $trim_range = range		;save for possible use
 r0 = range(0)		r1 = range(1)
 ;;ty,'r0, r1 =', r0, r1
 if pow ne 1.0 then return,pscale(xin, r0, r1, pow)
 ;conditional below to bypass bug with F*8 arrays, remove when distributions debugged
 ;;if symdtype(xin)  eq 4 then return, scale(float(xin), r0, r1)
 ;;  else return, scale(xin, r0, r1)
 return, scale(xin, r0, r1)
 endfunc
 ;===========================================================================
subr setup_convert,xin,mode,cutb,cutt
 ;
 ;intended for integer images, we set up a lookup table to allow scaling
 ;images given a "sample" image in xin
 ;xin is used only to establish the bottom and top cutoff points
 ;the lookup table array is put in  $CONVERT, the offset in  $OFFSET,
 ;and a scaled image is obtained using the sister function CONVERT
 ;
 ;mode controls 3 types of mapping, 0 is linear, 1 is sqrt, and 2 is log
 ;note that using CONVERT with these lookup tables is much faster than
 ;arithmetic conversions
 ;
 ;if CUTT is not given, a value of 1.0-CUTB is used; if neither is given,
 ;then .01 (1%) is trimmed off both ends
 ;
 ncase !narg-1
 { cq=0.001	cq2=0.999	mode=0 }
 { cq=0.001	cq2=0.999 }
 { cq=cutb	cq2=1.0-cq }
 { cq=cutb	cq2=cutt }
 else { ty,'SETUP_CONVERT doesn''t accept',!narg,' arguments'	return,0 }
 endcase
 xq=min(xin)  if symdtype(xq) eq 0 then xq=word(xq)
 if symdtype(xin) ne 0 then hin=histr(xin-xq) else hin=histr(word(xin)-xq)
 hinx=indgen(hin)+xq
 range=polate(hin,hinx,[cq,cq2])
  $offset=range(0)
  $convert=indgen(fltarr(range(1)-range(0)+1))
 ncase mode-1
 {  $convert=sqrt( $convert) }
 {  $convert=alog( $convert+1.0) }
 endcase
  $convert=scale( $convert)
 endsubr
 ;===========================================================================
func convert,x
 ;setup_convert must have been run first
 ;uses the lookup table and offset in  $CONVERT and  $OFFSET
 xq= $convert(x- $offset)
 return,xq
 endfunc
 ;===========================================================================
subr setup_pconvert,xin,pow,cutb,cutt
 ;
 ;intended for integer images, we set up a lookup table to allow scaling
 ;images given a "sample" image in xin
 ;xin is used only to establish the bottom and top cutoff points
 ;the lookup table array is put in  $CONVERT, the offset in  $OFFSET,
 ;and a scaled image is obtained using the sister function CONVERT
 ;
 ;pow is the power (or gamma) of the image, linear is 1
 ;note that using CONVERT with these lookup tables is much faster than
 ;arithmetic conversions
 ;
 ;if CUTT is not given, a value of 1.0-CUTB is used; if neither is given,
 ;then .01 (1%) is trimmed off both ends
 ;
 ncase !narg-1
 { cq=0.01	cq2=0.99	pow=1.0 }
 { cq=0.01	cq2=0.99 }
 { cq=cutb	cq2=1.0-cq }
 { cq=cutb	cq2=cutt }
 else { ty,'SETUP_PCONVERT doesn''t accept',!narg,' arguments'	return,0 }
 endcase
 xq=min(xin)
 hin=histr(xin-xq)	hinx=indgen(hin)+xq
 range=polate(hin,hinx,[cq,cq2])
  $offset=range(0)
  $convert=indgen(fltarr(range(1)-range(0)+1))
 if pow ne 1.0 then  $convert= $convert^pow
  $convert=scale( $convert)
 endsubr
 ;===========================================================================
func pconvert,x
 ;setup_pconvert must have been run first
 ;uses the lookup table and offset in  $CONVERT and  $OFFSET
 ;note that the older CONVERT uses the same variables so don't use both at the
 ;same time !
 xq= $convert(x- $offset)
 return,xq
 endfunc
 ;===========================================================================