#! /usr/bin/python3 # Last edited on 2025-08-23 16:22:49 by stolfi import sys, re, os, numpy from math import floor, inf from PIL import Image, PngImagePlugin from error_funcs import file_line_error, prog_error, arg_error from process_funcs import bash, run_command def read_mask_image(mfile): # Reads a mask image from file {mfile}. The result has shape # {(ny,nx)}, where {ny} and {nx} are the height and width of the # image, and pixel values either 0 or 255. M = numpy.asarray(Image.open(mfile)) Msh = numpy.shape(M) assert len(Msh) == 2 or (len(Msh) == 3 and Msh[2] == 1), f"mask {mfile} is not grayscale" # debug("M.shape", numpy.shape(M)) if len(Msh) == 3: # Reduce image from 3 to2 channels: M = M[:,:,0]; Msh = numpy.shape(M) # debug("M.shape", numpy.shape(M)) assert len(numpy.shape(M)) == 2 ny, nx = numpy.shape(M) smp_min = numpy.min(M, axis=(0,1)) smp_max = numpy.max(M, axis=(0,1)) smp_avg = numpy.mean(M, axis=(0,1)) sys.stderr.write(f"mask sample range = [{smp_min} _ {smp_max}] avg = {smp_avg}\n") # debug("M[0,0]", M[0,0]) return M # ---------------------------------------------------------------------- def read_clip_png_image(cfile, maxval): # Reads a clip image from file {cfile}. Converts the pixels from # integers to floats by dividing them by the specified {maxval}, without # any gamma correction, and clipping the result to {[0_1]}. The result # has shape {(ny,nx)}, where {ny} and {nx} are the height and width of # the image. # # Pixels greater than {maxval} are counted and a report is written to # {stderr} at the end. sys.stderr.write(f"reading {cfile} ...\n") P = numpy.asarray(Image.open(cfile)) # debug("P.shape", numpy.shape(P)) if len(numpy.shape(P)) == 3: # Reduce image from 3 to2 channels: nc = numpy.size(P,2) P = P[:,:,nc-1] # debug("P.shape", numpy.shape(P)) assert len(numpy.shape(P)) == 2 ny, nx = numpy.shape(P) # Convert pixels to floats and count outliers: C = numpy.zeros((ny, nx)) nbig = 0; true_vlo = 999999999; true_vhi = 0 for iy in range(ny): for ix in range(nx): ismp = P[iy,ix] assert ismp >= 0 true_vlo = min(true_vlo, ismp); true_vhi = max(true_vhi, ismp); if ismp > maxval: nbig += 1; ismp = maxval fval = float(ismp)/float(maxval) C[iy,ix] = max(0.0, min(1.0, fval)) if nbig > 0: sys.stderr.write(f"float range = {{{true_vlo:5d}..{true_vhi:5d}}}") fbig = float(nbig) / (nx*ny) sys.stderr.write(f" {nbig:5d} samples over {maxval:5d} ({100*fbig:6.2f}%)\n") return C # ---------------------------------------------------------------------- def read_clip_npy_images(page, clip, bands): # Reads the images of the clip {clip} in the spectral bands listed # in {bands}. The result is a {numpy} array {C} of floats with shape # {(ny,nx,nb)} with float samples in {[0 _ 1]}. C = None # For now. Later, the clips in array form. nb = len(bands) for ib in range(nb): band = bands[ib] cfile = f"pages/{page}/clips/{clip}/npy/{band}.npy" Ci = numpy.load(cfile) assert len(numpy.shape(Ci)) == 2, "clip image should be monochrome" if C is None: ny, nx = numpy.shape(Ci) C = numpy.zeros((ny,nx,nb)) else: assert numpy.shape(Ci) == (ny, nx), "bug inconsistent clip shape" C[:,:,ib] = Ci return C # ---------------------------------------------------------------------- def quantize_image(C, maxval): # Expects {C} to be a float array with samples in {[0_1]}. # Returns {C} converted to a matrix # of integers in {0..maxval} by linear scale. # # The data type of the resutlt will be {numpy.uint8} if {maxval} # is in {1..255}, and {numpy.uint16} if it is in {256..65535} assert isinstance(maxval,int) and maxval >= 1 and maxval <= 65535, \ f"bad {{maxval}} = {maxval}" sh = numpy.shape(C); dim = len(sh) if dim == 2: ny, nx = sh; nc = 1 fval_min = numpy.min(C, axis=(0,1)) fval_max = numpy.max(C, axis=(0,1)) elif dim == 3: ny, nx, nc = sh fval_min = numpy.min(C, axis=(0,1,2)) fval_max = numpy.max(C, axis=(0,1,2)) else: assert False, f"invalid image shape {sh}" dtype = numpy.uint8 if maxval <= 255 else numpy.uint16 P = numpy.zeros(sh, dtype=dtype) fval_min = +inf; fval_max = -inf ival_min = 9999999; ival_max = -1 nsma = 0; nbig = 0 for iy in range(ny): for ix in range(nx): for ic in range(nc): fval = C[iy,ix] if dim == 2 else C[iy,ix,ic] assert fval >= 0, f"invalid sample {fval:24.16e}" fval_min = min(fval_min, fval) fval_max = max(fval_max, fval) if fval < 0: nsma += 1; fval = 0.0 elif fval > 1.0: nbig += 1; fval = 1.0 ival = int(floor(maxval*fval + 0.5)) assert 0 <= ival and ival <= maxval, "bug quantization" ival_min = min(ival_min, ival) ival_max = max(ival_max, ival) if dim == 2: P[iy,ix] = ival else: P[iy,ix,ic] = ival sys.stderr.write(f"quantize_image: float range [{fval_min:8.6f} _ {fval_max:8.6f}]") sys.stderr.write(f" int range {{{ival_min:5d} .. {ival_max:5d}}}\n") if nsma > 0 or nbig > 0: sys.stderr.write(f" {nsma:6d} input samples below 0.00") sys.stderr.write(f" {nbig:6d} input samples above 1.00") return P # ---------------------------------------------------------------------- def write_image_as_gray_png_scaled(pfile, C, clo, chi, plo, phi): # Writes {C} as an grayscale linear PNG file, # mapping {[clo _ chi]} to {[plo _ phi]} and clipping # the result to {[0 _ 1]}. sys.stderr.write(f"writing {pfile} ...\n") ny, nx = numpy.shape(C) P = numpy.zeros((ny,nx)) nbig = 0; cval_max = -inf; pval_max = -inf nsma = 0; cval_min = +inf; pval_min = +inf nbad = 0; # Infinities or NANs. for iy in range(ny): for ix in range(nx): cval = C[iy,ix] if isfinite(cval): cval_min = min(cval_min, cval) cval_max = max(cval_max, cval) pval = plo + (cval - clo)*(phi - plo)/(chi - clo); pval_min = min(pval_min, pval) pval_max = max(pval_max, pval) if pval < 0.0: nsma += 1; pval = 0.0 elif pval > 1.0: nbig += 1; pval = 1.0 else: pval = 0.5 nbad += 1 P[iy,ix] = pval sys.stderr.write("input range [%+8.5f _ %+8.5f]" % (cval_min, cval_max)) sys.stderr.write(" mapped to [%+8.5f _ %+8.5f]\n" % (pval_min, pval_max)) if nsma > 0 or nbig > 0 or nbad > 0: sys.stderr.write(f"{nsma:6d} mapped samples were below 0.00\n") sys.stderr.write(f"{nbig:6d} mapped samples were above 1.00\n") sys.stderr.write(f"{nbad:6d} input samples were INF or NAN\n") write_image_as_gray_png(pfile, P) return # ---------------------------------------------------------------------- def write_image_as_gray_png(pfile, C): # Assumes that {C} is a grayscale image represented as a float # matrix with samples in {[0_1]}. Writes {C} to file # "{pfile}" as a linear grayscale PNG image file. sh = numpy.shape(C); dim = len(sh) if dim == 2: ny, nx = sh elif dim == 3: ny, nx, nc = sh assert nc == 1, f"invalid channel count {nc}" else: assert False, f"invalid image shape {sh}" maxval = 65535 sys.stderr.write(f"{pfile}: ") P = quantize_image(C, maxval) write_uint_array_as_png(pfile, P) return # ---------------------------------------------------------------------- def write_uint_array_as_png(pfile, P): # Expects {P} to be a {numpy} array with type {numpy.uint8} # or {numpy.uint16}. Writes it to file "{pfile}" as a # linear (grayscale or RGB) PNG image. # sh = numpy.shape(P); dim = len(sh) if dim == 2: ny, nx = sh; nc = 1 elif dim == 3: ny, nx, nc = sh assert nc >= 1 and nc <= 4, f"invalid channel count {nc}" else: assert False, f"invalid image shape {sh}" assert P.dtype == numpy.uint8 or P.dtype == numpy.uint16, \ "invalid matrix type" img = Image.fromarray(P, mode="I;16") img.save(pfile, pnginfo=PngImagePlugin.PngInfo()) return # ---------------------------------------------------------------------- def normalize_gray_image(G, M): # Assumes that {G} is a {numpy} array representng an image with # positive and negative samples and arbitrary mean, with shape # {(ny,nx,1)} # # The procedure scales and shifts all samples of G so that those # samples where the mask image {M} is nonzero span the range {[0_1]}, # apart from a small percentage. Mapped values (anywhere) that fall # outside that range are clipped to it. # # The image {M} must be an array with shape {(ny,nx)}. ny, nx, nc = numpy.shape(G) assert numpy.shape(M) == (ny, nx) # Get the : frin = 0.95 # Fraction of masked pixels to keep in range. smin, smax = get_channel_quantiles(G, 0, M, frin) sys.stderr.write(f"quantile range:") sys.stderr.write(f" [{smin:7.4f} _ {smax:7.4f}]") smid = (smin + smax)/2 srad = (smax - smin)/2 # Scale and shift each channel: nbig = 0 fsmp_min = +inf fsmp_max = -inf for iy in range(ny): for ix in range(nx): fsmp = (G[iy,ix,0] - smid)/srad fsmp = 0.5*(1 + fsmp) if fsmp < 0.0 or fsmp > 1.0: nbig += 1 fsmp_min = min(fsmp_min, fsmp) fsmp_max = max(fsmp_max, fsmp) fsmp = min(1.0, max(0.0, fsmp)) G[iy,ix,0] = fsmp sys.stderr.write(f"channel {0}: {nbig:4d} out of range") sys.stderr.write(f" {frin:0.3f} masked quantiles = ") sys.stderr.write(f" [{smin:7.4f} _ {smax:7.4f}]") sys.stderr.write(f" {nbig:4d} out of range") sys.stderr.write(f" actual range [{fsmp_min:7.4f} _ {fsmp_max:7.4f}]\n") return # ---------------------------------------------------------------------- def normalize_YUV_image(YUV, M): # Assumes that {YUV} is a {numpy} array representing image in the # linear YUV color space, with positive and negative samples and # arbitrary mean, with shape {(ny,nx,3)}. # # The procedure scales and shifts the samples in each channel, so # that they are normalized. The normalization leaves the Y channel # (channel 0) with mean 0.5 and deviation 0.25. while the U and V # channels (1 and 2) are left with zero mean and deviation 0.5. # # If {M} is not {None}, it must be a {numpy} array with # shape {(ny,nx)}. In this case the scaling and shifting will # be determined only by the channel {ic} samples of those pixels # where {M} is nonzero. # ny, nx, nc = numpy.shape(YUV) assert nc == 3, "image {YUV} msut have three channels" assert numpy.shape(M) == (ny, nx) # Scale and shift each channel: for ic in range(nc): cvals = [ \ YUV[iy,ix,ic] \ for ix in range(nx) \ for iy in range(ny) \ if M == None or M[iy,ix] != 0 \ ] savg = numpy.mean(cvals) sdev = numpy.std(cvals) # New {numpy} allows {mean = savg} arg. Yrad = 0.30 # Max Y variation from 0.5. vavg = 0.50 if ic == 0 else 0.00 # Desired mean. vdev = Yrad/2 if ic == 0 else 0.50 # Desired deviation. for iy in range(ny): for ix in range(nx): fsmp = vavg + vdev*(YUV[iy,ix,ic] - savg)/sdev if ic == 0: # Clip the {Y} value to a proper sub-range of {[0 _ 1]}: fsmp = min(0.5+Yrad, max(0.5-Yrad, fsmp)) YUV[iy,ix,ic] = fsmp # Rescale the U,V coordinates to fit in the RGB cube: sUV = choose_UV_scale(YUV, M) YUV[:,:,1] = sUV*YUV[:,:,1] YUV[:,:,2] = sUV*YUV[:,:,2] return # ---------------------------------------------------------------------- def get_channel_quantiles(G, ic, M, frin): # Assumes that {G} is a {numpy} array with shape {(ny,nx,nc)} # representing an image with {nc} channels, and {ic} is # a channel index in {0..nc-1}. # # Returns the tightest values {smin,smax} so that the fraction of the # samples of {G} that fall in the {[smin _ # smax]}, is at least {frin}, equal numbers falling below {smin} and # above {smax}. # # If {M} is not {None}, it must be a {numpy} array with shape # {(ny,nx)}. The quantiles above will be computed considering only # those pixels of {G} where {M} is nonzero. # ny, nx, nc = numpy.shape(G) assert numpy.shape(M) == (ny, nx), "bad mask shape" assert 0 <= ic and ic < nc, "bad {ic}" cvals = [ \ G[iy,ix,ic] \ for ix in range(nx) \ for iy in range(ny) \ if M == None or M[iy,ix] != 0 \ ] # Normalize to {95%} of samples spanning {[0_1]}, with clipping: assert 0.0 <= frin and frin <= 1.0, "bad {frin}" smin = numpy.quantile(cvals, 0.5 - frin/2) smax = numpy.quantile(cvals, 0.5 + frin/2) smid = (smin + smax)/2 srad = (smax - smin)/2 return smin, smax # ---------------------------------------------------------------------- def choose_UV_scale(YUV, M): # Assumes that {YUV} is a {numpy} array representing an image in the # linear YUV color space, with positive and negative samples and # arbitrary mean, with shape {(ny,nx,3)}. # # Chooses a scaling factor for the U and V coordinates of {YUV} so # that most of the pixels will be inside the {[0 _ 1]} cube when the # image is converted to RGB. # # Requires that the Y values in {YUV} be in the range {[0.2 _ 0.8]}. # # If {M} is not {None}, it must be a {numpy} array with shape # {(ny,nx)}. In this case, the scaing factor will be computed # considering only those pixels of {YUV} where {M} is nonzero. ny, nx, nc = numpy.shape(YUV) assert nc == 3, "image {YUV} msut have three channels" assert numpy.shape(M) == (ny, nx) svals = [] for ix in range(nx): for iy in range(ny): if M == None or M[iy,ix] != 0: Y, U, V = YUV[iy,ix,:] assert Y >= 0.19999 and Y <= 0.80001, "bug: Y range" RGB_UV = RGB_pixel_from_YUV_pixel((0,U,V)) sUV = 1.0e100 for kc in range(3): if RGB_UV[kc] > 1.0e-100: sUV = min (sUV, (1.0-Y)/RGB_UV[kc]) elif RGB_UV[kc] < -1.0e-100: sUV = min(sUV, (0.0-Y)/RGB_UV[kc]) svals.append(sUV) sUV_max = numpy.quantile(svals, 0.05) return sUV_max # ---------------------------------------------------------------------- def RGB_image_from_YUV_image(YUV): # Assumes that {YUV} a {numpy} array representing an image in YUV # color space, with shape {(ny,nx,3)}. Converts each pixel by linear # transformation from YUV space to linear RGB, still with positive and # negative samples. Clips the pixels to the unit cube. # ny, nx, nc = numpy.shape(YUV) assert nc == 3, "bug {nc}" nout = 0; RGB_min = +inf; RGB_max = -inf RGB = numpy.zeros((ny, nx, 3)) for iy in range(ny): for ix in range(nx): pix = YUV[iy,ix,:] pix = RGB_pixel_from_YUV_pixel(pix) for ic in range(3): if pix[ic] < 0: nout += 1; RGB_min = min(RGB_min, pix[ic]); pix[ic] = 0 elif pix[ic] > 1: nout += 1; RGB_max = max(RGB_max, pix[ic]); pix[ic] = 1 RGB[iy,ix,:] = pix if nout > 0: sys.stderr.write(f"RGB_image_from_YUV_image: {nout:5d} pixels outside cube\n") sys.stderr.write(f"Actual range [{RGB_min:7.3f} _ {RGB_max:7.3f}]\n") return RGB # ---------------------------------------------------------------------- def RGB_pixel_from_YUV_pixel(pix): Y, U, V = pix; R = +1.000000*Y -0.001164*U +1.139704*V; G = +1.000000*Y -0.395735*U -0.580624*V; B = +1.000000*Y +2.030875*U -0.000606*V; return [R, G, B] # ----------------------------------------------------------------------