proj/mag/harp/libs/matlab/mfile-mex/hmi-mask-patch/smoothsphere.c

#include "mex.h"  /* must appear first */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <pthread.h>
#include <sys/time.h> // for timing
#include "mexhead.h"
#include "Doc/smoothsphere_docstring.h"  /* autogenerated from this file */

/**************************************************************

%smoothsphere: smooth a projected sphere with a kernel
% 
% [y,s,p]=smoothsphere(x,geom,k,kparam,kwt,bws)
% * Given a solar image of location given by `geom', smooth it 
% using the radially-symmetric kernel `k'.
% * This version is threaded, using 8 threads by default.
% Set MXT_NUM_THREADS in the environment to lower this number.
% * Off-disk values are indicated by NaN in the output.
% * The kernel used is dependent only on the weighted distance 
% between two positions, say P1 and P2, in three-dimensional 
% coordinates, normalized to live on the unit sphere:
%   d = P1 - P2
%   dist = d' W d  (>= 0)
% for a diagonal weight matrix W.
% The kernel values in `k' are specified for values of `dist'
% in a linear range from 0 to kparam(1), typically 0.015.  
% Given a certain value of `dist' found between an image pixel 
% and the kernel center, the associated weight is interpolated 
% linearly using the table.  If dot > kparam(1), a zero weight 
% is used.
% * We have typically let k be a gaussian kernel of width
% (standard deviation) = 0.0325, i.e., 
%   k(dist) = C * exp[ -0.5 * dist / sqr(0.0325) ], 
% since dist is a squared quantity already, and C is a constant
% chosen to make k have unit norm as a spherical kernel.
% * As a shortcut, if k is a scalar, the kernel is taken to 
% be the above function, but with width 0.0325 multiplied by
% k(1).  In this case, the LUT consists of 256 points evenly
% spaced over [0,kparam(1)].
% * The diagonal portion of the distance weight matrix W may be 
% supplied as a triple kwt.  This weights distances between 
% P1 and P2 in the x, y, and z directions respectively.   The
% P-angle (rotation in the x-y plane) is taken into account in
% this weighting, so that y is cross-track, and x and z are
% always along-track.
% * Typically the kernel falls to zero rather quickly.  As a
% computational shortcut, it is assumed that the kernel extends 
% only kparam(2) pixels on each side of its center.  For W = I,
% this implies that an "on" pixel in x extends to influence at 
% most a "swath" 2*kparam(2)+1 pixels on a side.  Both the P-angle 
% and the W matrix are taken into account in finding the swath.
% For instance, a W_y > 1 will cause the y portion of the swath
% to shrink, because distance decreases faster.  The swath is 
% always parallel to the (i,j) image axes, so P-angles not a 
% multiple of 90 degrees will cause the swath to be enlarged
% to contain a rotated rectangle.
% * To decrease computational load, nearby pixels may be grouped
% into meta-pixels called blocks.  This is especially trouble-free
% away from the limb, where local geometry is nearly planar, so
% the blocking is given as a table that depends on z (in [0,1]).
% A row in the table of (z,bw) means: above the value z, use
% a blocking of bw (>1).  For z=0 to bws(1,1), no blocking is used.
% * For example, the small-image default means to use single pixels 
% below 0.4, 2x2 blocks above 0.4, and 4x4 above 0.6.  For typical 
% MDI images, this is 16%, 20%, and 64% of on-disk pixels, respectively.
% Using block sizes that do not pack together is legal, but causes 
% inefficiency due to poor fits between abutting blocks; diagnose
% bad packing with the p output.
% * Supply bws=[] to treat all pixels singly (turn off blocking).
% * The optional output s is the z-coordinate (in [0,R]).
% The optional p output is the block (patch) number of each
% pixel.
% 
% Inputs:
%   real x[m,n];
%   real geom[5];  -- [x0 y0 rsun b0 p0]
%   real k[p] or k[1];
%   real kparam[2];  -- [top window]
%   opt real kwt = [1 1 1];
%   opt real bws[bwnum,2] = [0.4 2;0.6 4];         -- m <  2048
%                         = [0.3 2;0.5 4; 0.7 8];  -- m >= 2048
% 
% Outputs:
%   real y[m,n];
%   opt real s[m,n];
%   opt int  p[m,n];
% 
% See also: 

% implemented as a mex file

****************************************************************/

/****************************************************************
 Notes October 2010:

 * Pixel ordering.  I have started using a "mode" input that can be given
 as SESW for HMI and SENE for the old, transposed-MDI image ordering.  See
 roi_stats_mag for how this works.  Doing this allows this one code to be
 easily used for both cases.  It is implemented by adding a stride
 parameter for both x and y to the underlying computational routines.  The
 stride, and the x/y loop boundaries, are switched depending on "mode".
 For now, it's easier to just leave it alone; the user has to monkey with
 x0 and y0 to get the right behavior.

 * Beta.  I pass in beta-angle, but I don't use it.  I should.
 Perceived computational burden, as well as lack of technical clarity,
 has been the issue.  
 Technical clarity can come from the following observation.  
 Consider two points, P1 (the center point of the convolution) and P2 
 (the point where the kernel is evaluated).  Now, the distance is a 
 function of:

   D = P1 - P2 (in 3d coordinates: x, y, and z)

 In 3D coordinates, it is hard to account for beta, while still performing
 the right weighting according to kwt.  One way out is to realize there are
 two angles (distances along the sphere) that are important: longitudinal
 and latitudinal.  You can always talk about the offset of P2 from P1 in
 terms of a delta-lat and a delta-lon.  The kwt penalty should have two
 terms in this case, the penalty for a given delta-lat and that for a given
 delta-lon.  (We don't measure delta-lat and delta-lon in degrees, because
 the polar singularity would cause delta-lon to be unstable.)

 One way to find delta-lat and delta-lon is to inner-product the 
 difference D with an along-track unit vector and a cross-track unit
 vector.  (These are per-pixel triplets, found once at P1 and used for all
 P2.)  It might be best to just use a cross-track vector, and then 
 dump the entire residual into the along-track vector, so that the entire
 difference D is penalized using either the along-track or the cross-track
 weighting.  (I.e., D will not lie in the space spanned by just the
 along-track and cross-track vectors; there will be a third degree of 
 freedom because the planar approximation is not exact.)

 Depending on how it was done, this could correspond to a chord 
 length or an arc length.  It probably doesn't much matter which, because
 the kernel is a merely a computational device to obtain smoothing, and
 a simple scaling rule relates the two.

 This does not seem that burdensome computationally; it's about the same
 amount of computation in the inner loop as now, just a handful of multiplies
 and adds.  It should be prototyped to see how the results look for high beta.

****************************************************************/

/* constants and globals used for error checking */

#define NARGIN_MIN  4      /* min number of inputs */
#define NARGIN_MAX  6      /* max number of inputs */
#define NARGOUT_MIN 0      /* min number of outputs */
#define NARGOUT_MAX 3      /* max number of outputs */

#define ARG_X      0
#define ARG_GEOM   1
#define ARG_K      2
#define ARG_KParam 3
#define ARG_KWt    4
#define ARG_BWs    5

#define ARG_Y      0
#define ARG_S      1
#define ARG_P      2

static const char *progname = "smoothsphere";
#define PROGNAME smoothsphere
static const char *in_specs[NARGIN_MAX] = {
  "RM",
  "RV(5)",
  "RV",
  "RV(2)",
  "RV(0)|RV(3)",
  "RM"};
static const char *in_names[NARGIN_MAX] = {
  "x",
  "geom",
  "k",
  "kparam",
  "kwt",
  "bws"};
static const char *out_names[NARGOUT_MAX] = {
  "y", "s", "p"};


// short form of PROGNAME for #include identifiers
#define SHORTNAME ssp

#define KWt_count 3
#define Default_KWt_count 3
static double Default_KWt[Default_KWt_count] = { 1.0, 1.0, 1.0 };
#define Default_BWs_count_small 4
static double Default_BWs_small[Default_BWs_count_small] = { 0.4, 0.6, 2.0, 4.0 };
#define Default_BWs_count_large 6
static double Default_BWs_large[Default_BWs_count_large] = { 0.3, 0.5, 0.7, 2.0, 4.0, 8.0 };

/********************************************************************************
 *
 * Utility for making kernels
 *
 ********************************************************************************/

