/* Last edited on 2024-11-06 06:08:08 by stolfi */

void spectrum_table_exact_append_all
  ( float_image_t *P, 
    bool_t center,
    int32_t c,
    spectrum_table_exact_t *tx, 
    bool_t verbose
  )
  { /* To reduce the cost of {spectrum_table_exact_sort}, we try to 
      generate the entries in approx sorting order, and pre-combine
      entries with the same absolute natural frequency. 
      
      Specifically, we enumerate the non-negative natural frequency
      vectors {fn[0],fn[1]} in order of increasing sum, and, for each
      vector, we collet and combine all terms that differ from those
      only in sign. This strategy enumerates the terms in order of
      increasing L1 norm of their integer frequency vectors {fn},
      where {fn[0]} is the number of wave cycles per image row, and
      {fn[1]} is the number of wave cycles per image column.
      
      Since {spectrum_table_exact_sort} uses the L2 metric on the
      fractional frequency vector {fn[0]/fd[0],fn[1]/fd[1]} (waves
      per pixel), this strategy is only an heuristic, even for square
      images, and becomes less effective as the aspect ratio moves
      away from 1:1. In any case, it usually saves 75% of the entries
      that would be used by the naive approach (one entry for each
      term of the spectrum). */

    int32_t cols = (int32_t)P->sz[1];
    int32_t rows = (int32_t)P->sz[2];
    int32_t fd[2] = { cols, rows };  /* Denominators of int32_t freq vectors. */
    int32_t fnMax[2] = { fd[0]/2, fd[1]/2 }; /* Max value of numerators {fn[0],fn[1]}. */
    int32_t s; /* Sum of numerators {fn[0]+fn[1]}. */
    int32_t sMax = fnMax[0] + fnMax[1];  /* Max value of {s}. */
    int32_t ntx = tx->ne; /* Table entries in use are {tx.e[0..ntx-1]}. */
    for (s = 0; s <= sMax; s++)
      { int32_t d; /* Difference of numerators {fn[0] - fn[1]}. */
        int32_t dMin = (int32_t)imax(s - 2*fnMax[1], -s);
        int32_t dMax = (int32_t)imin(2*fnMax[0] - s, +s);
        for (d = dMin; d <= dMax; d += 2)
          { /* Compute the numerators {fxp,fyp} from {s,d}: */
            assert((s + d) % 2 == 0);
            
            /* Main frequency vector corodinates: */
            int32_t fxp = (s + d)/2; assert((fxp >= 0) && (fxp <= fnMax[0]));
            int32_t fyp = (s - d)/2; assert((fyp >= 0) && (fyp <= fnMax[1]));
            
            /* Mirrored frequency vector corodinates: */
            int32_t fxn = (fd[0] - fxp) % fd[0]; /* {-fxp mod fd[0]} */
            int32_t fyn = (fd[1] - fyp) % fd[1]; /* {-fyp mod fd[1]} */ 
            
            /* Add the terms {fxn,fxp}×{fyn,fyp} but avoiding duplicates: */
            double nTerms = 0;
            double power = 0; 
            
            auto void add_term(int32_t fx, int32_t fy);
              /* Adds the Hartley spectrum term with indices {fx,fy} to {nTerms,power}. */
              
            void add_term(int32_t fx, int32_t fy)
              { int32_t rx = (center ? (fx + cols/2) % cols : fx);
                int32_t ry = (center ? (fy + rows/2) % rows : fy);

                double val = float_image_get_sample(P, c, rx, ry);
                if (verbose) { fprintf(stderr, "    adding %2d %4d %4d = %18.10f\n", c, fx, fy, val); }
                nTerms += 1;
                power += val;
              }
              
            /* Enumerate the four (or less) terms with freqs {±fn[0],±fn[1]}: */
            
            add_term(fxp, fyp);
            if (fxn != fxp) { add_term(fxn, fyp); }
            if (fyn != fyp) { add_term(fxp, fyn); }
            if ((fxn != fxp) && (fyn != fyp)) { add_term(fxn, fyn); }
            
            /* Append term to spectrum table: */
            int32_t fn[2] = { fxp, fyp };
            spectrum_table_exact_append_term(tx, &ntx, fn, fd, nTerms, power);
          }
      }

    spectrum_table_exact_trim(tx, ntx);
  }