## -*-SM-*- Set SM-mode in emacs
entropy	13	# Find the entropy of a vector. Usage: entropy vec [min max]
		# If min and max are provided, take them as the minimum
		# and maximum values of the vector (if omitted, they
		# will be found for you). E.g. calc entropy(image[*,*])
		define min local  define max local  
		if($?2) {
		   if(!$?3) {
		      user abort "$0: if you specify min you must specify max"
		   }
		   define min $2  define max $3
		} else {
		   vecminmax $1 min max
		}
		local set h=histogram($1: do($min, $max))/dimen($1)
		set $0=-sum(h*(h == 0 ? 0 : ln(h)/ln(2)))
		#
hist_equalise 2	# Histogram equalise an image. 
		# Usage: hist_equalise min max
		set ii local  set h local
		set ii=image[*,*]
		set h=histogram(ii:do($1-1,$2+1))
		set h[0] = 0  set h[dimen(h)-1] = 0
		set ii=(ii < $1) ? $1 : (ii > $2 ? $2 : ii)
		set h=cumulate(h)/sum(h)
		set image[*,*] = h[ii-$1]
		#
gammaCorrect	3	# Gamma correct an image, for values between $1 and $2;
		# $3 is the value of gamma
		set ii local
		set ii = image[*,*]
		set ii = (ii < $1) ? 0 : (ii > $2 ? 1 : (ii - $1)/($2 - $1))
		set ii = $1 + ($2 - $1)*(1 + $3)*ii/(1 + $3*ii)
		set image[*,*] = ii
		#
contour_dot  contour_points
contour_points 37  #A macro to make contour and dot plots from color-color diagrams.
		# Usage: contour_dot flags x y [binx biny mincontour logical]
		# If flags & 0x1    Don't plot points
		# If flags & 0x2    Don't plot contours
	        #Arguments are quantities in x and y, and optionally,
		#binsize in x and y (0 => default), base contour level
		# in units of the contour interval, and logical
		#condition for selecting the sources.  Having the
		#points drawn ptype 0 0 is the default; if you want a
		#different ptype, set the variable set_ptype to 1 and
		#your chosen value will be respected
		local define flags ($1)

		local define noPoints   ($flags & 0x1)
		local define noContours ($flags & 0x2)

		local set _x = $2
		local set _y = $3

		local define binsizex 0.02
		local define binsizey 0.02
		if ($?4 AND $4 != 0) {
		   define binsizex $4
		   if ($?5 AND $5 != 0) {
		      define binsizey $5
		   } else {
		      define binsizey $binsizex
		   }
		} 
		local define mincontour 0.2
		if ($?6) {
		   define mincontour $6
		}

		if ($?7) {
		   local set selection = $7
		   set _x = _x if (selection)
		   set _y = _y if (selection)
		}

		local define nbinx (int(abs((($fx2 - $fx1)/$binsizex))))
		local define nbiny (int(abs((($fy2 - $fy1)/$binsizey))))

		local set temp = 0*do(1,$nbinx*$nbiny)
		image ($nbinx,$nbiny) $fx1 $fx2 $fy1 $fy2
		set image[*,*] = temp 
		
		local set binx = abs(int((_x - $fx1)/$binsizex))
		local set biny = abs(int((_y - $fy1)/$binsizey))
		#This should be a one-dimensional quantity unique to
		#each binx, biny pair.
		local set onedindex = binx*$nbiny + biny
		local set ii = 0, $nbinx*$nbiny - 1
		local set number_in_bin = histogram(onedindex:ii)
		#Now we just have to unpack this
		local set j = 0, $nbiny - 1
		define i local
		do i = 0, $nbinx - 1 {
		   set image[$i,j] = number_in_bin[$i*$nbiny + j]*$nbinx*$nbiny/dimen(binx)     
		}

		#The histogram routine puts all objects outside the range in the
		#very last entry.  We can fix this as follows:	
		set image[$nbinx-1,$nbiny-1] = image[$nbinx-2,$nbiny-2]

		#Can we smooth the contours a bit?
		local set x = 0, $nbinx - 1
		local set y = 0, $nbiny - 1
		define x local
		set cut local  set smoothcut local
		do x = 1, $nbinx - 1 {
		   set cut = image[$x,y]
		   smooth cut smoothcut 3
		   set image[$x,y] = smoothcut
		}

		define y local
		do y = 1, $nbiny - 1 {
		   set cut = image[x,$y]
		   smooth cut smoothcut 3
		   set image[x,$y] = smoothcut
		}
		
		#Next, we need to assign to each point a contour level. 
		set binx = (binx < 0)?(0):(binx)
		set binx = (binx > $nbinx -1)?($nbinx-1):(binx)
		set biny = (biny < 0)?(0):(biny)
		set biny = (biny > $nbiny -1)?($nbiny-1):(biny)
		local set local_density = image[binx,biny]

		define min local  define max local
		minmax min max
		local define interval (($max-$min)/9)

		if (!$noContours) {
		   #Draw contours
		   local set levs = $min + $interval*$mincontour, $max, $interval
		   levels levs
		   contour
		}
		
		if (!$noPoints) {
		   if (!$?set_ptype OR !$set_ptype) {
		      local define ptype |
		      ptype 0 0
		   }

		   #Dither a bit to get rid of discretization
		   local set temp1 = _x + 0.01*(random(dimen(_x)) - 0.5)
		   local set temp2 = _y + 0.01*(random(dimen(_y)) - 0.5)
		   local set draw_points = (local_density < $min + $interval*$mincontour)
		   points temp1 temp2 if(draw_points)
		   if ($verbose) {
		      echo The total number of points plotted is $(sum(draw_points))
		   }
		   
		   if (!$?set_ptype OR !$set_ptype) {
		      local set pt = <$ptype>
		      ptype $(pt[0]) $(pt[1])
		   }
		}

