#! /usr/bin/python3
# Last edited on 2025-12-12 06:44:05 by stolfi

  
def split_component(arg):
  # The {arg} must be a string "{band}/{maxval}".
  # Returns a tuple {({band}, {maxval})} with {maxval} as {int}.
  flds = re.split(r"/", arg)
  if len(flds) != 2: arg_error(f"bad component {arg}")
  band = flds[0]
  maxval = int(flds[1])
  return (band, maxval,)
  # ----------------------------------------------------------------------


######################################################################

  # Add the minimum luminosity and clip the chroma:
  sys.stderr.write(f"adding the minimum luminosity ...\n")
  B = add_luminosity(D, luma_fg)
  
  # Shift and rotate {\RR^{nd}} so that {bg_avg} goes to the origin
  # {(0,0,...,0)} and {bg_pc0} becomes the luma axis {(1,1,...,1)}:
  sys.stderr.write(f"rotating the luma axis ...\n")
  B = rotate_main_axis(B, bg_avg, bg_pc0) # ???
  
  # Normalize chroma coordinates to unit dev along each princ axis:
  sys.stderr.write(f"normalizing chroma coordinates ...\n")
  mchr = normalize_pure_chroma(D, M) # ???

def rotate_main_axis(B, bg_avg, bg_pc0):
  
  # Create the matrix {V} that projects {\RR^{nd}} on the subspace
  # orthogonal to the "luminance" axis {(1,1,1,...,1)}:
  
  V = numpy.full((nb,1,), 1.0/sqrt(nb))
  Q = numpy.identity(nb) - (V @ numpy.transpose(V))
  assert numpy.shape(Q) == (nb,nb,), "bug Q shape"
  PQ = P @ Q
  
  assert abs(numpy.sum(dpix, 0)) < 1.0e-6, "bug Q projection"
  assert False
  return
  # ----------------------------------------------------------------------

######################################################################

  # Determine the "maximum" relative saturation of any pixel in {M}:
  sys.stderr.write("extracting chromas of relevant pixels ...\n")
  sats = []
  for iy in range(ny):
    for ix in range(nx):
      if M[iy,ix] != 0:
        fpix = F[iy,ix,:]
        # Get luma {flum} and pure chroma vector {fchr} of {fpix}:
        flum = numpy.mean(fpix)
        assert flum > 0.1 and flum < 0.9*luma_bg
        fdif = fpix - numpy.full((nb,), flum)
        fchr = numpy.linalg.norm(fdif)
        if fchr > 1.0e-6: sats.append(fchr/flum)
  ns = len(sats)
  sys.stderr.write(f"collected {ns} rel saturations\n")
  assert ns >= 100, "bug too few non-gray pixels"

  frac_sat = 0.002
  msat = numpy.quantile(numpy.array(sats), 1-frac_sat, axis=0)
  sys.stderr.write(f"assumed max rel saturation before mixing = {msat:10.6}\n")
  

######################################################################
  
  sys.stderr.write(f"adjusting luminance ...\n")
  # Ranges of spectral averages after mapping:
  mlum_frac = 0.01 # Fraction of hi and lo outliers in selected pixels.
  mlum_min = 0.15  # Absolute min.
  mlum_lo = 0.25   # Low fractile of selected pixels.
  mlum_hi = 0.65   # High fractile of selected pixels.
  mlum_max = 0.75  # Absolute max.
  D = adjust_luminance(D, M, mlum_frac, mlum_lo, mlum_hi, mlum_min, mlum_max)


######################################################################

  # Collect the relevant spectra and average them:
  smin = numpy.zeros((nb,)); nmin = 0
  smax = numpy.zeros((nb,)); nmax = 0
  for ks in range(ns):
    spix = S[ks,:]
    tpix = numpy.sum(spix)
    if tpix >= talb_min_clip and tpix <= talb_min_take:
      smin += spix; nmin += 1
    if tpix >= talb_max_take and tpix <= talb_max_clip:
      smax += spix; nmax += 1
  debug("nmin", nmin)
  debug("nmax", nmax)
  assert nmin >= 1 and nmax >= 1, "no pixels in range"
  smin /= nmin;
  smax /= nmax;