/*
 * make_kernel_lut: make a lookup table for a gaussian kernel
 * lut has Nker points, from 0..hi_ker (hi_ker typically = 0.015)
 * width of the kernel is kwid, which is given in scaled units,
 * so that kwid is order unity.
 * The lut, of length Nker, must be already allocated.
 */
static
int
make_kernel_lut(double *ker, int Nker, double hi_ker, double kwid)
{
  const double del = hi_ker/(Nker-1); // Nker fenceposts, Nker-1 gaps
  double x;
  double scalefac;
  double sum;
  double norm;
  int i;

  if (kwid <= 0   || isnan(kwid) || 
      hi_ker <= 0 || isnan(hi_ker) ||
      Nker <= 1   || ker == NULL)
    return 0;
  kwid *= 0.0325; // = the rescaled kernel width
  scalefac = 1/(2*kwid*kwid); // = the multiplier in the exp() call
  // fill in exp(- x / (2*kwid*kwid)) = exp(- x * scalefac)
  for (sum = i = 0; i < Nker; i++) {
    x = i*del;
    ker[i] = exp(- x * scalefac);
    sum += ker[i];
  }
  // find the normalization factor for unit mass
  // this should perhaps use del, but we use hi_ker/Nker because
  // that is what we have historically used.
  norm = M_PI * (hi_ker/Nker) * sum;
  scalefac = 1/norm;
  // re-normalize ker
  for (i = 0; i < Nker; i++)
    ker[i] *= scalefac;
  return 1;
}


/********************************************************************************
 *
 * Utilities for blocked convolution
 *
 ********************************************************************************/

/*
 * Idea: Divide source pixels into classes.  The class is a blocknum_t (like an
 * int) and it is a tag for a "meta-pixel" that subsumes a whole block of 
 * source-image pixels.  Each block is summarized by a blockstats_t, and is 
 * convolved all at once.
 *
 * Blocking is done because the main source of computational cost in this 2d convolution
 * is large AR's.  Formerly, each pixel in a large AR was mapped individually to the
 * output.  If the kernel width is KW, this one pixel results in KW^2 updates to the
 * output, each of which requires many flops for computation of a distance, and 
 * kernel lookup of the distance.  (The issue is *not* tiny one or two-pixel AR's; this
 * has been checked on real data; the typical AR pixel has 5-6 neighbors.)  Much
 * of this computation was wasted, because the smoothed activity is thresholded.
 *
 * To cut this time, we collect AR pixels that are relatively far from the limb into
 * blocks of size BWxBW, which are handled together.  The scheme below is more flexible
 * than this; it allows any blocking scheme (not just uniform square blocks).  Any mapping
 * of pixels to blocks is OK.  To change it, see block_create.
 *
 * The scheme assigns each pixel a blocknum_t.  If this is <0, it indicates the
 * pixel is handled individually.  If >0, it is a block number; all pixels with
 * that block number are handled in a single convolution inner loop over KW*KW
 * pixels.  The blocks are assigned by one routine, and then another loop over
 * the image computes the overall weight (sum of the weights in the block), and 
 * the net position (weighted by the convolved input).  The summary of each block
 * is put into a blockstat_t structure.  It would even be possible to combine 
 * blocks by manipulating just the blockstat_t summaries.
 *
 * Fundamentally, because considerable processing is devoted to each nonzero, ondisk
 * pixel, considerable intelligence can be used to optimize blocks.  Note that no
 * blocks are created for singletons (they are handled individually), and that
 * each block must contain some activity; the whole away-from-limb disk is not blocked
 * in its entirety.  It's typical to have a few thousand blocks (1024^2 image).
 * The current scheme uses blocks of width BW that are aligned on multiples of BW
 * in the input image, this is more efficient than fully ad-hoc block placement,
 * which tends to cause interference.
 *
 * We tried simpler combination schemes, but it turns out to be hard that every 
 * value in the source is being used exactly once, because of the interaction of
 * the blocks, the limb, and the offdisk stuff.
 *
 * In the final convolution loop, we look at the blocknum_t of each pixel.
 * If it is a block (and not a singleton), we look up its block in the list,
 * and use its average position and cumulative weight as a "metapixel".  We
 * also mark the block as used.  When we come to next pixel that is a member 
 * of that block, we do not count it again.
 */

// legal block values: 
// offdisk, skip, singleton, clear; or anything from 1 to #blocks
// (thus, any negative values are OK for the explicit block constants)
//  offdisk = outside visible disk
//  skip = visible, but zero weight.  skip during convolution
//  singleton = visible, nonzero weight, treated singly during the convolution
//  clear = un-initialized.  after block_create runs, *all* clears are gone
//  first, first+1, ... = part of the given block
// Thus, block 0 is un-used.  The valid blocks are 1..BlockNum, inclusive
typedef enum {
  Block_Min       = -3,      // always holds the minimum value, for external saves
  Block_Offdisk   = -3,      // not on the visible disk
  Block_Skip      = -2,      // on-disk, zero weight (skippable)
  Block_Singleton = -1,      // on-disk, nonzero weight, treated individually
  Block_Clear     = 0,       // un-initialized
  Block_First     = 1,       // this and following numbers: members of a given block
  Block_MAX       = INT_MAX  // ensures it's as big as an int
} blocknum_t;
// is the given block number part of a group
#define BLOCK_GROUPED(b) (b>=Block_First)

/*
 * statistics summarizing a block: its mass-weighted location (x,y,z),
 * its weight (wt>0), and its function value (f).  If the function
 * value is zero, the block is effectively turned off.  This is done
 * during the iteration once the block has been smoothed into the output.
 */
typedef struct {
  double x;  // location
  double y;
  double z;
  double wt; // weight (sum of abs(f))
  double f;  // function value
} blockstat_t;


/*
 * block_lut_validate: determine if a bwlen x 2 lookup table
 * supplied as an input argument is valid.  returns an
 * error message.  LUT order is:
 *
 *    tau1 block1
 *    tau2 block2
 *    tau3 block3
 *
 * There is an implicit 0 above the tau's, and an implicit 1 below
 * them.  For 0..tau1, no blocking.  For tau1..tau2, block1.  Etc.
 * Above tau3 and up to 1, block3.
 */
static
char *
block_lut_validate(double *bws, int bwlen)
{
  int i;

  for (i = 0; i < bwlen; i++) {
    // bin boundaries
    if (bws[i] < 0 || bws[i] > 1)
      return "bin boundaries must be in [0,1]";
    if (i > 0 && (bws[i-1] >= bws[i]))
      return "bin boundaries must be monotone increasing";
    // block sizes
    if (bws[bwlen+i] < 2)
      return "block sizes must be at least 2";
    if (((int) bws[bwlen+i]) != bws[bwlen+i])
      return "block sizes must be integers";
  }
  return NULL;
}

/*
 * block_lut_reformat: alter the lut as supplied slightly
 * so we can more readily use it.
 */
static
void
block_lut_reformat(double *tau1, double *taus, int *bws, 
           double *bws_in, int bwlen, double R)
{
  int i;

  // empty case
  if (bwlen == 0) {
    *tau1 = R + 1; // always use singleton rule
    return;
  }
  // non-empty
  *tau1 = bws_in[0] * R;
  // bin boundaries
  for (i = 1; i < bwlen; i++)
    taus[i-1] = bws_in[i] * R;
  taus[bwlen-1] = R + 1; // last bin > largest z
  // block sizes
  for (i = 0; i < bwlen; i++)
    bws[i] = (int) bws_in[bwlen+i];
}

/* 
 * block_create: make matrix of block-tags from (z,tau), the input map,
 * and a blockwidth parameter.  
 * Returns the number of distinct tags needed, or -1
 * for error (should never happen).
 */
