#! /usr/bin/python3 # Last edited on 2025-10-15 11:50:50 by stolfi # Reads a set of multispectral images and removes from them the spectrum # of each pixel those components that could be due to variation in the # substrate (parchment) spectrum, as deermined in a previosu spectral # analysis of a substrate sample. Replaces them with a neutral ("gray") # spectrum of uniform brightness. Then writes the resulting # "background-cleaned" images, with all pixel spectra and scaled to fit # in the unist spectrum hypercube. # # Command line arguments are the {page} and {clip}, a tag {run_tag} for # this run of the program, the tag {mask_bg} of the mask used to determine # the spectral features of the substrate, and the tag {mask_re} of a # second mask that selects the most iportant pixels which are considered # when adjusting luminance and saturation of the background-cleaned # spectra. # # The pixels {mask_bg} must have been analyzed with # {analyze_clip_mask_pixel_spectra.py} and the result should be in the # file "pages/{page}/clips/{clip}/masks/{mask_bg}/analysis.txt". That # file specifies the number {nb} and names {bands[0:nb]} of the spectral # bands used in the analysis. The spectrum of the background is supposed # to vary over the 2D suspace of {\RR^{nd}} generated by the mean # spectrum {bg_avg} and the first two eigenvevtors {P[0:nb,0]} and # {P[0:nb,1]} found by the analysis. # # The original clip images are read from the files # "pages/{page}/clips/{clip}/npy/{bands[ib]}.npy", for {ib} in # {0..nb-1}, as a numpy array {P[0:ny,0:nx,0:nb]}. # # The background-cleaned images are written to files # "pages/{page}/clips/{clip}/rembgs/{run_tag}/nobgs/{bands[ib]}.png", for # {ib} in {0..nb-1}. # # Also writes to "pages/{page}/clips/{clip}/rembgs/{run_tag}/parms.txt" a # small file with the main parameters used in the run. # import sys, re, os from math import sqrt, sin, cos, log, exp, floor, pi, inf, isfinite import numpy from PIL import Image from error_funcs import prog_error, arg_error, debug from process_funcs import bash, run_command import argparser import vms_linear_gray_image_funcs as gfn import vms_multispectral_image_funcs as mfn import vms_spectral_analysis_funcs as afn def main(): page, clip, run_tag, mask_bg, mask_re = parse_options() # Read the spectral bands tables: pfile = f"MS/davis/{page}/bands.txt" sys.stderr.write(f"reading bands table {pfile} ...\n") bands_tb = mfn.read_page_bands_tables(pfile) # Read the spectral analysis results for the mask {mask_bg}: sys.stderr.write(f"reading analysis results for mask {mask_bg} ...\n") bands, bg_avg, bg_pc0, bg_pc1 = read_analysis_data(page, clip, mask_bg) nb = len(bands) debug("nb", nb) assert numpy.shape(bg_avg) == (nb,) assert numpy.shape(bg_pc0) == (nb,) assert numpy.shape(bg_pc1) == (nb,) # Read the clip images: sys.stderr.write(f"reading multispectral '.npy' clip images ...\n") C = mfn.read_clip_npy_images(page, clip, bands) ny, nx = numpy.shape(C[:,:,0]) debug("ny", ny) debug("nx", nx) assert numpy.shape(C) == (ny,nx,nb,) # Write false color image before bg removal: bdir_inp = make_output_directory(page, clip, run_tag, "inp", -1, -1) sys.stderr.write(f"writing initial false-color image to {bdir_inp}/ ...\n") write_false_color_image(bdir_inp, bands, bands_tb, C) # Read the important pixels mask image: Mfile = f"pages/{page}/clips/{clip}/masks/{mask_re}/mask.png" sys.stderr.write(f"reading mask {Mfile} ...\n") M = gfn.read_mask_image(Mfile) # Create background-less images: sys.stderr.write(f"removing background spectra ...\n") D = remove_background_from_images(C, bg_avg, bg_pc0, bg_pc1) # Add back the mean background color: mchr_man = 5.0 B = mfn.scale_shift_and_clip_images(D, mchr_man, bg_avg) # Write raw background-less images: bdir_raw = make_output_directory(page, clip, run_tag, "clean", -1, -1) sys.stderr.write(f"writing background-less images to {bdir_raw}/ ...\n") write_multispectral_clip_images(bdir_raw, bands, B) write_false_color_image(bdir_raw, bands, bands_tb, B) write_parameters(bdir_raw, mask_bg, mask_re, mchr_man, -1.0, -1.0) # Again, but with optimum chroma: sys.stderr.write(f"optimum rescaling and shifting ...\n") B, mchr_opt = mfn.scale_and_shift_to_fit_unit_cube(D, bg_avg, M) # Write chroma-stretched background-less images: bdir_str = make_output_directory(page, clip, run_tag, "stretched", -1, -1) sys.stderr.write(f"writing stretched background-less images to {bdir_str}/ ...\n") write_multispectral_clip_images(bdir_str, bands, B) write_false_color_image(bdir_str, bands, bands_tb, B) write_parameters(bdir_str, mask_bg, mask_re, mchr_opt, -1.0, -1.0) return 0 # Add the new substrate spectrum: sys.stderr.write(f"mixing in the new artificial substrate ...\n") luma_fg = 0.30 # Luminance for fully saturated pixels. luma_bg = 0.95 # Luminance for pure substrate pixels. B = mix_new_substrate(B, M, luma_fg, luma_bg) # No need to rescale chroma to fit cube since it is a convex op. # Write cleaned RGB image with uniform luma: bdir_mix = make_output_directory(page, clip, run_tag, "remixed", luma_fg, luma_bg) sys.stderr.write(f"writing artificial substrate images to {bdir_mix}/ ...\n") write_multispectral_clip_images(bdir_mix, bands, B) write_false_color_image(bdir_mix, bands, bands_tb, B) sys.stderr.write(f"writing parameter file to {bdir_mix}/ ...\n") write_parameters(bdir_mix, mask_bg, mask_re, mchr_opt, luma_fg, luma_bg) return 0 # ---------------------------------------------------------------------- def read_analysis_data(page, clip, mask): # Reads the bands list {bands} and the analysis data for the given # {page}, {clip}, and {mask}. Namely, reads and returns the list # {bands[0:nb]} of band images used in the analysis, the average # spectrum {bg_avg[0:nb]}, and the first two peincipal components # (eigenvectors) {bg_pc0=P[:,0]} and {bg_pc1=P[:,1]}. # nb = None # For now, reset later. ne = None # For now, reset later. bands, illums, wlens, bg_avg, smin, smax, rdev, tdev, P = \ afn.read_mask_analysis_data(page, clip, mask) bg_pc0 = P[:,0] bg_pc1 = P[:,1] return bands, bg_avg, bg_pc0, bg_pc1 # ---------------------------------------------------------------------- def remove_background_from_images(C, bg_avg, bg_pc0, bg_pc1): # Receives a multispectral clip image set {C}, a 2-dimensional affine # subspace of {\RR^{nd}} defined by the mean spectrum {bg_avg} and # principal components {bg_pc0,bg_pc1}. # # Subtracts the assumed mean color {bg_avg} from every pixel spectrum {C[iy,ix,:]} and projects # the result on the {nd-2}-dimensional linear supspace of {\RR^{nd}} that # is orthogonal to {bg_pc0} and {bg_pc1}. # ny, nx, nb = numpy.shape(C) debug("nb", nb) assert numpy.shape(bg_avg) == (nb,), "bug bad bg_avg shape" assert numpy.shape(bg_pc0) == (nb,), "bug bad bg_pc0 shape" assert numpy.shape(bg_pc1) == (nb,), "bug bad bg_pc1 shape" # Create the matrix {U} that projects {\RR^{nd}} on the subsspace # orthogonal to {bg_pc0} and {bg_pc1}: U = numpy.zeros((nb,2,)) U[:,0] = bg_pc0/numpy.linalg.norm(bg_pc0) U[:,1] = bg_pc1/numpy.linalg.norm(bg_pc1) P = numpy.identity(nb) - (U @ numpy.transpose(U)) assert numpy.shape(P) == (nb,nb,), "bug P shape" # Testing: apix = 0.7*bg_pc0 + 0.3*bg_pc1 bpix = apix @ P for pc in ( bg_pc0, bg_pc1,): dot = numpy.dot(bpix, pc) sys.stderr.write("test: dot(bpix, pc) = %+10.8f\n" % dot) assert abs(dot) <= 2.0e-6, "bug P projection 1" # Apply projection to each pixel: sys.stderr.write("applying projection to all pixels ...\n") D = numpy.zeros((ny,nx,nb,)) luma_min = +inf luma_max = -inf for iy in range(ny): for ix in range(nx): cpix = C[iy,ix,:] dpix = (cpix - bg_avg) @ P assert numpy.shape(dpix) == (nb,), "bug dpix shape" luma = numpy.sum(dpix)/nb luma_min = min(luma_min, luma) luma_max = max(luma_max, luma) D[iy,ix,:] = dpix sys.stderr.write(f"luma range = [{luma_min:.6f} _ {luma_max:.6f}]\n") for ib in range(nb): val_min = numpy.min(D[:,:,ib]) val_max = numpy.max(D[:,:,ib]) sys.stderr.write(f"band {ib:2} range = [{val_min:.6f} _ {val_max:.6f}]\n") return D # ---------------------------------------------------------------------- def mix_new_substrate(F, M, luma_fg, luma_bg): # Receives a set {F} of multispectral images with shape {(ny,nx,nb,)}. # ??? {M} # # Shifts every spectrum in {F}, preserving its chroma vector, so that # its luminance is {luma_fg}. The norm of the chroma vector must not # exceeed {mchr}. # # Then scales up each shifted spectrum {fpix} so that it becomes a # mixture of the new substrate spectrum {L = luma_bg*(1,1,...,1)} and # a spectrum with same hue, chroma norm {mchr}, and luminosity # {luma_fg}. # ny, nx, nb = numpy.shape(F) assert numpy.shape(M) == (ny,nx,) # project spectra onto cone defined by {msat,fg,bg}: sys.stderr.write("mixing neutral background spectrum ...\n") G = numpy.zeros((ny,nx,nb,)) gpix_fg = numpy.full((nb,), luma_fg) gpix_bg = numpy.full((nb,), luma_bg) for iy in range(ny): for ix in range(nx): fpix = F[iy,ix,:] flum = numpy.mean(fpix) fchr = fpix - numpy.full((nb,), flum) fchr_norm = numpy.linalg.norm(fchr) assert fchr_norm < 1.0001*mchr, "bug excessive chroma" fsat = fchr_norm/luma_fg # Compute background mixing coeffs {alfa,beta}: alfa = fsat/(fsat + (msat - fsat)*luma_fg/luma_bg) beta = 1 - alfa gpix = alfa*(gpix_fg + fchr*msat/fsat) + beta*gpix_bg G[iy,ix,:] = gpix # Checking: glum = numpy.mean(gpix) gchr = gpix - numpy.full((nb,), glum) gsat = numpy.linalg.norm(gchr)/glum assert abs(gsat - fsat) < 0.0001, "bug sat changed" return G # ---------------------------------------------------------------------- def make_output_directory(page, clip, run_tag, tag, luma0, luma1): # Constructs the name {bdir} of the output directory, namely # "pages/{page}/clips/{clip}/rembgs/{run_tag}/{tag}-{MM}-{NN}" # where {MM} and {NN} are the values {lum0} and {lum1} scaled # by 100 with format '%02d'. If either is {-1}, skips that part. # # Also creates that diectory. subdir = tag for luma in (luma0, luma1): if luma >= 0: subdir += ("-%02d" % int(100*luma + 0.5)) bdir = f"pages/{page}/clips/{clip}/rembgs/{run_tag}/{subdir}" sys.stderr.write(f"creating directory {bdir}/ ...\n") bash(f"mkdir -p {bdir}") return bdir # ---------------------------------------------------------------------- def write_multispectral_clip_images(bdir, bands, B): # Given a multispectral set {B} of images, with shape {(ny,nx,nb,)} # writes each image {B:,:,ib]} to disk with nams # "{bdir}/{band[ib]}.png" for {ib} in {0:nb}. # ny, nx, nb = numpy.shape(B) assert nb == len(bands), "bug nb" for ib in range(nb): band = bands[ib] bfile = f"{bdir}/{band}.png" sys.stderr.write(f"writing {bfile} ...\n") gfn.write_image_as_gray_png(bfile, B[:,:,ib]) return # ---------------------------------------------------------------------- def write_false_color_image(bdir, bands, bands_tb, B): # Given a multispectral set {B} of images, with shape {(ny,nx,nb,)} # writes a false RGB image "{bdir}/comp.png" derived from {B} here # the long, medium, and short wavelength images above are # mapped to the red,green, and blue channels, respectively. # ny, nx, nb = numpy.shape(B) assert nb == len(bands), "bug nb" sys.stderr.write(f"creating the pseudocolor image ...\n") RfromS = mfn.spectrum_to_color_matrix(bands, bands_tb) assert numpy.shape(RfromS) == (nb,3,) P = numpy.zeros((ny,nx,3,)) for iy in range(ny): for ix in range(nx): P[iy,ix,:] = B[iy,ix,:] @ RfromS # Write the false-color composite: pfile = f"{bdir}/comp.png" sys.stderr.write(f"writing {pfile} ...\n") P = gfn.quantize_image(P, 255) gfn.write_uint_array_as_png(pfile, P) dtag = re.sub(r"^.*[/]", "", bdir) bash(f"display -title '{dtag}/comp.png' {pfile}") return # ---------------------------------------------------------------------- def write_image(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") gfn.write_image_as_gray_png(pfile, P) return # ---------------------------------------------------------------------- def parse_options(): na = len(sys.argv) assert na >= 6, "insuff args" page = sys.argv[1] clip = sys.argv[2] # name of clip to take pixels from. run_tag = sys.argv[3] # Name of rembgration run. mask_bg = sys.argv[4] # Name of background mask. mask_re = sys.argv[5] # Name of relevant pixels mask. return page, clip, run_tag, mask_bg, mask_re # ---------------------------------------------------------------------- def write_parameters(bdir, mask_bg, mask_re, mchr, luma_fg, luma_bg): pfile = f"{bdir}/parms.txt" with open(pfile, "w") as wr: wr.write(f"mask_bg = {mask_bg}\n") wr.write(f"mask_re = {mask_re}\n") wr.write(f"mchr = {mchr:10.6}\n") wr.write(f"luma_fg = {luma_fg:8.6}\n") wr.write(f"luma_bg = {luma_bg:8.6}\n") # More parameters? wr.close() return # ---------------------------------------------------------------------- main()