text: #include "math.h"
  #include "memory.h"
  #include "pthread.h"
  #include "stdarg.h"
  #include "stdint.h"
  #include "stdio.h"
  #include "stdlib.h"
  #include "unistd.h"
  #include "c_common.h"
  #include "fft.h"
  #include "_sigproc_v3.h"
  
  
  /*
  
  During calibration:
      We assume:
       - that all channels are rotationally aligned within tolerance
          What that tolerance is is not clear yet.
       - We do NOT assume that channels are translationally aligned.
       - The PSF has some axis of symmetry (like an Elliptical Gaussian)
  
      We obtain:
       - RegPSF per channel
          These regions are the same as the sub-division we will make of the field
          Compute FFT for these psfs along with scale factors.
       - RegIllum per channel
          Which includes balance, bias, threshold determined by sweeping
          over these parameters to maximize recall
       - ChannelAlignment
          A sub-pixel translation constant for each channel relative to channel 0
  
  During analysis:
      Cycle alignment:
          if cycle == 0:
              for several small "alignment regions" across the channel 0 image:
                  FFT those regions and store
          else:
              for several small "alignment regions" across the channel 0 image:
                  FFT region and multiply by the FFT from cycle 0: result is 1 pix alignment
                  Gradient descent to sub-pixel align
              if all alignment regions have approx same answer, then take the mean and proceed
              else resort to a slow "low density" aligner where the size of the
              "alignment regions" are much larger.
  
      For each regional subdivision:
          if cycle == 0:
              multichannel_peak_find using regillum parameters from calib
                  This uses 1-pixel channel alignment and binary masks
                  so that any peak in any channel is found
                  Store
  
          for each peak (cached from cycle 0):
              for each channel:
                  extract a peak_im (small, ~ 16x16 pixels)
                      (Offset by the 1-pix channel alignment)
  
                  fft2 peak_im
  
                  if cycle == 0:
                      gradient-descend to sub-pixel center and save
                      (What if there is no signal in this channel for this peak?!)
  
                  add in:
                      the sub-pix channel alignment
                      the sub-pix cycle alignment (zero if cycle == 0)
                      the sub-pix peak alignment
  
                  eval signal
                  eval noise (or similar for collision detection)
                  store sub-pix x, y, sig, noise
  
  Naming conventions:
      f_: A Fourier transformed variable
      pix_: A one-pixel location of some kind (may be either float or int but has no fractional part)
      spx_: A sub-pixel location, always float
      ch_: An per-channel array
      reg_: some sort of regional variable
  
      (Field, channel, cycle are always in that order when an order is required)
      fl_i: Field index
      ch_i: Channel index
      cy_i: Cycle index
  
      im_mea: The width/height of the full frame image (typically 1024+)
      reg_mea: The width/height of the a regional sub-image of im
      peak_mea: The width/height of a peak which is the size of the PSF plus some margin.
  
      ch_ims: An array (n_channels, im_mea, im_mea)
      im: One frame of ch_ims (im_mea, im_mea).
      reg_im: Some (reg_mea, reg_mea) sub-image of an im.
      peak_im: Some (peak_mea, peak_mea) sub-image of a reg_im
  */
  
  
  char *peak_find_multichannel(SigprocV3Context *ctx, Index reg_i) {
      /*
      Combines channels with binary thresholds, flood-fill and get approximate centers.
  
      Context requires:
          ch_ims: (n_channels, im_mea, im_mea)
          ch_aln: Offsets of each channel relative to channel 0
          ch_bal_im: int(im_mea, im_mea)
          ch_bias: int
          ch_thresh: int
          out_warning_n_locs_overflow: bool
  
      Outputs:
          By adding rows to out_sublocs with reg_i
  
      TODO:
          Angela has a feature request that this return a
          flag so we know WHICH channels peaks were found in
          In the current implementation I "or" the per-channel
          masks together before returning.
      */
      ensure_only_in_debug(ctx->sub_locs_tab.base != NULL, "Bad out_sublocs");
      const Size n_channels = ctx->n_channels;
      Index ch_i;
  
      const Size im_mea = ctx->im_mea;
      Size mask_shape[2] = { im_mea, im_mea };
      U8Arr msk_arr = u8arr_malloc(2, mask_shape);
      Uint8 *msk = u8arr_ptr1(&msk_arr, 0);
  
      Uint16 *ch_ims[N_MAX_CHANNELS];
      Uint16 *ch_bal_ims[N_MAX_CHANNELS];
      Sint64 ch_bias[N_MAX_CHANNELS];
      Sint64 ch_thresh[N_MAX_CHANNELS];
      check_and_return(n_channels < N_MAX_CHANNELS, "n_channels > N_MAX_CHANNELS");
      for(ch_i=0; ch_i<n_channels; ch_i++) {
          ch_ims[ch_i] = u16arr_ptr1(&ctx->ch_ims, ch_i);
          ch_bal_ims[ch_i] = u16arr_ptr1(&ctx->ch_bal_ims, ch_i);
          ch_bias[ch_i] = (Sint64)*u64arr_ptr1(&ctx->ch_bias, ch_i);
          ch_thresh[ch_i] = (Sint64)*u64arr_ptr1(&ctx->ch_thresh, ch_i);
      }
  
      ensure(is_pow_2(im_mea), "im_mea must be a power of 2");
      Uint64 y_shift = pow_2_to_shift(im_mea);
  
  	Index stack_top = 0;
  	const Size stack_size = 256;
  	const Size stack_size_minus_4 = 256 - 4;
  	Uint32 stack[stack_size];
  	#define PUSH(a,b) stack[stack_top] = ((a & 0xFFFF) << 16) | (b & 0xFFFF); stack_top++;
  	#define POP(a,b) stack_top--; a=(stack[stack_top] & 0xFFFF0000) >> 16; b=(stack[stack_top] & 0xFFFF);
  
      Index off = 0;
      Size hits = 0;
  
  	#define CMP(_off) \
  	    hits = 0; \
  	    for(ch_i=0; ch_i < n_channels; ch_i++) { \
  	        hits += (Sint64)ch_bal_ims[ch_i][_off] * ( (Sint64)ch_ims[ch_i][_off] - ch_bias[ch_i] ) > ch_thresh[ch_i]; \
  	    } \
  
      for(Index y=0; y<im_mea; y++) {
          for(Index x=0; x<im_mea; x++) {
              if(msk[off] == 0) {
                  CMP(off);  // Sets hits
  
                  // debug
                  // *u64arr_ptr2(&ctx->out_debug, y, x) = hits > 0 ? 1 : 0;
  
                  if(hits > 0) {
                      // Seed found, flood
                      Uint64 accum_loc_y = 0;
                      Uint64 accum_loc_x = 0;
                      Uint64 loc_n_pixels = 0;
  
                      PUSH(y, x);
                      while(stack_top > 0) {
                          // POP and EVALUATE the predicate of the pixel on the stack
                          Uint16 _y, _x;
                          POP(_y, _x);
                          if(0 <= _x && _x < im_mea && 0 <= _y && _y < im_mea) {
                              Uint8 *_msk = u8arr_ptr2(&msk_arr, _y, _x);
                              if(_msk[0] == 0 && stack_top < stack_size_minus_4) {
                                  CMP( (_y << y_shift) + _x );
                                  if(hits > 0) {
                                      *u64arr_ptr2(&ctx->out_debug, _y, _x) = 1;
                                      accum_loc_y += _y;
                                      accum_loc_x += _x;
                                      loc_n_pixels += 1;
  
                                      _msk[0] = 1;
                                      PUSH(_y+1, _x);
                                      PUSH(_y-1, _x);
                                      PUSH(_y, _x+1);
                                      PUSH(_y, _x-1);
                                  }
                              }
                          }
                      }
  
                      Size n_rows_free = tab_n_rows_remain(&ctx->sub_locs_tab);
                      if(n_rows_free > 0) {
                          const Float64 loc_n_pixels_f = (Float64)loc_n_pixels;
                          SubLoc sub_loc;
                          sub_loc.reg_i = reg_i;
                          sub_loc.sub_y = (Float64)accum_loc_y / loc_n_pixels_f;
                          sub_loc.sub_x = (Float64)accum_loc_x / loc_n_pixels_f;
                          tab_add(&ctx->sub_locs_tab, &sub_loc, NULL);
                      }
                      else {
                          Uint64 *out_warning = u64arr_ptr1(&ctx->out_warning_n_locs_overflow, 0);
                          *out_warning = 1;
                          goto break_scan;
                      }
                  }
              }
              off++;
          }
      }
  
      break_scan:
  
      u8arr_free(&msk_arr);
  
      return NULL;
  }
  
  char *context_init(SigprocV3Context *ctx) {
      return NULL;
  }
  
  char *context_free(SigprocV3Context *ctx) {
      return NULL;
  }
  
  // TODO. Lots of work here!
  