static
blocknum_t
block_create(blocknum_t *blocks, double *s, double *a, 
         double tau1, double *taus, int *bws, int m, int n)
{
  int k;
  int x1, y1;
  int x2, y2;
  int xlo, xhi, ylo, yhi;     // loop bounds
  blocknum_t bn = Block_First;
  double s1;
  int bw;

  /* loop over image pixels */
  for (x1 = 0; x1 < n; x1++)
    for (y1 = 0; y1 < m; y1++) {
      // ensure blocks(x1,y1) is assigned a block number
      s1 = s[m*x1+y1];
      if (isnan(s1)) {
    blocks[m*x1+y1] = Block_Offdisk;
      } else if (a[m*x1+y1] == 0) {
    blocks[m*x1+y1] = Block_Skip; // zero contribution
      } else if (s1 < tau1) {
    // note, strict < allows tau1 = 0 to suppress all singletons
    blocks[m*x1+y1] = Block_Singleton; // singleton
      } else if (blocks[m*x1+y1] == Block_Clear) {
    // This code block fully controls the discretization
    // scheme.  Another scheme would just fill in the
    // block number "bn" in a different pattern, being sure
    // to not over-write blocks that are already assigned.
    //
    // determine the block width using s1 and the LUT
    // - if the LUT was empty, tau1 will be > R, and
    // this case will not be reached.
    // - otherwise, the last tau[] will be > R, and will 
    // cause exit of the loop
    for (k = 0; s1 > taus[k]; k++)
      ;
    bw = bws[k];
    // make a new block numbered bn
    xlo = (x1/bw)*bw;  // round down
    xhi = xlo + bw;
    if (xhi > n) xhi = n;
    ylo = (y1/bw)*bw;
    yhi = ylo + bw;
    if (yhi > m) yhi = m;
    for (x2 = xlo; x2 < xhi; x2++)
      for (y2 = ylo; y2 < yhi; y2++)
        if (blocks[m*x2+y2] == Block_Clear)
          blocks[m*x2+y2] = bn;
    bn++;
      } else if (BLOCK_GROUPED(blocks[m*x1+y1])) {
    continue; // already in a block: this is OK
      } else {
    // should not happen
    mexWarnMsgTxt("create: bad blocknum.\n");
    return Block_Min; // anything < 0
      }
    }
  return (blocknum_t) ((int) bn - 1); // the last valid block number
}

/*
 * block_stats: compute statistics over all blocks
 * Valid blocks run from Block_First (=1) to btot, inclusive of both endpoints.
 * Returns 0 for out-of-range block number.
 */
static
int
block_stats(blockstat_t *stats, 
        blocknum_t *blocks, 
        double *z,             // z-coord
        double *a,             // the image to be convolved
        int m, 
        int n, 
        blocknum_t btot)
{
  int x1, y1;
  blocknum_t bn = Block_First;
  double wt, norm;

  // stats[].* is already zeroed out
  /* loop over image pixels, accumulating statistics */
  for (x1 = 0; x1 < n; x1++)
    for (y1 = 0; y1 < m; y1++) {
      bn = blocks[m*x1+y1];
      if (BLOCK_GROUPED(bn)) {
    if (bn > btot)
      return 0; // error: out-of-range block number
    wt = fabs(a[m*x1+y1]);
    stats[bn].x  += x1 * wt; // un-centered
    stats[bn].y  += y1 * wt; // un-centered
    stats[bn].z  += z[m*x1+y1] * wt;
    stats[bn].wt += wt;
    stats[bn].f  += a[m*x1+y1];
      }
    } // end for x1, y1
  /* loop over blocks, computing average */
  for (bn = Block_First; bn <= btot; bn++) {
    if (stats[bn].wt == 0)
      return 0; // error: empty block
    norm = 1.0 / stats[bn].wt;
    stats[bn].x *= norm;
    stats[bn].y *= norm;
    stats[bn].z *= norm;
  }
  return 1;
}

/*
 * block_check: check that blocks have been cleared
 */
static
int
block_check(blockstat_t *stats, int btot)
{
  int bn;
  int bn_nc = 0; // count blocks not cleared

  /* loop over blocks, checking for uncleared block */
  for (bn = Block_First; bn <= btot; bn++)
    if (stats[bn].f != 0)
      bn_nc++;
  return bn_nc;
}

/*
 * convert blocks for output to matlab:  int-to-double
 */
static
void
block_export(double *p, int *blocks, int m, int n)
{
  int inx;

  for (inx = 0; inx < m*n; inx++)
    p[inx] = (double) blocks[inx];
}


/******************************************************************
 *
 * Threading Support
 *
 * Structures to hold state passed down to computational threads,
 * a routine for thread-number selection
 *
 ******************************************************************
 *
 * Note on threading strategy:
 * The key problem with threading this, is that the convolution is
 * pretty wide (typically 200 pixels wide) and so the footprint of
 * pixels in the output `b' that are modified by looping over a 
 * rectangle of pixels in the input `a' is pretty large.
 *
 * We do the convolution by splitting the input into strips of size
 * Rwid along the x coordinate, each thread smoothing one strip at a 
 * time.  These partial results, which are wider than Rwid, are 
 * accumulated into the final output array.  E.g., consider:
 *
 *   Input: X = [x0|x1|x2|x3|x4|x5]
 *
 * The thread corresponding to x2 = X(:,Rwid*2+[1:Rwid]), say, computes 
 * a result y2 which is roughly of width Rwid + Kswath, where Kswath is
 * the full width of the kernel.  This partial result must be summed in
 * to the full output starting at the correct offset.  It will extend
 * Kswath/2 rows "below" Rwid*2, and Kswath/2 rows "above" Rwid*3.
 * (The partial result is located in a thread-specific buffer, but the 
 * full output array is shared across threads.)
 * While it's being accumulated, other threads must not alter the 
 * output array.
 *
 * We pay a penalty for separating the input into chunks, because 
 * the corresponding output chunks are enlarged.  Smaller chunks
 * means more granularity in the threading, so all threads are always 
 * working (good), but too small means that we spend too much time 
 * accumulating partial results into a full output (bad).  In general,
 * each output array location will be accumulated-into 
 * 
 *    1 + Kswath/Rwid times,
 *
 * because in total, if there are k strips, there will be k*Rwid output
 * rows, but k*(Rwid+Kswath) result-rows found.  Therefore each output
 * row will be written (Rwid+Kswath)/Rwid = 1 + Kswath/Rwid times.
 * For HMI, Kswath = 200
 *
 * We put most of the logic for partitioning into the individual threads. 
 * We maintain an index of the largest row still remaining unfinished, and
 * let each thread check that index to determine whether to carve off 
 * another rectangle.
 * 
 ******************************************************************/

/*
 * these structures are needed because the pthread interface passes
 * only one item to the thread as an argument, so we put all the
 * convolution arguments in that one argument
 */

// one instance of this flavor of structure holds common info
struct thread_info {
  // unchanging info, like geometry and overall parameters
  // (see smooth_rect for descriptions)
  // sizes
  int m; // m = y
  int n; // n = x
  // strides
  int strIx;
  int strIy;
  // disk geometry
  double xcen;
  double ycen;
  double R;
  double cos_a;
  double sin_a;
  // input matrices
  double *a;
  double *s;
  // output matrix
  double *b;
  // kernel
  double *kernel;
  int Klen;
  double Ktop;
  int Kswath;
  double *Kwt;
  // block structure of input
  blocknum_t *blocks;
  blockstat_t *stats;
  // swath rectangle sizes
  int xwid; // max size along x direction
  int ywid; // max size along y direction
  int Rwid; // target size for each thread's rectangle
  // next row to be assigned to a thread
  int smooth_bb_top;
};

// MAX_THREAD instances of this flavor of structure holds per-thread info
struct thread_args {
  // per-thread info
  int tnum;  // my own thread number, just informational
  // output error status (zero is OK, nonzero is not)
  int status;
};

/*
 * thread info, declared as local to this file for simplicity
 * (in particular, SA could be local to a function, and it could contain
 * pointers to the mutex and to SI, but this way keeps the declaration
 * close to the structure definition)
 */
// this is a hard max, but it can be lowered by an environment variable
#define MAX_THREAD 8
static struct thread_info SI; // one slot for common info
static struct thread_args SA[MAX_THREAD]; // slots for per-thread info
static pthread_mutex_t blockstat_mutex;
static pthread_mutex_t accum_mutex;
static pthread_mutex_t smooth_bb_top_mutex;

/*
 * get the requested number of threads, either from environment
 * (MXT_NUM_THREAD) or from the compiled-in default
 */
