/* Last edited on 2024-10-12 13:27:26 by stolfi */

    argparser_get_keyword(pp, "-winSize");
    o->winSize = (int32_t)argparser_get_next_int(pp, 3, 99);  
    if (o->winSize % 2 != 1) { argparser_error(pp, "window size must be odd"); }

void mfdo_read_term_index_table
  ( char *fname, 
    int32_t NB,
    char *belName[],
    int32_t *NP_P,
    multifok_term_prod_t **prix_P,
    int32_t *NT_P,
    char ***termName_P
  );
  /* Calls {multifok_term_read_index_table} on the file "{fname}". */

void mfdo_read_term_index_table
  ( char *fname, 
    int32_t NB,
    char *belName[],
    int32_t *NP_P,
    multifok_term_prod_t **prix_P,
    int32_t *NT_P,
    char ***termName_P
  )
  { bool_t verbose = TRUE;
    FILE *rd = open_read(fname, verbose);
    multifok_term_read_index_table(rd, NB, belName, NP_P, prix_P, NT_P, termName_P, verbose);
    fclose(rd);
  }

    fprintf(stderr, "ignored %d pixels for low sharpness\n", NL_blur); 

int32_t NL_blur = 0;    /* Number of pixels discarded because of low {sharp}. */
    
            double sharp_min = 0.2;      /* Ignore pixels with sharpness below this. */

 if (sharp < sharp_min)
              { NL_blur++;
                mkval = 0.25;
              }
            else

    
    o->term = string_vec_new(20);
    int32_t NT = 0; /* Number of terms to plot. */
    while (argparser_keyword_present(pp, "-term"))
      { char *tnk = argparser_get_next_non_keyword(pp);
        string_vec_expand(&(o->term), NT);
        o->term.e[NT] = tnk;
        NT++;
      }
    string_vec_trim(&(o->term), NT);
    if (NT == 0) { argparser_error(pp, "must specify at least one \"-term\""); }

    
    
    mfdo_frame_spec_vec_t frameSpec;            /* Specifications of the input image files. */
 
typedef struct mfdo_frame_spec_t
  { 
  } mfdo_frame_spec_t;
  /* Data that determines the name of each input image file. */
 
vec_typedef(mfdo_frame_spec_vec_t,mfdo_frame_spec_vec,mfdo_frame_spec_t);
    mfdo_frame_spec_vec_t *imo

vec_typeimpl(mfdo_frame_spec_vec_t,mfdo_frame_spec_vec,mfdo_frame_spec_t);

  ( int32_t NI,
    mfdo_frame_t *frame,
    int32_t NW,
    double ws[],
    double noise,
    int32_t NB,
    double **bas,
    int32_t NP,
    multifok_term_prod_t prix[],
    int32_t NT,
    char *outPrefix
  )
    FILE *wr_dat,
    char *inDir,
    char *outPrefix,

typedef struct mfdo_terms_t 
  { /* Basis of linear window operators: */
    int32_t NB;                  /* Number of basis elements. */
    double **bas;                /* Coefficients window pixels in each basis element. */
    char **belName;              /* Names of basis elements. */
    /* Set of quadratic terms derived from basis elements: */
    int32_t NT;                  /* Number of quadratic terms to analyze. */
    char **termName;             /* Term names. */
    int32_t NP;                  /* Number of basis coeff products in the quadratic terms. */
    multifok_term_prod_t *prix;  /* Mapping of basis element index pairs to terms. */
  } mfdo_terms_t ;
  /* Tables describing an orthogonal basis of linear window operators,
    and a set of quadratic window terms that may be
    combined into a focus estimator.
    
    The basis elements are {bas[0..NB-1][0..NS-1]}, and their names
    are {belName[0..NB-1]}.  Specifically, {bas[ip][ks]} is the 
    value of window sample {ks} in element {ib}.
    
    The {NT} quadratic terms are described by the table {prix[0..NP-1]}.
    Specifically, {prix[it]} specifies the indices {jb1,jb2} of two 
    basis elements, and the index {kt} of the quadratic term 
    to which that product is to be added. */

/* Auxiliary op basis and op term functions for {test_mfok_diff_ops.c}. */
/* Last edited on 2024-10-10 18:39:37 by stolfi */

#ifndef test_mfok_diff_ops_terms_H
#define test_mfok_diff_ops_terms_H

#define _GNU_SOURCE
#include <stdint.h>
#include <stdio.h>

#include <bool.h>

#include <multifok_basis.h>

mfdo_terms_t *mfdo_get_basis_and_terms(int32_t NW, double ws[], multifok_basis_type_t bType);
  /* First, creates an operator basis of the specified {bType}
    for an {NW} by {NW} window with with window sample weights {ws[0..NS-1]}, 
    where {NS=NW*NW}.  See {multifok_basis_make}.  The basis will be orthonormal.
    
    Also prints the basis and writes its names out to file with
    {mfdo_write_basis_elem_names(outPrefix,NP,belName)}. */

#endif
/* See {test_mfok_diff_ops_basis.h}. */
/* Last edited on 2024-10-10 13:07:03 by stolfi */

#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <string.h>
#include <assert.h>

#include <vec.h>
#include <wt_table.h>
#include <wt_table_binomial.h>
#include <affirm.h>
#include <bool.h>
#include <multifok_window.h>
  
#include <test_mfok_diff_ops_basis.h>

/* INTERNAL PROCS */

void mfdo_read_term_names(char *fname, int32_t *NT_P, char ***termName_P);
  /* Reads term names (formulas) from file {fname}.
    See {multifok_term_read_names} for the file format. */

