# Last edited on 2025-12-21 01:00:21 by stolfi

          
  
  
def least_squares_system(nb, basis, nd, sample, weight):
   

  if M is not None: assert numpy.shape(M) == (ny,nx,), "mask/image shape mismatch"
  
  for iy in range(ny):
    for ix in range(nx):
      smps = []  # List of pairs {(dx, dy)}.
      for q in range(4):
        smps_q = collect_pixels_in_quadrant(M, ix, iy, q, rmax, nmin_q)
        smpsw_q = add_weights_to_displacements(smps_q, rmax)
        
def add_apodizing_weights_to_displacements(smps, rmax):
  # Takes a list of index displacement vectors {(dx,dy)}.
  # Appends to each pair a weight
  # that is 1 for the nearest pixels and decreases smoothly to zero 
  # at distance {rmax}.
  smpsw = []
  rmin = hypot(smps[0][0], smps[0][1]) 
  rmax = hypot(smps[-1][0], smps[-1][1]) 
  if len(smps > 0):
    # Get distance of nearest and farthest pixels.
    # Note: these may be equal.
    assert rmax > 0 and rmin <= rmax
    rmid = (rmin + rmax)/2;
    # Compute edge-smoothing weights for this quadrant:
    for (dx, dy) in smps:
      if rmin == rmax:
        w = 0
      else:
        r = hypot(dx,dy)
        assert r >= rmin and r <= rmax
        if r >= rmax:
          w = 0.0
        elif r <= rmid:
          w = 1.0
        else:
          t = (r - rmid)/(rmax - rmid)
          w = (1-t) + 2*t*(t-0.5)*(t-1)
      if w > 0:
        smpsw.append((dx,dy,w))
  return
  # ----------------------------------------------------------------------

def collect_pixels_in_quadrant(M, ix, iy, q, rmax, nmin_q):     
  # Enumerates the pixels in quadrant {q} around pixel {M[iy,ix]}
  ny, nx = numpy.shape(M)
  smps = []
  for dv in range(1,rmax+1):
    for du in range(rmax+1):
        dx = (+du, +dv, -du, -dv)[q]
        dy = (+dv, -du, -dv, +du)[q]
        jx = ix+dx; jy = iy+dy
        if jx >= 0 and jx < nx and jy >= 0 and jy < ny and M[jy,jx] != 0:
          smps.append((dx, dy))
  smps.sort(key = lambda d : hypot(d[0],d[1]))
  if len(smps) > nmin_q: smps = smps[:nmin_q]
  return smps
  # ----------------------------------------------------------------------
  
  created = "created on " + str(datetime.now(tz = timezone.utc)) + "html_report_funcs.html_subpage_link "

  imgfile = f"{clip_dir}/annotated.png"  # Image file.
  titfile = f"{clip_dir}/title.txt"      # File with title of figure.
  legfile = f"{clip_dir}/legend.txt"     # File with plain text caption.
  disfile = f"{clip_dir}/discussion.txt" # File with plain text discussion.

  # Obtain the figure's dimensions:
  imgsize = h.get_image_size(imgfile)
  
  title = h.get_text_from_file(titfile)
  title = h.protect_html(title)
  title = h.simple_markup(title)
  
  if thumb:
    thumb_width = 80*st['text_width']/100
  else:
    thumb_width = 0
    assert link_text != None, "must specify a link text" 
  link_to_sub = h.make_link(st, subfile, link_text, thumb_width, imgfile)

  # Create the subsidiary page in memory:
  st_sub = h.new_doc(title, "#eeffdd", 800)
  
  caption = h.get_text_from_file(legfile)
  caption = h.protect_html(caption)
  caption = h.simple_markup(caption)
  caption = h.make_parags(caption, align = "left", width = "80%")

  imgfile_sub = f"annotated.png"  # Name of image file rel to subpage.
  imgtag_sub = f"<img src=\"{imgfile_sub}\"/>"
  h.figure(st_sub, imgtag_sub, caption, centered = True)
  
  discussion = h.get_text_from_file(legfile)
  if discussion != "":
    discussion = h.protect_html(discussion)
    discussion = h.simple_markup(discussion)
    discussion = h.make_parags(discussion)
    h.section(st_sub, 2, "Discussion")
    h.parags(st_sub, discussion)

  # Write out the subpage:
  wr_sub = open(subfile, "w")
  h.output_doc(st_sub, wr_sub, 0, created)
  wr_sub.close()

  # Scale each spectral component to standard deviation 1.0,
  # to give it equal chances of contributing to the result:

  # sys.stderr.error(f"scaling band images to unit deviation ...\n")
  # C = normalize_band_intensities_to_unit_std(C)
  


  # tfile = f"/tmp/{os.getpid()}-read_mask_image.pnm"
  # bash(f"convert {fname} -set colorspace LinearGray -colorspace LinearGray {tfile}")
  # M = read_pnm_image(tfile, 0, 0)