static
int
get_thread_num()
{
  int tn, tn1;
  char *thread_env;

  if ((thread_env = getenv("MXT_NUM_THREAD")) != NULL) {
    // get from environment
    tn = strtol(thread_env, NULL, 10);
    if (tn == 0)
      tn = MAX_THREAD; // error in conversion
  } else {
    // use compiled-in default
    tn = MAX_THREAD;
  }
  // clip to 1..thread_max
  if (tn > MAX_THREAD)
    tn = MAX_THREAD;
  else if (tn <= 0)
    tn = 1;
  return tn;
}

/********************************************************************************
 *
 * Spherical convolution utility
 *
 ********************************************************************************/

/* swath_extent: establish loop boundaries
 * 
 * move Kswath in each direction away from (x0,y0), 
 * adjusting for both rotation angle (cos_a, sin_a) and
 * distance weighting (Kwt), to determine the four corners
 * of a box containing the region to include in the 
 * convolution.  Returns the largest and smallest x and y
 * indices which contain the whole box.
 */
static
void
swath_extent(int *xlo_out, int *xhi_out, 
         int *ylo_out, int *yhi_out, 
         double x0, double y0,
         double xcen, double ycen, double R, 
         double cos_a, double sin_a, 
         double *Kwt, int Kswath, int m, int n)
{
  int xB, yB;
  double x1, y1;
  double xp1, yp1;
  double xlo, ylo, xhi, yhi;

  // high and low markers
  xlo = ylo = m+n; // higher than either alone
  xhi = yhi = 0;   // lower than either
  // loop over 4 corners
  for (xB = -1; xB <= 1; xB += 2)
    for (yB = -1; yB <= 1; yB += 2) {
      // scaled point
      x1  = xB * Kswath / sqrt(Kwt[0]);
      y1  = yB * Kswath / sqrt(Kwt[1]);
      // rotate by p-angle, place in image coords
      xp1 = x0 +  x1 * cos_a + y1 * sin_a;
      yp1 = y0 + -x1 * sin_a + y1 * cos_a;
      // save low and high
      if (xp1 < xlo) xlo = xp1;
      if (xp1 > xhi) xhi = xp1;
      if (yp1 < ylo) ylo = yp1;
      if (yp1 > yhi) yhi = yp1;
    }
  // clip to Rsun
  if (xlo < xcen-R) xlo = xcen-R;
  if (xhi > xcen+R) xhi = xcen+R;
  if (ylo < ycen-R) ylo = ycen-R;
  if (yhi > ycen+R) yhi = ycen+R;
  // round outward to integers
  xlo = floor(xlo);  ylo = floor(ylo);
  xhi = ceil (xhi);  yhi = ceil (yhi);
  // clip to image
  if (xlo < 0)   xlo = 0;
  if (xhi > n-1) xhi = n-1;
  if (ylo < 0)   ylo = 0;
  if (yhi > m-1) yhi = m-1;
  // make output
  *xlo_out = (int) xlo;
  *xhi_out = (int) xhi;
  *ylo_out = (int) ylo;
  *yhi_out = (int) yhi;
}

/* 
 * 
 * swath_max_size: determine an upper size limit on the footprint
 * of the kernel.  In other words, this returns an upper bound on the
 * length of the x and y dimensions, generated by swath_extent above,
 * needed to hold a convolution result array.  We want to be able to 
 * alloc xwid*ywid doubles to hold the result from convolving one input pixel.
 *
 * This routine assumes that the size of the swath can be upper-bounded 
 * independent of position (x0,y0).  Currently, the only need for an 
 * upper bound (rather than an exact figure) is that differences in 
 * rounding can cause +/-1 changes in size at each end of the convolution 
 * footprint.  Gross differences are not possible.  
 */
static
void
swath_max_size(int *xwid, 
           int *ywid,
           double xcen, double ycen, double R, 
           double cos_a, double sin_a, 
           double *Kwt, int Kswath, int m, int n)
{
  int xlo, xhi, ylo, yhi;

  swath_extent(&xlo, &xhi, &ylo, &yhi, 
           xcen, ycen,  // use image center as the center of the swath
           xcen, ycen, R, 
           cos_a, sin_a, 
           Kwt, Kswath, m, n);
  // first, +1 because a range of lo = hi still includes 1 pixel
  // then,  +2 because an unpredictable roundoff at either the lower
  // or upper endpoint, or both, might result in a larger range.
  *xwid = (xhi - xlo + 1) + 2;
  *ywid = (yhi - ylo + 1) + 2;
}

/********************************************************************************
 *
 * Spherical convolution core
 *
 ********************************************************************************/