void mfdo_read_term_index_table
  ( char *fname, 
    int32_t NB,
    char *belName[],
    int32_t *NP_P,
    multifok_term_prod_t **prix_P,
    int32_t *NT_P,
    char ***termName_P
  );
  /* Calls {multifok_term_read_index_table} on the file "{fname}". */

/* IMPLEMENTATIONS */

mfdo_basis_t *mfdo_get_basis_terms
  ( char *termSetFile,
    int32_t NW,
    double ws[],
    multifok_basis_type_t bType
  )
  {
    mfdo_basis_t *basis = talloc(1, mfdo_basis_t);
    
    /* Generate the op basis {bas[0..NB-1]}: */
    basis->bas = NULL;
    basis->belName = NULL;
    bool_t ortho = TRUE;
    multifok_basis_make(bType, NW, ws, ortho, &(basis->NB), &(basis->bas), &(basis->belName));
    fprintf(stderr, "obtained NB = %d independent basis elements\n", NB);
    multifok_basis_print(stderr, NW, NB, bas, belName);
    multifok_basis_ortho_check(stderr, NW, NB, bas);
    
    /* Get quadratic term descriptions: */
    int32_t NT;   /* Number of quadratic terms to analyze. */
    char **termName; /* Term names. */
    int32_t NP;   /* Number of coeff products in the quadratic terms. */
    multifok_term_prod_t *prix;    /* Mapping of basis element index pairs to terms. */
    multifok_test_read_term_weights_names_get_indices
      ( termSetFile,
        basis->NB, basis->belName,
        &(basis->NT), NULL, &(basis->termName),
        &(basis->NP), &(basis->prix),
        TRUE
      ); 
   
      
    return basis;
  }

typedef struct multifok_frame_t 
  { double zFoc;            /* Height of frame's focal plane. */
    double zDep;            /* Depth of focus. */
    float_image_t *csimg;   /* Sythetic image with simulated focus blurring, converted to grayscale. */
    float_image_t *azimg;   /* The average {Z}-coord of the scene within each pixel. */
    float_image_t *dzimg;   /* The deviation of that {Z} coordinate within each pixel. */
    float_image_t *shimg;   /* The actual sharpness of {csimg[k]} at each pixel, from the ray tracing. */
  } multifok_frame_t;
  /* The four images that comprise a frame of the multifocus stack. */
  
typedef struct multifok_stack_t 
  { /* Image dimensions: */
    int32_t NX;            /* Pixels per row in each frame. */
    int32_t NY;            /* Pixels per column in each frame. */
    /* Frame data: */
    int32_t NI;            /* Number of frames in stack. */ 
    multifok_frame_t *frame;   /* The {NI} frames. */
    /* Window weights and noise: */
    int32_t NW;            /* Window height and width (must be odd). */
    double *ws;            /* Weights of pixels in window, lienarized by rows. */
    double noise;          /* Assumed noise level in images. */
    /* Contrast normalization: */
    float_image_t *nrimg;  /* Constrast normalization factor per pixel. */
  } multifok_stack_t ;
  /* A stack of synthetic multifocal frames, {frame[0.NI-1]}, all with the same 
    dimensions ({NX} by {NY} pixels). 
    
    The window pixel weights are {ws[0..NS-1]} where {NS=NW*NW} and {NW=2*HW+1}. 
    The weight of pixel in column {ix} and row {iy} of the window is {ws[ix + NW*iy]}. */

    stack->frame = talloc(NI, multifok_frame_t);
    
    stack->NX = NX;
    stack->NY = NY;
    
    /* Get the window size and weight mask: */

    for (int32_t ki = 0; ki < NI; ki++)
      {
        multifok_frame_t *fri = &(stack->frame[ki]);
        
        fri->zDep = zDep[ki];
        fri->zFoc = zFoc[ki];
        
        /* Assemble the frame name: */
        char *inPrefix = NULL;
        asprintf(&inPrefix, "%s/st%s-%04dx%04d-%s/frame", inDir, sceneType, NX, NY, pattern);
        char *frameTag = NULL;
        asprintf(&frameTag, "-fd%05.2f-zf%05.2f", fri->zDep, fri->zFoc);

        /* Read the simulated scene view image: */
        int32_t NC_csimg = 3;
        float_image_t *csimg_RGB = multifok_image_scene_color_read(inPrefix, frameTag);
        float_image_check_size(csimg_RGB, NC_csimg, NX, NY);

        /* Convert to grayscale: */
        fri->csimg = float_image_new(1, NX, NY);
        float_image_map_channels_RGB_to_YUV(csimg_RGB, fri->csimg); /* Discard {U,V} channels. */
        float_image_free(csimg_RGB);
        
        double zMax = multifok_scene_ZMAX;

        /* Read the actual sharpness image: */
        fri->shimg = multifok_image_sharpness_read(inPrefix, frameTag);
        float_image_check_size(fri->shimg, 1, NX, NY);

        /* Read the scene {Z} average image: */
        fri->azimg = multifok_image_zave_read(inPrefix, frameTag, zMax);
        float_image_check_size(fri->azimg, 1, NX, NY);

        /* Read the scene {Z} deviation image: */
        fri->dzimg = multifok_image_zdev_read(inPrefix, frameTag, zMax);
        float_image_check_size(fri->dzimg, 1, NX, NY);

      } 
"  Before computing the basis coefficients, the samples in each" \
  " window of each frame are adjusted to have zero weighted average and gradient." \
  " Thus the basis coefficients and the quadratic terms are independent of local image brightness and gradient.\n" \
  "\n" \
  
    /* Prefix for outopt image files: */
    char *stackPrefix = NULL;
    asprintf(&stackPrefix, "%s-img-st%s-%04dx%04d-%s", outDir, imo->sceneType, NX, NY, imo->pattern);