#!/bin/sh

############################################################################
#
# MODULE:	i.landsat.rgb
#
# AUTHOR(S):	Markus Neteler. <neteler itc it>
#		Hamish Bowman, scripting enhancements
#
# PURPOSE:      create pretty LANDSAT RGBs: the trick is to remove outliers 
#               using percentiles (area under the histogram curve)
#
# COPYRIGHT:	(C) 2006-2012 by the GRASS Development Team
#
#		This program is free software under the GNU General Public
#		License (>=v2). Read the file COPYING that comes with GRASS
#		for details.
#
# TODO: implement better brightness control
#############################################################################

#%Module
#%  description: Performs auto-balancing of colors for LANDSAT images.
#%  keywords: raster, imagery, colors
#%End
#%option
#% key: red
#% type: string
#% gisprompt: old,cell,raster
#% description: LANDSAT red channel
#% required : yes
#%end
#%option
#% key: green
#% type: string
#% gisprompt: old,cell,raster
#% description: LANDSAT green channel
#% required : yes
#%end
#%option
#% key: blue
#% type: string
#% gisprompt: old,cell,raster
#% description: LANDSAT blue channel
#% required : yes
#%end
#%option
#% key: strength
#% type: double
#% description: Cropping intensity (upper brightness level)
#% options: 0-100
#% answer : 98
#% required : no
#%end
#%flag
#% key: f
#% description: Extend colors to full range of data on each channel
#%end
#%flag
#% key: p
#% description: Preserve relative colors, adjust brightness only
#%end
#%flag
#% key: r
#% description: Reset to standard color range
#%end
#%flag
#% key: s
#% description: Process bands serially (default: run in parallel)
#%end


if [ -z "$GISBASE" ] ; then
 echo "You must be in GRASS GIS to run this program." >&2
 exit 1
fi

# save command line
if [ "$1" != "@ARGS_PARSED@" ] ; then
    CMDLINE=`basename "$0"`
    for arg in "$@" ; do
        CMDLINE="$CMDLINE \"$arg\""
    done
    export CMDLINE
    exec g.parser "$0" "$@"
fi

# check if we have awk
if [ ! -x "`which awk`" ] ; then
    g.message -e "awk required, please install awk or gawk first"
    exit 1
fi 

# setting environment, so that awk works properly in all languages
unset LC_ALL
LC_NUMERIC=C
export LC_NUMERIC

PROG=`basename "$0"`

BLUE="$GIS_OPT_BLUE"
GREEN="$GIS_OPT_GREEN"
RED="$GIS_OPT_RED"

# 90 or 98? MAX value controls brightness
# think of percent (0-100), must be positive or 0
# must be more than "2" ?
BRIGHTNESS="$GIS_OPT_STRENGTH"


if [ 1 -eq $GIS_FLAG_F ] ; then
   for i in $RED $GREEN $BLUE ; do
      r.colors $i col=grey
   done
   exit 0
fi
if [ 1 -eq $GIS_FLAG_R ] ; then
   for i in $RED $GREEN $BLUE ; do
      r.colors $i col=rules << EOF
	0 black
	255 white
EOF
   done
   exit 0
fi

BRI_VAR_STR=percentile_`echo "$BRIGHTNESS" | awk '{printf("%.15g", $1)}' | tr '.' '_'`

if [ 0 -eq $GIS_FLAG_P ] ; then

   if [ 1 -eq $GIS_FLAG_S ] ; then
      # Unlocked, run serially
      for i in $RED $GREEN $BLUE ; do
   	 g.message "Processing <$i>..."
   	 eval `r.univar -ge "$i" perc=2,"$BRIGHTNESS" | grep '^percentile_'`
   	 MIN="$percentile_2"
   	 eval `echo "MAX=\"\\$\$BRI_VAR_STR\""`
   	 g.message -d message="<$i>:  min=$MIN   max=$MAX"
   	 r.colors $i col=rules << EOF
   	   0% black
   	   $MIN black
   	   $MAX white
   	   100% white
EOF
         done
   else
      # Unlocked, run in parallel
      g.message "Processing ..."
      eval `(
        r.univar -ge "$RED"   perc=2,"$BRIGHTNESS" | grep '^percentile_' | sed -e 's/^/R_/' &
        r.univar -ge "$GREEN" perc=2,"$BRIGHTNESS" | grep '^percentile_' | sed -e 's/^/G_/' &
        r.univar -ge "$BLUE"  perc=2,"$BRIGHTNESS" | grep '^percentile_' | sed -e 's/^/B_/'
        wait
      )`

      R_MIN="$R_percentile_2"
      eval `echo "R_MAX=\"\\$R_\$BRI_VAR_STR\""`
      G_MIN="$G_percentile_2"
      eval `echo "G_MAX=\"\\$G_\$BRI_VAR_STR\""`
      B_MIN="$B_percentile_2"
      eval `echo "B_MAX=\"\\$B_\$BRI_VAR_STR\""`

      r.colors "$RED" col=rules << EOF
	0% black
	$R_MIN black
	$R_MAX white
	100% white