/*
 * Background on the convolution
 *
 * There are two aspects of this code that are different from standard
 * convolution (say, in 1 or 2 dimensions, in a periodic setting):
 *  a - the kernel can be viewed as spatially varying
 *  b - the discretization is not uniform
 * 
 * Point (a) is debatable.  If we imagine the functions x, y, k defined
 * on the sphere, then k can be expressed as depending only on s-t where
 * we are convolving x(t) and k(s).  That is, we can write
 *   y(s) = int_S x(t) k(s-t) dt
 * where t, s are 3-vectors of unit norm (on the sphere).  In this case,
 *   ||s - t||^2 = (s-t)'(s-t) = 2(1-t's), or
 *   t's = 1 - 0.5 ||s-t||^2
 * and the dot product is a function of only s-t.
 * In this interpretation, the kernel is not spatially varying.  The 
 * convolution is computed on the sphere S and then, in a separate 
 * operation, projected down to the unit disk.
 * 
 * As an aside, if an exponential kernel k is used as input this function,
 *   k(dot) = const * exp(lambda * dot), where dot = t's,
 * then this is equivalent to a gaussian function of the difference
 *   delta = s - t
 * since
 *   k(dot) = k(t's) = const * exp(lambda * 0.5 * ||s - t||^2)
 * This is a pretty good rationale for using a kernel that is of exponential
 * form, i.e. k = const1 * exp(const2 * linspace(-1, 0, kparam(2))), where
 * const2 sets the fwhm and const1 is chosen so k is of unit norm.
 * 
 * Point (b) is tricky.  It arises because one pixel at disk-center 
 * corresponds to seeing a region of (say) "one unit", but at the limb,
 * the pixel could correspond to seeing many units.  The conversion factor
 * is 1/z, where z is the distance from the x-y (image) plane outward to
 * the sphere boundary (z = sqrt(1-x^2-y^2)).  The limb pixels carry
 * information about more area, so they will influence other pixels 
 * proportionally more strongly.  We just need to quantify the effect.
 *
 * Before we go on to do that, a word about the (1/z) factor, which is
 * the area of a pixel.  As we go closer to the limb, z goes (slowly)
 * to zero, and 1/z, which approximates the area, goes (slowly) to 
 * infinity.  Especially for high-pixel-density images, this will 
 * make the integral blow up.  The area ~= 1/z approximation thus
 * breaks down there.  One way out is to say that a pixel with a center
 * precisely at the limb (z=0) may as well be regarded as having area
 * equal to that of its maximally-on-disk part, which is half a pixel 
 * farther onto the disk.  Drawing the right triangle, we see that the 
 * z of a point 0.5 pixel from the limb is sqrt(R) (where the max z, at
 * disc center, is R).  To tame the normalization factor, we take 
 * sqrt(R) as the lowest allowed z for the area computation.  This
 * typically does not arise in practice -- even for HMI, at 3 pixels in,
 * z = sqrt(6*R) ~ 100, the min-z = sqrt(R) ~ 44.
 * 
 * It's natural to compute the convolution as a discrete approximation 
 * to a continuous convolution:
 * 
 *   I-1:  y(s) = \int_S  x(t) k(s,t) dt         (using k(s,t) for k(s-t))
 *         y(i) = \sum_j  x(j) k(i,j) |S_j|      [*]
 * 
 * where the whole sphere S is partitioned into {S_j}, and we assume that 
 * to first order, the function x(t) equals a number x(j) on the patch S_j, 
 * which has volume |S_j|.
 * The problem is that one can expand the integral the other way:
 * 
 *   I-2:  y(s) = \int_S  x(s-t) k(t,0) dt       (using k(t,0) for k(t-0))
 *         y(i) = \sum_j  x(i-j) k(j,0) |S_j|
 *
 * and the discretized version weights x(i-j) differently.  There are some
 * subtleties here, because I-2 assumes the kernel is translation-invariant.  
 * But, in any case, the two forms are not symmetric.  And one then 
 * wonders about the exact interpretation of the terms within [*].  
 * (E.g., does k(i,j) have to be k(s,t) weighted by |S_i| or |S_j|?)
 * 
 * These questions become focused when you wonder how to normalize k() so
 * that the convolution integral has unit norm.
 *
 * It's possible to clarify these questions, which have to do with 
 * weighting of nonnegative functions x and y, by analogy with 
 * probability. The analogy is that x(t) and y(s) are both 
 * densities with respect to a reference measure A(.), and k(s,t) 
 * is a conditional density tying them together:
 * 
 *    x(t): pdf of an r.v. X relative to measure A(t)
 *    y(s): pdf of an r.v. Y   "      "      "    "
 *    k(s,t): conditional density p(Y=s|X=t) "    "
 * 
 * If the density x(t) is relatively smooth over patches S_j, we can 
 * discretize everything, where |S_j| = A(S_j), the area of the patch:
 * 
 *    x_j: pmf of the r.v. X.  Thus Pr(X in S_j) = x_j = x(t) |S_j| 
 *         for any t in S_j
 *    y_i: pmf of the r.v. Y.  Thus Pr(Y in S_i) = y_i = y(s) |S_i| 
 *         for any s in S_i
 *    k_(i,j): conditional density p(Y|X), so k_(i,j) = k(s,t) |S_i|
 *            (t and s as above).  This is so because we are *given* 
 *            the value of X, so the average is over Y (S_i).
 *    
 * To relate this back to the original problem, at large (limb) patches,
 * this says we are equally sure, throughout the patch, that the
 * patch has a given label (x(t) = const on the patch).  But, over
 * the whole patch, there is a greater mass than the at disk center
 * (since the |S_j| term is larger).  So, the limb patch as a whole
 * will have greater influence than disk-center patches as a whole.
 *
 * In short: the x(t), y(s) are the observed and smoothed values.
 * The x_j, y_i capture the discretization effects that make some
 * pixels "more strong" than others.
 * 
 * To relate Y to X, note
 * 
 *   Pr(Y in S_i) = \sum_j  Pr(X in S_j) * Pr(Y in S_i | X in S_j)
 *                = \sum_j  [x(t) |S_j|] * [k(s,t) |S_i|]
 * 
 * and, noting Pr(Y in S_i) = y_i = y(s) |S_i|,
 *
 *   y(s) = \sum_j  x(t) k(s,t) |S_j|      [t from S_j & s from S_i always]
 * 
 * This recovers the formula [*]. 
 *
 * It also tells us how to interpret "unit norm".  For y(s) and  x(t),
 * it means integrating to one against the area measure A(.).  This is
 * identical to the discretized y_i and x_j summing to one against
 * |S_i|.  In short: the discretized versions sum to unity, not the
 * continuous ones.  In matlab, this would mean:
 *   >> z = sqrt(1 - x^2 - y^2)  or
 *   >> z = rsun - ringmatrix([m  n], [], [], center)
 *   >> z(z < 0) = Inf;  % make them drop out
 * and then, since the patch areas |S_j| are proportional to 1/z:
 *   >> sum(sum(x .* (1./z))) = 1
 *   >> sum(sum(y .* (1./z))) = 1
 * 
 * More important, for k(s,t), it means 
 * 
 *   (forall t)  \int_S  k(s,t) dA(s) = 1
 *
 * Because k(s,t) depends only on the radial distance between 
 * s and t, we can take t = (0 0 1), and then s runs from 0..1 with
 * dot = s't = sqrt(1 - s^2).  However, this simplifies considerably.
 * As we have parameterized it, k(s,t) = k(dot), which is given 
 * explicitly in the input k(.).  We need only that
 *
 *   \int_{dot=0..1}  k(dot) dA(dot) = 1
 *
 * In fact, the minimum value of dot is kparam(1).  The area element
 * corresponding to s't in [dot, dot+d_dot] is just 2 pi d_dot.  Thus
 * we want
 *
 *   2 pi \int_{dot=kparam(1)..1}  k(dot) d_dot = 1
 *
 * The integral, when discretized, needs only the width of a length
 * element, which is delta_dot (roughly corresponding to d_dot above):
 * 
 *   >> delta_dot = (1 - kparam(1))/length(k);
 *   >> normer = 2*pi*delta_dot*sum(k);
 *
 * Letting k = k/normer makes k have unit norm.
 */

// error codes for smooth(), must be < 0
#define SMOOTH_NOMEM    -1
#define SMOOTH_INTERNAL -2
#define SMOOTH_NOTHREAD -3

/*
 * smooth_rect: convolve input a with kernel and put into output b
 *
 * Does convolution only in an input rectangle given by {x,y}R{lo,hi}.
 * The output will be written into a larger range, depending on the kernel 
 * width Kswath, the kernel weights Kwt, and the tilt angle (cos_a, sin_a).
 * There are separate parameters for input and output strides, so the two
 * can be differently-sized.
 * This routine operates in a thread, but is not highly aware that it's 
 * threaded (we need one mutex reference).
 * 
 * returns 1 if error, 0 otherwise
 * (the only error condition is an internal logic error)
 */
