proj/mag/harp/libs/matlab/mfile-mex/standalone/roi_stats_mag.c

#include "mex.h"  /* must appear first */
#include <stdio.h>
#include <strings.h>
#include <string.h>
#include <math.h>
#include "mexhead.h" /* my mex defines */
#include "Doc/roi_stats_mag_docstring.h"  /* autogenerated */


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

%roi_stats_mag: accumulate statistics on regions
% 
% [s,names,combo]=roi_stats_mag(x,y,mag,geom,nroi,mode)
% * A set of per-region statistics is gathered based on regions 
% encoded in the image inputs x (containing region tags, 1..Nr, or
% 0 for no tag) and y (containing region indicators, 0/1).  The
% statistics are functions of region configuration, and of 
% line-of-sight magnetic field `mag'.
% * x should be 0/NaN, or a nonnegative integer.  Pixels in the 
% range 1..Nr, where Nr is the number of regions, are treated as
% indicators for being within the numbered region.
% Inputs x of 0 or NaN are treated as not belonging to any
% labeled region.
% * The number of regions, Nr, is deduced from x.  Often you want
% the statistics to cover a known range of classes regardless of x.  
% If so, specify the nroi input (an integer).  If you give nroi < 0,
% or empty, it is taken as equal to Nr.  If nroi < Nr, pixels with 
% x > nroi are ignored.
% * y that is finite and nonzero (e.g., 1) indicates activity present.  
% Not all tagged pixels will be active; typically the tagged pixels 
% are large blobs, and the active pixels are finer details within 
% each blob.
% * For each region, a row of statistics is computed:
%    1: rgnnum = # pixels tagged (all may not be active)
%    2: rgnsize = projected (flat) area in microhemispheres (0..1e6)
%    3: rgnarea = un-projected (solid-angle) area in microhemispheres (0..1e6)
%    4: arnum = # active pixels (x == 1)
%    5: arsize = projected (flat) active area in microhemispheres (0..1e6)
%    6: ararea = unprojected (solid-angle) active area in microhemis (0..1e6)
%    7: arlat,lon = average location, weighted by un-projected area
%    9: arminlat,lon  lower corner of (lat,lon) bounding box
%   11: armaxlat,lon  upper corner of (lat,lon) bounding box
%   13: rgnbtot = sum of absolute LoS flux within the identified region
%   14: rgnbnet = net LoS flux within the identified region
%   15: rgnbpos = absolute value of total positive LoS flux
%   16: rgnbneg = absolute value of total negative LoS flux
%   17: arfwtlat,lon = flux-weighted center of active pixels
%   19: arfwtpos_lat,lon = flux-weighted center of positive flux
%   21: arfwtneg_lat,lon = flux-weighted center of negative flux
%   23: daysgone = #days until bounding box vanishes from front of disk
%   24: daysback = #days until bounding box first reappears on front
% These definitions are in a machine-readable include file,
% roi_stats_mag_def.h
% * The solar disk is at location given by center, observed at tip 
% angle beta, in degrees.
% * Optionally returns the standard text name of each of the
% computed statistics. 
% * Optionally returns the means of combination of each of the 
% computed statistics. 
% * This is implemented as a MEX file.
% 
% Inputs:
%   real x(m,n);   -- 1..Nr, or 0/NaN
%   int  y(m,n);   -- 0/NaN, or otherwise
%   real mag(m,n);
%   real geom(5);  -- [x0 y0 r_sun b p]
%   int nroi;      -- [] same as -1
%   string mode;
% 
% Outputs:
%   real s(nr,ns)  -- nr = (nroi < 0) ? Nr : nroi
%   opt string names(ns)
%   opt string combo(ns)
% 
% See Also:

% turmon june 2010, sep 2010

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


/* standard boilerplate */
static const char *progname = "roi_stats_mag";
#define PROGNAME roi_stats_mag
#define SHORTNAME rsm

#define NARGIN_MIN  6      /* 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_Y     1
#define ARG_MAG   2
#define ARG_GEOM  3
#define ARG_NROI  4
#define ARG_MODE  5

#define ARG_S     0
#define ARG_NAMES 1
#define ARG_COMBO 2

static const char *in_specs[NARGIN_MAX] = {
  "RM",
  "RM",
  "RM",
  "RV(5)",
  "IS",
  "SV"};
static const char *in_names[NARGIN_MAX] = {
  "x",
  "y",
  "mag",
  "geom",
  "nroi",
  "mode"};
static const char *out_names[NARGOUT_MAX] = {
  "s",
  "names",
  "combo"};

#include "roi_stats_mag_defs.h"

/*
 * Solar rotation constants
 * all in deg/day
 * first has sidereal rate, less Earth's rotational rate, to make synodic rate
 * if you're sensitive to these values, you're probably doing something wrong...
 */

#define ROTATE_CONST_K0 (14.6-0.9865) 
#define ROTATE_CONST_K1 (-2.2)
#define ROTATE_CONST_K2 0.0


/*************************************************************
 *
 * HELPER ROUTINES
 *
 *************************************************************/


/*
 * extracts various (OK, 1) mode parameters from a free-form string
 *
 * returns 0 if the string was OK, an explanatory tag otherwise
 */
static
char *
extract_mode(char *mode, int *Msesw) 
{
  char *word;
  char *sep = ", \t";

  if (!mode) return "null"; /* should not happen */
  for (word = strtok(mode, sep); word; word = strtok(NULL, sep)) {
    /* sesw, or transposed */
    if      (strcasecmp(word, "sesw"  ) == 0) *Msesw = 1;
    else if (strcasecmp(word, "sene"  ) == 0) *Msesw = 0;
    /* unrecognized word: not OK */
    else return word; 
  }
  return 0; /* OK */
}

/* simple utilities */ 
#define WithinROI(y)   ((!isnan(y)) && ((y) > 0)) /* not NaN and not 0 */
#define ActiveLabel(y) ((!isnan(y)) && ((y) > 0)) /* not NaN and not 0 */

/*
 * stat_compute: calculate statistics for each region
 */
static
void
stat_compute(double *roi,   // input roi image ("x")
         double *ar,    // input acr image ("y")
         double *mag,   // input mgram image, or NULL
         double **s,    // result: statistics per region (2d)
         int Nr,        // number of regions
         int maxx,      // image size 
         int maxy,
         int strx,      // stride along x, units of doubles
         int stry,      // stride along y
         // geom
         double x0,     // x center, in C coords, origin at 0 
         double y0,     // y center, in C coords, origin at 0 
         double R,      // radius, in pixels 
         double B,      // tilt angle, radians 
         double P)      // twist angle, radians

{
  int x, y;               // loop counters
  double r_dbl;           // region number as a double
  int r;                  // region number as an int
  int offset;             // offset into arrays
  double lat, lon;        // hold current lat,lon
  double p1x, p1y, p1z;   // point P1 coords in 3D
  double p2x, p2y, p2z;   // point P2 coords in 3D
  double temp;            // workspace slot
  double area1;           // un-projected area
  double tau;             // a threshold for finding the limb
  double mag1, magM1;     // current mgram value
  const double nan = mxt_getnand(); // cache nan
  const double cosB = cos(B);
  const double sinB = sin(B);
  const double cosP = cos(P);
  const double sinP = sin(P);
  const double Rinv = 1/R;

  /* 
   * STEP 0: set initial values for accumulators 
   */
  // most stuff starts at zero, so do it all now
  bzero(&(s[0][0]), ((size_t)RS_num_stats)*Nr*sizeof(s[0][0]));
  for (r = 0; r < Nr; r++) {
    /* (first, rgn stuff) */
    // size
    // 0: rgn_num, rgn_size, rgn_area
    // extent
    // set greater than +2pi: rgn_minlat, rgn_minlon
    // set lower than -2pi: rgn_maxlat, rgn_maxlon
    s[RS_rgn_min_lat][r] = s[RS_rgn_min_lon][r] =  10;
    s[RS_rgn_max_lat][r] = s[RS_rgn_max_lon][r] = -10;
    // return time
    // 0: daysgone, daysback (they are not accumulators)
    // flux
    // 0: rgn_btot, rgn_bnet, rgn_bpos, rgn_bneg,    
    // flux moments
    // 0: rgn_bsum1, 2, 3, 4
    // flux moments, standardized
    // 0: rgn_bmean, rgn_bsdev, rgn_bskew, rgn_bkurt
    /* (next, ar stuff) */
    // size
    // 0: ar_num,  ar_size,  ar_area     
    // flux
    // 0: ar_btot, ar_bnet, ar_bpos, ar_bneg,    
    // mean location
    // 0: ar_area_lat,  ar_area_lon
    // 0: ar_fwt_lat,    ar_fwt_lon
    // 0: ar_fwtpos_lat, ar_fwtpos_lon
    // 0: ar_fwtneg_lat, ar_fwtneg_lon
  }

  /* 
   * STEP 1: loop over all pixels and accumulate region info 
   */
  for (x = 0; x < maxx; x++) 
    for (y = 0; y < maxy; y++) {
      offset = x*strx + y*stry;  // still in units of doubles
      r_dbl = roi[offset]; // a double
      if (!WithinROI(r_dbl)) continue; // ensure the NaN check works
      r = (int) floor(r_dbl);
      if (r > Nr || r < 1) continue;  // ensure it's in-range
      r--; /* make region number (1..Nr) into a region index (0..Nr-1) */
      mag1  = mag ? mag[offset] : nan; // abbreviate the magnetic field...
      magM1 = fabs(mag1); // ...and its magnitude
      /*
       * Find the coordinates, including lat/lon
       */
      /* find original vector coordinates P2 = (p2x, p2y, p2z) */
      p2x = x - x0;
      p2y = y - y0;
      temp = (R + p2x) * (R - p2x) - p2y*p2y; /* R*R-x*x = (R-x)*(R+x) */
      /* temp is the residual which goes into the p2z component */
      if (temp <= 0) continue; /* off-disk (should not happen) */
      p2z = sqrt(temp);   /* totally OK -- on visible side */
      area1 = R/p2z;      /* note: the original, unrotated, 1/z */
      if (area1 > 1e4) area1 = 1e4; // cap it -- could do better here
      /* rotate by p in the xy plane: P1 = rot(P2) (the p-angle) */
      p1x =  cosP * p2x + sinP * p2y;
      p1y = -sinP * p2x + cosP * p2y;
      p1z =  p2z;
      /* rotate by +beta in yz plane (x axis fixed): P2 = rot(P1) */
      p2x =  p1x;  /* no change to x */
      p2y =  cosB * p1y + sinB * p1z;
      p2z = -sinB * p1y + cosB * p1z;
      /* find (lat,lon) of the point */
      lon = atan2(p2x, p2z);  /* p2=(1,0,0), at limb, has lon=pi/2 */
      lat = asin(p2y * Rinv);

      /*
       * Update all the accumulators 
       * (listed in order)
       * RS_rgn_* are whole-region accumulators (all over the ROI)
       * RS_ar_*  are ar-only accumulators
       */
      // basic sizes
      s[RS_rgn_num ][r]++;               // count pixels
      s[RS_rgn_size][r]++;               // count projected area
      s[RS_rgn_area][r] += area1;        // unprojected area, proportional to 1/z
      // lat/lon bounding box, whole-region
      if (lat < s[RS_rgn_min_lat][r]) s[RS_rgn_min_lat][r] = lat;
      if (lat > s[RS_rgn_max_lat][r]) s[RS_rgn_max_lat][r] = lat;
      if (lon < s[RS_rgn_min_lon][r]) s[RS_rgn_min_lon][r] = lon;
      if (lon > s[RS_rgn_max_lon][r]) s[RS_rgn_max_lon][r] = lon;
      // (return time is not found by accumulation)
      // line-of-sight flux, whole-region
      s[RS_rgn_btot][r] += magM1;
      s[RS_rgn_bnet][r] += mag1;
      s[RS_rgn_bpos][r] += (mag1 > 0) ? magM1 : 0.0;
      s[RS_rgn_bneg][r] += (mag1 < 0) ? magM1 : 0.0;
      // line-of-sight flux, moments
      // TODO: check if number-of-pixels normalization is correct
      s[RS_rgn_bsum1][r] += mag1;
      s[RS_rgn_bsum2][r] += mag1*mag1;
      s[RS_rgn_bsum3][r] += mag1*mag1*mag1;
      s[RS_rgn_bsum4][r] += mag1*mag1*mag1*mag1;
      if (ActiveLabel(ar[offset])) {
    // basic sizes
    s[RS_ar_num ][r]++;              // count pixels
    s[RS_ar_size][r]++;              // count projected area
    s[RS_ar_area][r] += area1;       // unprojected area, proportional to 1/z
    // line-of-sight flux, ar-only
    s[RS_ar_btot][r] += magM1;
    s[RS_ar_bnet][r] += mag1;
    s[RS_ar_bpos][r] += (mag1 > 0) ? magM1 : 0.0;
    s[RS_ar_bneg][r] += (mag1 < 0) ? magM1 : 0.0;
    // ar-only area-weighted lat/lon
    s[RS_ar_area_lat][r] += lat * area1;
    s[RS_ar_area_lon][r] += lon * area1;
    // flux-weighted lat/lon
    s[RS_ar_fwt_lat][r] += lat * magM1;
    s[RS_ar_fwt_lon][r] += lon * magM1;
    // flux-weighted centers of + and - flux
    s[RS_ar_fwtpos_lat][r] += (mag1 > 0) ? lat*magM1 : 0.0;
    s[RS_ar_fwtpos_lon][r] += (mag1 > 0) ? lon*magM1 : 0.0;
    s[RS_ar_fwtneg_lat][r] += (mag1 < 0) ? lat*magM1 : 0.0;
    s[RS_ar_fwtneg_lon][r] += (mag1 < 0) ? lon*magM1 : 0.0;
      }
    } /* for (x,y) */

  /* 
   * STEP 2: normalize and do unit conversions on s array 
   */
  // 2a: convert accumulated sums for mean locations into averages
  //     in order of definition
  for (r = 0; r < Nr; r++) {
    // average location, weighted by un-projected area
    s[RS_ar_area_lat  ][r] /= s[RS_ar_area][r];  // the sum of the weights 1/z
    s[RS_ar_area_lon  ][r] /= s[RS_ar_area][r];
    // flux-weighted center of active pixels
    s[RS_ar_fwt_lat   ][r] /= s[RS_ar_btot][r];
    s[RS_ar_fwt_lon   ][r] /= s[RS_ar_btot][r];
    // flux-weighted center of positive flux
    s[RS_ar_fwtpos_lat][r] /= s[RS_ar_bpos][r];
    s[RS_ar_fwtpos_lon][r] /= s[RS_ar_bpos][r];
    // flux-weighted center of negative flux
    s[RS_ar_fwtneg_lat][r] /= s[RS_ar_bneg][r];
    s[RS_ar_fwtneg_lon][r] /= s[RS_ar_bneg][r];
  }
  // 2b: compute (daysgone, daysback) from lat-lon info.  
  /* Note, the visible disk is characterized by: 
   *    cos(lon) >= tan(beta) * tan(lat).  
   * Thus, find the values dt for which: 
   *    lon + (dt)*(vel) 
   * satisfies this. The above condition is independent of P.
   */
  for (r = 0; r < Nr; r++) {
    double ar_lat; // the latitude to use for the AR, in radians

    if (s[RS_ar_num][r] > 0)
      // use the flux-weighted AR-pixel lat, still in radians
      ar_lat = s[RS_ar_fwt_lat][r]; 
    else
      // fall back to the midpoint of the lat of the region
      ar_lat = (s[RS_rgn_min_lat][r] + s[RS_rgn_max_lat][r]) * 0.5;
    // find threshold for longitude (using average lat, still in radians)
    tau = (sinB/cosB) * tan(ar_lat);
    if (fabs(tau) >= 1) continue;  // beta, lat: all rotations visible!
    tau = acos(tau);     /* in the range 0..pi */
    // find angular velocity at our latitude (still in radians)
    temp = sin(ar_lat);              // sin of average latitude
    temp *= temp;                    // temp = sin(lat)^2
    temp =    ROTATE_CONST_K0 + 
      temp * (ROTATE_CONST_K1 + 
      temp * (ROTATE_CONST_K2));     // k0 + k1*temp + k2*temp^2, in deg/day
    temp *= (M_PI/180.0);            // convert deg/day -> radians/day
    // translate temp, an angular rate in radians/day, into a time in days
    // minlon: used for last passage behind limb
    // maxlon: used for first reappearance
    s[RS_rgn_daysgone][r] = (tau -            s[RS_rgn_min_lon][r])/temp;  
    s[RS_rgn_daysback][r] = ((2*M_PI - tau) - s[RS_rgn_max_lon][r])/temp; 
  }
  // 2c: compute (mean,sdev,skew,kurtosis) from moments
  for (r = 0; r < Nr; r++) {
    double N;        // normalizing constant
    double mu, sig2; // first moment, and second central moment (variance)
    double M_3, M_4; // third and fourth central moments

    // using number-of-pixels normalization, not area
    N    =  s[RS_rgn_num  ][r]; // currently, number-of-pixels
    mu   =  s[RS_rgn_bsum1][r] / N;
    sig2 = (s[RS_rgn_bsum2][r] / N) - mu*mu;
    if (sig2 < 0) sig2 = 0;  // numerical trouble could cause this
    if (N < 2) sig2 = nan;   // cannot compute sdev in this case
    // third central moment, where S_i is the i'th noncentral moment:
    //  (use binomial row [1 -3 3 -1]; last two entries combine into +2)
    // E (x-mu)^3 = S_3 - 3 mu S_2 + 2 mu^3
    M_3 = (s[RS_rgn_bsum3][r]/N) - 3*mu*(s[RS_rgn_bsum2][r]/N) + 2*mu*mu*mu;
    // fourth central moment:
    //  (use binomial row [1 -4 +6 -4 1]; last two entries combine into -3)
    // E (x-mu)^4 = S_4 - 4 mu S_3 + 6 mu^2 S_2 - 3 mu^4
    M_4 = (s[RS_rgn_bsum4][r]/N) - 4*mu*(s[RS_rgn_bsum3][r]/N) +
      6*mu*mu*(s[RS_rgn_bsum2][r]/N) - 3*mu*mu*mu*mu;

    s[RS_rgn_bmean][r] = mu;
    s[RS_rgn_bsdev][r] = sqrt(sig2);
    s[RS_rgn_bskew][r] = M_3 / (sig2 * sqrt(sig2));
    s[RS_rgn_bkurt][r] = (M_4 / (sig2*sig2)) - 3.0;
  }
  // 2d: unit conversions, in order of definition
  for (r = 0; r < Nr; r++) {
    const double rad2deg = 180.0/M_PI;
    const double size2uhemi = 1e6/(M_PI*R*R);   // 1 hemi has pi R^2 pixels
    const double area2uhemi = 1e6/(2*M_PI*R*R); // 1 hemi integrates to 2piR^2
    
    // area/sizes to microhemispheres
    // no conversion: RS_rgn_num
    s[RS_rgn_size][r] *= size2uhemi;
    s[RS_rgn_area][r] *= area2uhemi;
    // angles to degrees
    s[RS_rgn_min_lat][r] *= rad2deg;  
    s[RS_rgn_min_lon][r] *= rad2deg;
    s[RS_rgn_max_lat][r] *= rad2deg;  
    s[RS_rgn_max_lon][r] *= rad2deg;
    // no conversion: daysgone, daysback
    // no conversion: fluxes
    // no conversion: RS_ar_num
    s[RS_ar_size][r] *= size2uhemi;
    s[RS_ar_area][r] *= area2uhemi;
    // no conversion: fluxes
    // angles to degrees
    s[RS_ar_area_lat  ][r] *= rad2deg;  
    s[RS_ar_area_lon  ][r] *= rad2deg;
    s[RS_ar_fwt_lat   ][r] *= rad2deg;  
    s[RS_ar_fwt_lon   ][r] *= rad2deg;
    s[RS_ar_fwtpos_lat][r] *= rad2deg;  
    s[RS_ar_fwtpos_lon][r] *= rad2deg;
    s[RS_ar_fwtneg_lat][r] *= rad2deg;  
    s[RS_ar_fwtneg_lon][r] *= rad2deg;
  }
  // 2d: ensure NaNs for empty regions
  //     (Nr may be larger than range of roi array)
  for (r = 0; r < Nr; r++) {
    if (s[RS_rgn_num][r] == 0) {
      s[RS_rgn_min_lat ][r] = s[RS_rgn_min_lon ][r] = nan;
      s[RS_rgn_max_lat ][r] = s[RS_rgn_max_lon ][r] = nan;
      s[RS_rgn_daysgone][r] = s[RS_rgn_daysback][r] = nan;
    }
    if (s[RS_ar_num][r] == 0) {
      s[RS_ar_area_lat  ][r] = s[RS_ar_area_lon  ][r] = nan;
      s[RS_ar_fwt_lat   ][r] = s[RS_ar_fwt_lon   ][r] = nan;
      s[RS_ar_fwtpos_lat][r] = s[RS_ar_fwtpos_lon][r] = nan;
      s[RS_ar_fwtneg_lat][r] = s[RS_ar_fwtneg_lon][r] = nan;
    }
  }
}
   
/*
 * Are classes all integers, 0, or NaN?
 */
static
int
count_classes(double *y,
          int N)
{
  double Nr = 0; /* class count */
  int Nri;
  int n; 

  for (n = 0; n < N; n++)
    if (!isnan(y[n])) {
      if (y[n] < 0 || floor(y[n]) != y[n])
    return -1;
      if (y[n] > Nr)
    Nr = y[n];
    }
  Nri = (int) Nr;
  if (Nri != Nr)
    return -1; // trouble
  else
    return Nri;
}


/*************************************************************
 *
 * DRIVER ROUTINE
 *
 *************************************************************/

/*
 * Gateway routine
 */
#ifdef StaticP  /* undefined under mex */
StaticP
#endif
void 
mexFunction(
        int nlhs, 
        mxArray *plhs[], 
        int nrhs,
        const mxArray *prhs[])
{
  int m,  n;                 // size of x
  int maxx, maxy;            // sizes, transposed if needed
  int strx, stry;            // strides in images along x and y dims
  int Nr;                    // number of regions found in x
  int Nroi;                  // number of regions in the output
  int mode_sesw;             // sesw = 1, transposed = 0
  char *mode, *word;         // label for converted mode string and substring
  double *g;                 // convenience
  double x0, y0, r, b, p;    // geom
  double *mag;
  double **stat2;
  char errstr[256];

  /* 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 args
   */
  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);
  start_sizechecking();
  // size(x) == size(y)
  sizeinit(prhs[ARG_Y]);
  sizeagree(prhs[ARG_X]); 
  sizecheck_msg(progname, in_names, ARG_Y);
  // if size(mag) > 0, size(mag) == size(x)
  if (mxGetNumberOfElements(prhs[ARG_MAG]) > 0) {
    sizeinit(prhs[ARG_MAG]);
    sizeagree(prhs[ARG_X]); 
    sizecheck_msg(progname, in_names, ARG_MAG);
  }

  /* get size of input image x, needed below */
  m = (int) mxGetM(prhs[ARG_X]);
  n = (int) mxGetN(prhs[ARG_X]);

  /******************************************************************
   * Read and check arguments in order supplied
   ******************************************************************/

  /*
   * image geometry
   */
  // source image (x) geometry
  g   = mxGetPr(prhs[ARG_GEOM]); // abbreviate geom
  x0 = g[0] - 1;         // x0, in C coordinates
  y0 = g[1] - 1;         // y0, in C coordinates
  r  = g[2];             // disk radius
  b  = g[3]*M_PI/180.0;  // deg -> rad
  p  = g[4]*M_PI/180.0;  // deg -> rad

  /* 
   * `nroi' input
   * (required, can be empty by IS rules)
   */
  if (mxGetNumberOfElements(prhs[ARG_NROI])) 
    Nroi = mxGetScalar(prhs[ARG_NROI]);
  else
    Nroi = -1;

  /* 
   * `mode' input
   * (required)
   */
  if ((mode = mxArrayToString(prhs[ARG_MODE])) == NULL)
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: bad mode (non-string?). "
               "Could not convert mode arg to string.", 
               progname), errstr));
  /* no default values -- user must specify both explicitly */
  mode_sesw = -1; 
  if ((word = extract_mode(mode, &mode_sesw)) != NULL)
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: bad word <%s> in mode <%.80s>", 
               progname, word, /* need a fresh copy of the mode */
               mxArrayToString(prhs[ARG_MODE])), errstr));
  if (mode_sesw < 0)
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: need sesw/sene in mode <%.80s>", 
               progname, /* need a fresh copy of the mode */
               mxArrayToString(prhs[ARG_MODE])), errstr));
  mxFree(mode); // done with the mode string

  /*
   * create space for output
   */
  /* 1: find how big it must be */
  m  = mxGetM(prhs[ARG_X]);
  n  = mxGetN(prhs[ARG_X]);

  /* 2: make the space */
  // first, check x and find #regions
  if ((Nr = count_classes(mxGetPr(prhs[ARG_X]), m*n)) < 0)
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: Regions x must be NaN or integers >= 0", 
               progname), errstr));
  // this Nroi lets range of x determine the number of regions to use
  if (Nroi < 0)
    Nroi = Nr;
  // create space for the region info
  plhs[ARG_S] = mxCreateDoubleMatrix(Nroi, (mwSize) RS_num_stats, mxREAL); 
  if (!plhs[ARG_S])
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: Failed to allocate space for s output", 
               progname), errstr));

  /*
   * Depending on pixel-order mode, set up strides and bounding boxes in
   * (x,y) coordinates from their (m,n) settings.
   * Note that the (x0,y0) positions, in the geom1 and geom2 inputs, 
   * correspond to x and to y, so they do not need to be swapped.
   */
  if (mode_sesw) {
    // x (E-W axis) varies fastest -- x corresponds to matlab m coordinate
    strx = 1; stry = m; // x varies fastest, x has m points...
    maxx = m; maxy = n; // y has n points
  } else {
    // y (N-S axis) varies fastest -- x corresponds to matlab n coordinate
    strx = m; stry = 1; // y varies fastest, y has m points...
    maxx = n; maxy = m; // and x has n points
  }

  // simplify indexing into stat
  stat2 = mxt_make_matrix2(plhs[ARG_S], -1, -1, 0.0);
  if (!stat2)
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: Failed to allocate space for s pointers", 
               progname), errstr));
  // account for empty mag case
  mag = (mxGetNumberOfElements(prhs[ARG_MAG]) > 0) ? 
    mxGetPr(prhs[ARG_MAG]) : NULL;

  /*
   * do the computation
   */
  stat_compute(mxGetPr(prhs[ARG_X]), 
           mxGetPr(prhs[ARG_Y]), 
           mag,
           stat2,
           Nroi,                  // desired number of regions
           maxx, maxy,            // size
           strx, stry,            // strides in units of doubles
           x0, y0, r, b, p);      // geom
  // for debugging 
  if (0) {
    int r, sn;
    double *stats = mxGetPr(plhs[ARG_S]);

    for (sn = 0; sn < RS_num_stats; sn++)
      printf("\t%s", RS_index2name[sn]);
    printf("\n");
    for (r = 0; r < Nroi; r++) {
      printf("P%d: ", r);
      for (sn = 0; sn < RS_num_stats; sn++)
    printf("\t%.3g", stats[sn*Nroi+r]);
      printf("\n");
    }
  }
  // free up this indexing tool
  mxFree(stat2);

  /*
   * plug in names if wanted
   */
  if (nlhs > ARG_NAMES) {
    plhs[ARG_NAMES] = mxCreateCharMatrixFromStrings((int) RS_num_stats, 
                            RS_index2name);
    if (plhs[ARG_NAMES] == NULL)
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                 "%s: Failed to create `%s' output for %d names",
                 progname, out_names[ARG_NAMES], 
                 (int) RS_num_stats), 
            errstr));
  }

  /*
   * plug in combo if wanted
   */
  if (nlhs > ARG_COMBO) {
    plhs[ARG_COMBO] = mxCreateCharMatrixFromStrings((int) RS_num_stats, 
                            RS_index2combo);
    if (plhs[ARG_COMBO] == NULL)
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                 "%s: Failed to create `%s' output for %d combo-specs",
                 progname, out_names[ARG_COMBO], 
                 (int) RS_num_stats), 
            errstr));
  }
}


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

---