#! /usr/bin/python3
# Last edited on 2026-04-08 13:51:19 by stolfilocal 

import sys, os, re
from sys import stderr as err
from math import sqrt, hypot, exp, log, sin, cos, floor, ceil
from math import isnan, isinf, inf, nan
import numpy as np

def main(img_file_in, img_file_bg, img_file_fg):
  img_in = read_image(img_file_in)
  
  rad_bg = 20.0 # Fuzzy window radius for bg computation.
  img_bg = compute_background_image(img_in, rad_bg)
  
  rad_fg = 5.0 # Fuzzy window radius for fg computation.
  img_fg = compute_foregound_image(img_in, img_bg, rad_fg)
  
  write_image(img_fg, img_file_ot
  return
  # ------------------------------------------------------------------------
  
def compute_background_image(img_in, rad_bg):
  ny, nx = np.shape(img_in)
  img_bg = np.zeros((ny, nx))
  for iy in range(ny):
    for ix in range(nx):
      img_bg[iy,ix] = compute_background_pixel(img_in, ix, iy, rad_bg)
  return img_bg
  # ----------------------------------------------------------------------
  
def compute_background_pixel(img_in, ix, iy, rad_bg):
  ny, nx = np.shape(img_in)
  points = [ ]
  hw = int(ceil(3*rad_bg))
  wt = weight_table(hw)
  for ky in range(-hw,hw+1):
    jy = iy + ky; 
    if jy >= 0 and jy < ny:
      for kx in range(-hw, hw_1):
        jx = ix + kx;
        if jx >= 0 and jx < nx:
          wxy = wt[abs(ky)]*wt[abs(kx)]
          points.append((kx, ky, img_in[ky,kx], wxy))
  coef = fit_quadratic_robust(points)
  return eval_quadratic(coef, 0, 0)
  # ----------------------------------------------------------------------

def it_quadratic_robust(points):
  X, F, W = get_lsq_data(points)
  tol = 1.0e-8
  coef = coef = fit_quadratic_lsq(X, F, W)
  change = +inf
  while change > tol
    prev_coef = coef
    W = remove_outliers(X, F, coef, W)
    coef = fit_quadratic_lsq(X, F, W)
    change = rn_dist(prev_coef, coef)
    
  return coef
  # ----------------------------------------------------------------------
  
if len(sys.argv) > 1:
  main(sys.argv[1], sys.argv[2], sys.argv[3])