static
int
smooth_rect(int m,          // input image size (m = y)
        int n,          //                  (n = x)
        // rectangle and strides
        int xRlo,       // x rect range lo/hi (input images)
        int xRhi,
        int yRlo,       // y rect range lo/hi (input images)
        int yRhi,
        int xROlo,      // x rect offset (output image offset from input, typ. >= 0)
        int yROlo,      // y rect offset (output image offset from input, typ. >= 0)
        int strIx,      // stride in x, input (a, s arrays)
        int strIy,      // stride in y, input
        int strOx,      // stride in x, output (b array)
        int strOy,      // stride in y, output
        // disk geometry
        double xcen,    // center, n or x, in C coords, origin at 0
        double ycen,    // center, m or y, in C coords, origin at 0
        double R,       // radius, in pixels
        double cos_a,   // cos(alpha)
        double sin_a,   // sin(alpha)
        // input + output matrices
        double *a,      // input is mXn
        double *s,      // mXn (z coordinate, already found)
        double *b,      // output (convolved result)
        // kernel
        double *kernel, // kernel lookup table
        int Klen,       // length of the table
        double Ktop,    // upper limit of the table (0..Ktop)
        int Kswath,     // box enclosing kernel
        double *Kwt,    // weights on x, y, z
        // block structure of input
        blocknum_t *blocks,
        blockstat_t *stats)
{
  const double Kstep = (Klen-1)/Ktop;  // map [0,Ktop] -> [0,K-1]
  const double min_Zwt = sqrt(R); // sanity bound on area-weighting, see above
  int x1, y1;             // p1 coords in 2D, bounded by xRlo, xRhi, etc.
  int x2, y2;             // p2 coords in 2D, bounded by xlo, xhi, etc.
  int xlo, xhi, ylo, yhi; // inner loop boundaries
  double p1x, p1y, p1z;   // point P1 coords in 3D
  double p2x, p2y, p2z;   // point P2 coords in 3D
  int p1_off;             // array offset into a,s for the point (x1,y1)
  int x2_off, y2_off;     // array offsets into a,s for the point (x2,y2)
  int xO_off, yO_off;     // array offsets into b for the point (x2,y2)
  blocknum_t bn;          // block number of (x1,y1)
  double dist;            // distance from p1 to p2
  double contribution;    // point (x1,y1)'s contribution to integral
  double scale;           // scale factor for convolution integral
  double xx_term;
  double wt_xx, wt_yy, wt_zz, wt_xy;
  double inxC;            // real-valued index into kernel table
  int inxD;               // integer index into kernel table

  /* set up weights */
  wt_xx = Kwt[0] * cos_a * cos_a + Kwt[1] * sin_a * sin_a;
  wt_yy = Kwt[0] * sin_a * sin_a + Kwt[1] * cos_a * cos_a;
  wt_xy = 2.0 * (Kwt[0] - Kwt[1]) * cos_a * sin_a;
  wt_zz = Kwt[2];
  /* include 1/(R^2) factor in each */
  wt_xx /= R*R;
  wt_yy /= R*R;
  wt_zz /= R*R;
  wt_xy /= R*R;

  /* loop over image pixels */
  for (x1 = xRlo; x1 <= xRhi; x1++) 
    for (y1 = yRlo; y1 <= yRhi; y1++) {
      p1_off = x1*strIx + y1*strIy;
      bn = blocks[p1_off]; // block number
      if (bn == Block_Clear) 
    return 1; // should never happen
      if (bn == Block_Offdisk || bn == Block_Skip)
    continue; // off-limb or zero weight: no action needed
      if (bn == Block_Singleton) {
    // singleton
    /* find original vector coordinates P1 = (px1, py1, pz1) */
    p1x = x1;
    p1y = y1;
    p1z = s[p1_off];
    contribution = a[p1_off];
      } else {
    // look up coordinates and function value via per-block statistics
    p1x = stats[bn].x;  // un-centered
    p1y = stats[bn].y;  // un-centered
    p1z = stats[bn].z;
    // skip the mutex in the case where the block was already counted
    if (stats[bn].f == 0) continue; 
    // ensure only this thread is convolving with this block
    pthread_mutex_lock(&blockstat_mutex);
    contribution = stats[bn].f;
    stats[bn].f = 0; // indicate this thread has put it in the integral
    pthread_mutex_unlock(&blockstat_mutex);
      }
      // contrib == 0 when the block has been included once already,
      // in this thread or in another thread
      if (contribution == 0)
    continue; 
      /* establish loop boundaries */
      swath_extent(&xlo, &xhi, &ylo, &yhi, p1x, p1y,
           xcen, ycen, R, cos_a, sin_a, Kwt, Kswath, m, n);
      if (0)
    printf("(%.1f,%.1f) loop [%d,%d]x[%d,%d] (%dx%d)\n", 
           p1x, p1y, xlo, xhi, ylo, yhi, xhi-xlo+1, yhi-ylo+1); 
      // scale factor is the target of the convolution divided by an
      // area-based weight; the lowest weight allowed is min_Zwt
      scale = contribution / (((p1z < min_Zwt) ? min_Zwt : p1z) * R);
      /* loop over kernel */
      for (x2 = xlo, 
         x2_off =  xlo         *strIx, 
         xO_off = (xlo - xROlo)*strOx; 
       x2 <= xhi; 
       x2++, 
         x2_off += strIx,
         xO_off += strOx) {
    // optimizations for xx_term, offsets help a little (~5%): Apr 2010
    xx_term = (p1x - x2)*(p1x - x2)*wt_xx;
    for (y2 = ylo, 
           y2_off =  ylo         *strIy,
           yO_off = (ylo - yROlo)*strOy; 
         y2 <= yhi; 
         y2++, 
           y2_off += strIy,
           yO_off += strOy) {
      if (blocks[x2_off+y2_off] == Block_Offdisk) continue;
      p2z = s[x2_off+y2_off];
      dist = xx_term +
             (p1y -  y2) * (p1y -  y2) * wt_yy + 
             (p1y -  y2) * (p1x -  x2) * wt_xy + 
             (p1z - p2z) * (p1z - p2z) * wt_zz;  // distance ||P1-P2||
      /* do kernel lookup on distance */
      if (dist >= Ktop) continue; // dist >= Ktop: no contribution to b
      inxC = dist * Kstep;
      inxD = floor(inxC);
      if (inxD >= Klen-1) continue; // rounding puts dist on cusp: skip 
      // note, now 0 <= inxD < Klen-1
      dist = inxC - inxD;  /* in [0,1), smaller means closer to inxD */
      /*
      printf(" adding %g->%d, wt %g, interp %g\n", inxC, inxD, kernel[inxD], dist);
      */
      b[xO_off+yO_off] += 
        (kernel[inxD]*(1.0-dist) + kernel[inxD+1]*dist) * scale;
    } // for y2
      } // for x2
    } // for (x1,y1)
  return 0;
}

/*
 * smooth_thread: shell for Pthread invocation of the smoothing operation
 *
 * Each thread operates largely independently, selecting a rectangular region
 * of the input a, dividing along x and taking all of y, and sending that to
 * the core smoother.  Once the smoothed rectangle comes back, it is accumulated 
 * into the overall output matrix.
 * The thread control at this level needs two mutexes.  One to control the
 * division of the input among the threads (smooth_bb_top, which indicates
 * the next un-smoothed row in the input).  Another to control accumulation
 * into the result.  We make it simple, so that only one thread can touch the
 * output at a time.  By looking at timings, there seems to be little contention
 * at this point.
 *
 */
static
void *
smooth_thread(void *arg)
{
  const int verbose = 0;
  const int bb_max = SI.n-1; // the maximum value the rectangle needs to encompass
  double * const b = SI.b;   // points to output matrix
  double *b1;                // individual output matrix, lumped into b
  int rc;
  struct thread_args *sa;
  struct timeval tv1, tv2;
  int bb_top, bb_low, bb_hi; // bounding box boundary bookkeeping
  int bb_total, bb_count;    // internal counts for informational printouts
  int b_zero;                // offset to y origin of bb passed to smooth_rect
  int xsize, ysize;
  int x, y;
  int p1_off, p2_off;

  // arguments are taken from global arg list SI, and per-thread arg list *sa
  gettimeofday(&tv1, NULL);
  sa = (struct thread_args *) arg;
  // split the input along x (0...n-1)
  bb_total = bb_count = 0;
  while (1) {
    // carve off a new range of rows
    pthread_mutex_lock(&smooth_bb_top_mutex);
    bb_top = SI.smooth_bb_top; // NB, bb_top row itself is NOT assigned to a thread
    if (bb_top <= bb_max) {
      // more rows to do: find new range
      bb_low = bb_top;
      bb_hi  = bb_top + SI.Rwid - 1; // bb_hi is assigned to this thread
      if (bb_hi > bb_max) bb_hi = bb_max; // cap at the max x
      // install new top
      SI.smooth_bb_top = bb_hi + 1;
    }
    pthread_mutex_unlock(&smooth_bb_top_mutex);
    if (bb_top > bb_max)
      break; // there were no more rows to do
    bb_count++;
    bb_total += bb_hi - bb_low + 1;
    // new range of rows is [bb_low,bb_hi], inclusive
    // calloc a rectangle big enough to hold the convolved rectangle
    xsize = (bb_hi - bb_low + 1) + SI.xwid;
    b_zero = bb_low - SI.xwid/2; // b1(y==b_zero) will go into b(y==bb_low)
    ysize = SI.m;
    // fresh calloc each iteration is faster than clearing b1 with memset()
    b1 = calloc(xsize*ysize, sizeof(*b1)); 
    rc = smooth_rect(SI.m, SI.n, 
             // bounding rectangle (input images)
             bb_low, bb_hi,
             0, SI.m-1,
             b_zero, 0,
             // strides (input images)
             SI.strIx,  SI.strIy, 
             // strides (output image)
             ysize, 1,
             // disk geometry
             SI.xcen,  SI.ycen, SI.R, SI.cos_a, SI.sin_a, 
             // input matrices
             SI.a, SI.s, 
             // output matrix
             b1,
             // kernel
             SI.kernel, SI.Klen, SI.Ktop, SI.Kswath, SI.Kwt, 
             // block structure of input
             SI.blocks, SI.stats);
    sa->status = rc; // exit status, 0 is OK
    if (rc != 0) {
      free(b1);
      pthread_exit(NULL);  // trouble
    }
    //printf("    thread #%d: handled [%d,%d] * %d @ %d\n", sa->tnum, 
    //   bb_low, bb_hi, SI.xwid, b_zero);
    // accumulate into the final result
    pthread_mutex_lock(&accum_mutex);
    // x loop is over all of b1, some of which may be aligned with
    // an off-image part of the output b
    for (x = 0; x < xsize; x++) 
      for (y = 0; y < ysize; y++) {
    // off-image check: is (x,y) on-image for b?
    if ((x + b_zero < 0) || (x + b_zero >= SI.n)) continue;
    p1_off = x*ysize + y*1; // offset into b1
    p2_off = (x+b_zero)*SI.strIx + y*SI.strIy; // offset into b
    b[p2_off] += b1[p1_off];
      }
    pthread_mutex_unlock(&accum_mutex);
    free(b1);
  } 
  // thread is done
  gettimeofday(&tv2, NULL);
  int dt = (((long) tv2.tv_sec) - ((long) tv1.tv_sec)) * 1000000 + 
    (((int) tv2.tv_usec) - ((int) tv1.tv_usec));
  if (verbose)
    mexPrintf("    ST: exiting thread #%d, did %d/%d, dt = %.1f ms\n", 
          sa->tnum, bb_count, bb_total, dt/1000.0);
  pthread_exit(NULL);
}


