#! /bin/bash -e
# Last edited on 2026-01-18 22:55:27 by stolfi

# Computes the histogram of the grayscale PNG file ${ifile}, then plots
# it and its cumulative graph as a PNG file that is written to standard
# output.
#
# If ${hcrop} is not "NONE", it must be an ImageMagick "-crop"
# argument. In this case, considers only the rectangle of the image
# specified by that argment.
#

ifile="$1"; shift
hcrop="$1"; shift
maxval="$1"; shift
pct="$1"; shift

temp="/tmp/$$"

# Convert the file to PGM to compute the histogram:
echo "  computing histogram ..." 1>&2
if [[ "/${hcrop}" == "/NONE" ]]; then
  crop_opts=( )
else
  crop_opts=( "-crop" "${hcrop}" )
  echo "histogram crop_opts = ( ${crop_opts[@]} )" 1>&2
fi
tfile="${temp}-page.pgm" 
hfile="${temp}.hist"
convert ${ifile} \
    -set colorspace LinearGray \
    ${crop_opts[@]} \
    -colorspace LinearGray \
    ${tfile}
# display ${tfile}

pnmxhist ${tfile} > ${hfile}

export GDFONTPATH="${HOME}/ttf"

echo "  creating plots ..." 1>&2
for np in 0 1 ; do
  tfile_np="${temp}-plot-${np}.png"
  if [[ ${np} -eq 0 ]]; then
    sep=''; level=''
  else
    sep=','; level='(0.99) with lines lc rgb "#ff0000"'
  fi
  gnuplot <<EOF
    set term png size 4000,800 linewidth 2.0 font "Arial" 20
    set output "${tfile_np}"
    np = ${np}
    maxval = ${maxval}
    xoff = int(maxval/50 + 1)
    set xrange [(-xoff):(maxval+xoff)]
    
    # Choose major tics:
    xt0 = 100  # Major X tic spacing.
    xt1 = 5
    # while (10*xt0 < maxval) { xt0 = 10*xt0; }
    # 
    # if (maxval < 2*xt0) { 
    #   xt1 = 20
    # } else if (maxval < 5*xt0) { 
    #   xt1 = 10
    # } else if (maxval < 10*xt0) {
    #   xt1 = 5
    # } else {
    #   xt1 = 2
    # }
    set xtics xt0
    set mxtics xt1
    set grid xtics mxtics lc rgb "#44aaff", lc rgb "#aaddff"
    set grid ytics lc rgb "#44aaff", lc rgb "#aaddff"
    set lmargin 7
    xmin = 1.0e-14
    xmax = 1.0 - xmin
    xclip(x) = (x < xmin ? xmin : (x > xmax ? xmax : x)) 
    blog(x)= log(xclip(x)) - log(xclip(1-x))
    golb(y) = xclip(exp(y)/(1 + exp(y)))
    ymin = blog(1.0e-7)
    ymax = blog(1 - 1.0e-7)
    # printerr xmin, xmax, ymin, ymax, blog(0.001), blog(0.999)
    pct = ${pct}
    tcp = 1 - pct
    # printerr "pct =", pct, "tcp =", tcp
    if (np == 0) {
      # Absolute histogram
      set yrange [0.1:]; set logscale y;
      plot "${hfile}" using 1:2 title "count" with histeps lc rgb "#dd2200"
    } else {
      # Relative cumulative histograms:
      set yrange [(1.05*ymin):(1.05*ymax)]
      set ytics ( \
        "1e-9" (blog(1.0e-9)) 0, \
        "1e-8" (blog(1.0e-8)) 0, \
        "1e-7" (blog(1.0e-7)) 0, \
        "1e-6" (blog(1.0e-6)) 0, \
        "1e-5" (blog(1.0e-5)) 0, \
        "1e-4" (blog(1.0e-4)) 0, \
        "1e-3" (blog(1.0e-3)) 0, \
        "1e-2" (blog(1.0e-2)) 0, \
        "1e-1" (blog(1.0e-1)) 0, \
        "0.5" 0.0, 0, \
        "1-1e-1" (blog(1-1.0e-1)) 0, \
        "1-1e-2" (blog(1-1.0e-2)) 0, \
        "1-1e-3" (blog(1-1.0e-3)) 0, \
        "1-1e-4" (blog(1-1.0e-4)) 0, \
        "1-1e-5" (blog(1-1.0e-5)) 0, \
        "1-1e-6" (blog(1-1.0e-6)) 0, \
        "1-1e-7" (blog(1-1.0e-7)) 0, \
        "1-1e-8" (blog(1-1.0e-8)) 0, \
        "1-1e-9" (blog(1-1.0e-9)) 0  \
      )
      set yzeroaxis lt 1 lw 1.5 lc rgb "#2288ff"
      pct_title = sprintf(" %6.4f fractile", pct)
      tcp_title = sprintf(" %6.4f fractile", tcp)
      plot \
        (blog(pct))                   title (pct_title) with lines lc rgb "#aa9966", \
        (blog(tcp))                   title (tcp_title) with lines lc rgb "#aa9966", \
        "${hfile}" using 1:(blog(column(5))) title "accum fraction" with lines lc rgb "#dd2200"
    }
    quit
EOF
done

echo "  joining plots ..." 1>&2
convert \
  -append \
  ${temp}-plot-0.png \
  ${temp}-plot-1.png \
  -resize '50%' \
  PNG:-

rm -f ${temp}*.*