def compute_deviations_of_analysis_data(masks, bands, savg, smin, smax, rdev, P):
  # For each {mask} in {masks[0..nm]}, writes a file with the analysis
  # results for the subset of samples of clip {clip} in page {page} selected by the 
  # mask {mask}, namely {savg,smin,smax,rdev,P} and spectral bands with names
  # {bands[0..nb-1]}. See the top of this module for the format of this
  # file.
  #
  nb, ne, nm = numpy.shape(P)
  assert len(masks) == nm, "inconsitent {masks}"
  assert len(bands) == nb, "inconsitent {bands}"
  assert numpy.shape(savg) == (nb,nm,)
  assert numpy.shape(smin) == (nb,nm,)
  assert numpy.shape(smax) == (nb,nm,)
  assert numpy.shape(rdev) == (nb,nm,)

  for im in range(nm):
    mask = masks[im]
    sum2_savg = numpy.zeros((nm,))
    sum2_smin = numpy.zeros((nm,))
    sum2_smax = numpy.zeros((nm,))
    sum2_P = numpy.zeros((ne,nm,))
    sum2_rdev = numpy.zeros((nm,))
    for ib in range(nb+1):
      band_tag = bands[ib] if ib < nb else "DEVIATIONS"
      efmt = " %+7.4f" if ib < nb else " %7.4f"

      for svar, sum2_svar in \
        ( (savg, sum2_savg),
          (smin, sum2_smin),
          (smax, sum2_smax),
          (rdev, sum2_rdev)
        ):
        svar_val = svar[ib,im] if ib < nb else sqrt(sum2_svar[im])
        wr.write(efmt % svar_val)
        if ib < nb: sum2_svar[im] += svar_val**2



      for ie in range(ne):
        P_val = P[ib,ie,im] if ib < nb else sqrt(sum2_P[ie,im])
        if ib < nb: sum2_P[ie,im] += P_val**2

  return
  # ----------------------------------------------------------------------


########################################################################

def read_normalization_mask_image(page, clip, mask):
  mfile = f"pages/{page}/clips/{clip}/masks/{mask}.png"
  tfile = ".temp.pgm"
  bash(f"convert {mfile} -depth 16 -colorspace LinearGray {tfile}")
  M = read_pgm_image(tfile, 65535)
  return M
  # ----------------------------------------------------------------------

########################################################################
  
  #  Full color range minus zero:
  hblack = "#000100010001"
  hwhite = "#ffffffffffff"
  
  # Multiply raw clip image by mask, eliminate other zeros:
  pgdir_in = f"MS/davis/{page}"
  cfile = f"{pgdir_in}/{band}/clips/{clip}.pgm"
  mfile = f"pages/{page}/clips/{clip}/masks/{mask}.png"
  xfile = f".{clip}.pgm"
  vfile = f".{clip}.vals"
  bash \
    ( \
      f"convert " + \
      f" \\( {cfile} -set colorspace LinearGray +level-colors '{hblack},{hwhite}' \\)" + \
      f" \\( {mfile} -set colorspace LinearGray -depth 16 -colorspace LinearGray \\)" + \
      f" -compose multiply" + \
      f" -composite" + \
      f" -depth 16 -colorspace LinearGray" + \
      f" {xfile}"
    )
    
  # Extract nonzero sample values in scanline order:
  bash \
    ( \
      f"cat {xfile}" + \
      r" | pnmtoplainpnm" + \
      r" | tr ' ' '\012'" + \
      r" | egrep -e '.'" + \
      r" | gawk '// { if ((FNR > 4) && ($1 > 0)) { print $1; } }'" + \
      f" > {vfile}"
    )
    
  # Read samples, convert to float in {[0 _ 1]}, make a list:
  samples = []
  with open(vfile) as rd:
    nbig = 0; ismp_max = 0
    for line in rd:
      line = re.sub(r"[#].*$", "", line)
      if line.strip() != "":
        ismp = int(line.strip())
        if ismp != 0:
          if ismp > maxval:
            nbig += 1
            if ismp > ismp_max: ismp_max = ismp
            ismp = maxval
          fval = float(ismp)/float(maxval)
          samples.append(fval)
    rd.close()
    if nbig > 0:
      sys.stderr.write(f"{nbig} samples above {maxval}, max = {ismp_max}\n")


# ######################################################################

def save_pixels(M, bands, page, vset):
  odir = f"pages/{page}/vsets/{vset}"
  bash(f"mkdir -p {odir}")
  ofile = f"{odir}/spectra.npy"
  numpy.save(ofile, M)
  bfile = f"odir/bands.txt"
  nb = len(bands)
  assert nb == numpy.size(M,1)
  with open(bfile, "w") as wr:
    for ib in range(nb):
      wr.write(f"{bands[ib]}\n")
    wr.close()
  return 
  # ----------------------------------------------------------------------
      
def load_pixels(page, clip,mask):
  ifile = f"deriv/davis/{page}/{clip}_{mktag}.npy"
  M = numpy.load(ifile)
  debug("shape(M)", numpy.shape(M))
  ns, nb = numpy.shape(M)
  bfile = f"deriv/davis/{page}/{clip}_{mktag}.bands"
  bands = []
  with open(bfile, "r") as rd:
    for line in rd:
      bands.append(line.strip())
  if len(bands) != nb: prog_error("M/bands mismatch")
  return M, bands