/* control:  coordinates the smoothing operation
 *
 * sets up blocking data structures,
 * launches threads, combines results
 */
static
int
control(int m,          // image size
    int n,
    double *a,      // input is mXn
    double *b,      // output is mXn (convolved result)
    double *s,      // output is mXn (z coordinate)
    double *p,      // output, mXn (patch indexes)
    double xcen,    // center, n or x, in C coords, origin at 0
    double ycen,    // center, m or y, in C coords, origin at 0
    double R,       // radius, in pixels
    double p0,      // tilt angle
    double *kernel, // kernel lookup table
    int Klen,       // length of the table
    double Ktop,    // upper limit of the table (0..Ktop)
    int Kswath,     // box enclosing kernel
    double *Kwt,    // weights on x, y, z; angle
    double *bws_in, // block lookup table (length x 2)
    int bwlen)      // its length

{
  const int verbose = 0;
  const int mn = m*n;
  const double cos_a = cos(p0); // cos(alpha)
  const double sin_a = sin(p0); // sin(alpha)
  const double nan = mxt_getnand(); // cache nan
  int s_alloc_local;            // did we alloc s here?
  blocknum_t blockNum, bn;
  blocknum_t *blocks;
  blockstat_t *stats;
  double tau1, *taus;           // z-thresholds for blocking
  int *bws;                     // block widths (int's)
  int x1, y1, off;              // array position and offset
  double p1x, p1y;              // image position
  double temp;
  int xwid, ywid, Rwid;         // prototype rectangle widths
  int t, NT;                    // thread counts
  pthread_t threads[MAX_THREAD];// thread info
  struct timeval tv1, tv2;
  int dt;

  gettimeofday(&tv1, NULL);
  /*
   * 1: Allocations
   */
  // s, if necessary
  if (s_alloc_local = (s == NULL))
    s = calloc(mn, sizeof(double));
  if (!s) return SMOOTH_NOMEM;
  // block-number array
  blocks = calloc(mn, sizeof(blocknum_t)); /* must be cleared */
  if (!blocks) return SMOOTH_NOMEM;

  /*
   * 2: Find z
   */
  // loop over image pixels
  off = 0;  // in general, off = x1*m + y1
  for (x1 = 0; x1 < n; x1++)
    for (y1 = 0; y1 < m; y1++) {
      // find original vector coordinates P1 = (px1, py1, pz1)
      p1x = x1 - xcen;
      p1y = y1 - ycen;
      temp = (R + p1x) * (R - p1x) - p1y*p1y; /* R*R-x*x = (R-x)*(R+x) */
      if (temp < 0)
    b[off] = s[off] = nan;  // off-disk
      else
    s[off] = sqrt(temp);        // +sqrt: visible half-sphere
      off++;
    }
  gettimeofday(&tv2, NULL);
  dt = (((long) tv2.tv_sec) - ((long) tv1.tv_sec)) * 1000000 + 
    (((int) tv2.tv_usec) - ((int) tv1.tv_usec));
  if (verbose)
    mexPrintf("  do_SSa: dt = %.1f ms\n", dt/1000.0);

  /*
   * 3: Set up blocking info
   */
  // juggle the block table we were given
  taus = calloc(bwlen, sizeof(*taus));
  bws  = calloc(bwlen, sizeof(*bws));
  if (!taus || !bws) return SMOOTH_NOMEM;
  block_lut_reformat(&tau1, taus, bws, bws_in, bwlen, R);
  // sets up blocks, finds number of blocks
  blockNum = block_create(blocks, s, a, tau1, taus, bws, m, n);
  if (blockNum == Block_Min)
    return SMOOTH_INTERNAL; // logic error in code
  // don't need these any more
  free(taus); free(bws); taus = NULL; bws = NULL;
  // printf("Using %d blocks\n", blockNum);
  // compute aggregate block statistics
  stats = calloc(blockNum+1, sizeof(*stats));
  if (!stats) return SMOOTH_NOMEM;
  if (!block_stats(stats, blocks, s, a, m, n, blockNum))
    return SMOOTH_INTERNAL;  // out-of-range block number
  gettimeofday(&tv2, NULL);
  dt = (((long) tv2.tv_sec) - ((long) tv1.tv_sec)) * 1000000 + 
    (((int) tv2.tv_usec) - ((int) tv1.tv_usec));
  if (verbose)
    mexPrintf("  do_SSb: dt = %.1f ms\n", dt/1000.0);

  /*
   * 4: Set up thread info
   */
  NT = get_thread_num(); // check env for advice
  // determine maximum rectangle width 
  swath_max_size(&xwid, &ywid, xcen, ycen, R, cos_a, sin_a, Kwt, Kswath, m, n);
  // could make Rwid dependent on Kswath...
  // n=4096, NT=8, Kswath = 200 => Rwid = 128, 1+Ks/Rwid = 2.5
  Rwid = n / (NT * 4);
  if (Rwid < 4) Rwid = 4; // best not to go too small, certainly need > 0
  if (NT == 1) Rwid = n; // if one thread, do it all at once
  // record rectangle extent, common across threads
  SI.xwid = xwid; // max size along x direction
  SI.ywid = ywid; // max size along y
  SI.Rwid = Rwid; // target size for each thread's rectangle
  // set up geometry and param info valid for all calls to the smoother
  SI.m = m; SI.n = n;
  SI.strIx = m; SI.strIy = 1; // could generalize these strides
  SI.xcen = xcen; SI.ycen = ycen;
  SI.R = R; SI.cos_a = cos_a, SI.sin_a = sin_a;
  // input arrays (image size)
  SI.a = a; SI.s = s;
  // output array (image size)
  SI.b = b;
  // kernel
  SI.kernel = kernel; SI.Klen = Klen; SI.Ktop = Ktop;
  SI.Kswath = Kswath; SI.Kwt = Kwt;
  SI.blocks = blocks; SI.stats = stats;

  /*
   * 5: Do the smoothing in threads
   */
  pthread_mutex_init(&blockstat_mutex, NULL);
  pthread_mutex_init(&accum_mutex, NULL);
  pthread_mutex_init(&smooth_bb_top_mutex, NULL);
  // shared across threads, this marks our current position in filling 
  // in the output -- it indicates the "next row to convolve"
  SI.smooth_bb_top = 0; 
  // initiate NT threads
  for (t = 0; t < NT; t++) {
    // pack changing args
    SA[t].tnum = t;  // informational for the thread
    SA[t].status = 0; // status output
    // start the thread
    if (pthread_create(threads+t, NULL, smooth_thread, (void *) &SA[t]) != 0)
      return SMOOTH_NOTHREAD;
  }
  // join all threads
  for (t = 0; t < NT; t++) {
    if (pthread_join(threads[t], NULL) != 0)
      return SMOOTH_NOTHREAD;
    if (SA[t].status != 0)
      return SMOOTH_INTERNAL; // internal error in the rectangle smoother
  }
  pthread_mutex_destroy(&blockstat_mutex);
  pthread_mutex_destroy(&accum_mutex);
  pthread_mutex_destroy(&smooth_bb_top_mutex);
  gettimeofday(&tv2, NULL);
  dt = (((long) tv2.tv_sec) - ((long) tv1.tv_sec)) * 1000000 + 
    (((int) tv2.tv_usec) - ((int) tv1.tv_usec));
  if (verbose)
    mexPrintf("  do_SSc: dt = %.1f ms\n", dt/1000.0);

  /*
   * 6: Clean up and exit
   */
  if (block_check(stats, blockNum) != 0) {
    mexWarnMsgTxt("Error: did not integrate with all blocks\n");
    return SMOOTH_INTERNAL;
  }
  // optionally copy blocks into p
  if (p)
    block_export(p, blocks, m, n); 
  free(blocks);
  free(stats);
  if (s_alloc_local)
    free(s);
  return blockNum; // always >= 0
}