EOF
      r.colors "$GREEN" col=rules << EOF
	0% black
	$G_MIN black
	$G_MAX white
	100% white
EOF
       r.colors "$BLUE" col=rules << EOF
	0% black
	$B_MIN black
	$B_MAX white
	100% white
EOF
   fi

else
   g.message "Processing ..."
   if [ 1 -eq $GIS_FLAG_S ] ; then
      # Bands are locked, run serially

      # we can combine all 3 bands into one r.univar call:
      # [disabled: uses 3x the RAM and actually runs slightly slower than 3 in a row]
      #eval `r.univar -ge "$RED,$GREEN,$BLUE" perc=2,"$BRIGHTNESS" | grep '^percentile_'`
      #ALL_MIN="$percentile_2"
      #eval `echo "ALL_MAX=\"\\$\$BRI_VAR_STR\""`

      # using 3 separate r.univar calls:
      ALL_MAX=0
      ALL_MIN=999999
      for i in $RED $GREEN $BLUE ; do
   	 g.message "Processing <$i>..."
   	 eval `r.univar -ge "$i" perc=2,"$BRIGHTNESS" | grep '^percentile_'`
   	 MIN="$percentile_2"
   	 eval `echo "MAX=\"\\$\$BRI_VAR_STR\""`
   	 g.message -d message="<$i>:  min=$MIN   max=$MAX"
   	 ALL_MAX=`echo "$MAX $ALL_MAX" | awk '{if ($1 > $2) print $1; else print $2}'`
   	 ALL_MIN=`echo "$MIN $ALL_MIN" | awk '{if ($1 < $2) print $1; else print $2}'` 
      done
   else
      # Bands are locked, run in parallel
      eval `(
        r.univar -ge "$RED"   perc=2,"$BRIGHTNESS" | grep '^percentile_' | sed -e 's/^/R_/' &
        r.univar -ge "$GREEN" perc=2,"$BRIGHTNESS" | grep '^percentile_' | sed -e 's/^/G_/' &
        r.univar -ge "$BLUE"  perc=2,"$BRIGHTNESS" | grep '^percentile_' | sed -e 's/^/B_/'
        wait
      )`

      ALL_MAX=0
      ALL_MIN=999999
      for i in R G B ; do
         eval `echo "MIN=\"\\$\${i}_percentile_2\""`
         eval `echo "MAX=\"\\$\${i}_\$BRI_VAR_STR\""`
         g.message -d message="<$i>:  min=$MIN   max=$MAX"
         ALL_MAX=`echo "$MAX $ALL_MAX" | awk '{if ($1 > $2) print $1; else print $2}'`
         ALL_MIN=`echo "$MIN $ALL_MIN" | awk '{if ($1 < $2) print $1; else print $2}'` 
      done
   fi

   g.message -d message="all_min=$ALL_MIN   all_max=$ALL_MAX"
   for i in $RED $GREEN $BLUE ; do
      r.colors $i col=rules << EOF
	0% black
	$ALL_MIN black
	$ALL_MAX white
	100% white
EOF
   done
fi

# write cmd history (only if raster maps are located in the current mapset)
CURRENT_MAPSET=`g.gisenv MAPSET`
RED_MAPSET=`g.findfile -n elem=cell file="$RED" | grep mapset | cut -d'=' -f2`
BLUE_MAPSET=`g.findfile -n elem=cell file="$BLUE" | grep mapset | cut -d'=' -f2`
GREEN_MAPSET=`g.findfile -n elem=cell file="$GREEN" | grep mapset | cut -d'=' -f2`

if [ "$BLUE_MAPSET" = "$CURRENT_MAPSET" ] ; then
    r.support "$BLUE" history="${CMDLINE}"
fi

if [ "$GREEN_MAPSET" = "$CURRENT_MAPSET" ] ; then
    r.support "$GREEN" history="${CMDLINE}"
fi

if [ "$RED_MAPSET" = "$CURRENT_MAPSET" ] ; then
    r.support "$RED" history="${CMDLINE}"
fi

exit 0