#ifdef StaticP  /* undefined under mex */
StaticP
#endif
void
mexFunction(int nlhs, 
        mxArray *plhs[], 
        int nrhs,
        const mxArray *prhs[])
{
  const int verbose = 0;
  char errstr[160];
  char *msg;
  double datamin, datamax;
  double *s_data;             // place to put s in case we don't return it 
  double *p_data;             // place to put p in case we don't return it 
  mxArray *wt_holder;         // place to put KWt if it wasn't supplied 
  int m, n;                   // image size 
  double xcen, ycen, rsun, p0;// disk params 
  double *Klut;               // kernel lookup table
  int Nker;                   // and its length
  double Ktop;                // lut spans [0,Ktop]
  double KswathR; int KswathI;// box enclosing kernel 
  double *bws; int bwlen;     // block width lookup table
  int ok;
  struct timeval tv1, tv2;

  gettimeofday(&tv1, NULL);
  /* Hook for introspection (function signature, docstring) */
  if (nrhs < 0) { 
    plhs[0] = mxt_PackSignature((mxt_Signature) (-nrhs), 
                NARGIN_MIN, NARGIN_MAX, 
                NARGOUT_MIN, NARGOUT_MAX, 
                in_names, in_specs, out_names, docstring);
    return;
  }
  /* check number of arguments */
  if ((nrhs < NARGIN_MIN) || (nrhs > NARGIN_MAX))
     mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                "%s: Expect %d <= input args <= %d",
                progname, NARGIN_MIN, NARGIN_MAX), errstr));
  if ((nlhs < NARGOUT_MIN) || (nlhs > NARGOUT_MAX))
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                 "%s: Expect %d <= output args <= %d",
                 progname, NARGOUT_MIN, NARGOUT_MAX), errstr));
  mexargparse(nrhs, prhs, in_names, in_specs, NULL, progname);

  /* get size of input image x */
  m = (int) mxGetM(prhs[ARG_X]);  // m = y 
  n = (int) mxGetN(prhs[ARG_X]);  // n = x 
  /* get disk params */
  xcen = mxGetPr(prhs[ARG_GEOM])[0] - 1; // arg #1 = x0 (convert to C origin)
  ycen = mxGetPr(prhs[ARG_GEOM])[1] - 1; // arg #2 = y0 (convert to C origin)
  rsun = mxGetPr(prhs[ARG_GEOM])[2];     // arg #3 = radius
                                         // (do not use beta angle now)
  p0   = mxGetPr(prhs[ARG_GEOM])[4];     // p-angle
  p0  *= M_PI/180.0;                     // to radians
  // kernel params
  Ktop   = mxGetPr(prhs[ARG_KParam])[0];
  KswathR= mxGetPr(prhs[ARG_KParam])[1];
  KswathI= (int) KswathR;
  if (Ktop <= 0 || Ktop >= 1)
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                 "%s: Expect %g <= Kpar(1) < %g, got %g",
                 progname, 0.0, 1.0, Ktop), errstr));
  if (KswathI < 1 || KswathI != KswathR)
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                 "%s: Expect integer Kpar(2) >= %d, got %g",
                 progname, -1, KswathR), errstr));
  // kernel LUT - make or use supplied one
  if (mxGetNumberOfElements(prhs[ARG_K]) == 1) {
    // a Gaussian kernel width was given
    Nker = 256; // this is a reasonable default
    Klut = calloc(Nker, sizeof(*Klut));
    if (!Klut)
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                 "%s: Could not allocate kernel LUT (length %d)",
                 progname, Nker), errstr));
    // make_lut checks for negatives, NaNs
    if (!make_kernel_lut(Klut, Nker, Ktop, mxGetScalar(prhs[ARG_K]))) 
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                 "%s: Illegal spec for kernel LUT (check kparam, k(1))",
                 progname), errstr));
  } else {
    // explicit kernel given
    Klut = mxGetPr(prhs[ARG_K]);
    Nker = mxGetNumberOfElements(prhs[ARG_K]);
  }
  if (0)
    for (ok = 0; ok < Nker; ok++)
      printf("k[% 3d] = %.6f\n", ok, Klut[ok]);

  // KWt
  // if there's no KWt arg, make up one to fill in later
  if (nrhs > ARG_KWt)
    wt_holder = (mxArray *) prhs[ARG_KWt]; /* if uncast, the RHS is const */
  else
    wt_holder = mxCreateDoubleMatrix(0,0,mxREAL);

  // bin params
  if (nrhs > ARG_BWs) {
    bws = mxGetPr(prhs[ARG_BWs]);
    bwlen = mxGetM(prhs[ARG_BWs]);
    if (bwlen > 0 && mxGetN(prhs[ARG_BWs]) != 2)
      // (empty is OK)
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                 "%s: bws LUT should have 2 columns",
                 progname), errstr));
  } else {
    if (m < 2048) {
      // "small" image case
      bws   = Default_BWs_small;          // just a list of numbers
      bwlen = Default_BWs_count_small/2;  // always 2 columns
    } else {
      // "large" image case
      bws   = Default_BWs_large;          // just a list of numbers
      bwlen = Default_BWs_count_large/2;  // always 2 columns
    }
  }
  if ((msg = block_lut_validate(bws, bwlen)) != NULL)
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: bws LUT format error: %s",
               progname, msg), errstr));

  // output y
  plhs[ARG_Y] = mxCreateDoubleMatrix(m, n, mxREAL);  // always returned
  if (!plhs[ARG_Y])
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: Could not allocate y output",
               progname), errstr));
  // s, if we want it
  if (nlhs > ARG_S) {
    plhs[ARG_S] = mxCreateDoubleMatrix(m, n, mxREAL); // optionally returned
    if (!plhs[ARG_S])
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                 "%s: Could not allocate s output",
                 progname), errstr));
    s_data = mxGetPr(plhs[ARG_S]);
  } else {
    s_data = NULL;
  }
  // p, if we want it
  if (nlhs > ARG_P) {
    plhs[ARG_P] = mxCreateDoubleMatrix(m, n, mxREAL); // optionally returned
    if (!plhs[ARG_P])
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                 "%s: Could not allocate p output",
                 progname), errstr));
    p_data = mxGetPr(plhs[ARG_P]);
  } else {
    p_data = NULL;
  }

  /* call the function */
  ok = control(m, n,                 // sizes 
           mxGetPr(prhs[ARG_X]), // input is mXn
           mxGetPr(plhs[ARG_Y]), // convolved result is mXn 
           s_data,     // z coord, mXn (opt. returned if non-null)
           p_data,     // patch indexes, mXn (opt. returned if non-null)
           xcen,       // center, n or x, in C coords, origin at 0 
           ycen,       // center, m or y, in C coords, origin at 0 
           rsun,       // radius, in pixels 
           p0,         // p-angle, in radians
           Klut,       // kernel LUT
           Nker,       // and its size
           Ktop,       // lower range of the table 
           KswathI,    // box enclosing kernel 
           mxt_make_vector(wt_holder, KWt_count,Default_KWt, Default_KWt_count),
           bws,
           bwlen);
  // ok holds number of patches, or negative if trouble
  if (ok < 0)
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: Failed calloc, thread, or internal error in control(), code = %d",
               progname, ok), errstr));

  /* set up range of first output */
  // getrange(prhs[ARG_X], &datamin, &datamax);
  datamin = 0; datamax = 2; /* input is treated as an indicator */
  setrange(plhs[ARG_Y], datamin, datamax); 
  /* optional "s" output */
  if (nlhs > ARG_S) {
    setrange(plhs[ARG_S], 0.0, rsun); /* "z" index in pixels */
  } 
  /* optional "p" output */
  if (nlhs > ARG_P) {
    // ok is number of patches
    setrange(plhs[ARG_P], (double) Block_Min, (double) ok);
  } 
  // free up kernel lut if we generated it internally
  if (mxGetNumberOfElements(prhs[ARG_K]) == 1) 
    free(Klut);
  // free up weight matrix placeholder if it was used 
  if (nrhs <= ARG_KWt) mxDestroyArray(wt_holder);
  gettimeofday(&tv2, NULL);
  int dt = (((long) tv2.tv_sec) - ((long) tv1.tv_sec)) * 1000000 + 
    (((int) tv2.tv_usec) - ((int) tv1.tv_usec));
  if (verbose)
    mexPrintf("SS: dt = %.1f ms\n", dt/1000.0);
}



/* Hook for generic tail matter */
#ifdef MEX2C_TAIL_HOOK
#include "mex2c_tail.h"
#endif

---