proj/mag/ident/libs/mexfunctions/makemrfdiscwts.c

#include "mex.h"  /* must appear first */
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "mexhead.h" /* my mex defines */
#include "Doc/makemrfdiscwts_docstring.h"  /* autogenerated from this file */

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

%makemrfdiscwts make distance metric for mrf segmentation
% 
% dist=makemrfdiscwts(n,del,ctr,rho,mode,ell)
% * Given image sizes n, and neighborhood size del, make a distance 
% metric array dist, such that dist(:,m1,n1) is the distance-to-
% neighbors of the site m1,n1, where 1 <= m1 <= n(1) and 1 <= n1 <= n(2).
% * Size of del = 3 corresponds to a 3x3 neighborhood around each site.
% * The distance is computed for a disk of center ctr(1:2) and radius 
% ctr(3); ctr(1) corresponds to x or n1 and ctr(2) to y or m1.  
% ctr(1:2) = [1 1] corresponds to the first pixel.  
% * Thus, for instance, if Cx, Cy, r is the MDI center from mdidisk(), 
% we would specify ctr=[Cx Cy r].  Note that ctr(1) and n(2) correspond 
% to x, while ctr(2) and n(1) correspond to y.
% * Note that p-angle and b-angle are not needed by this routine.
% For simplicity, a 5-tuple containing these as entries 4 and 5 
% (a `geom' vector) can be given for ctr, and the tail will be ignored.
% * The distance scale factor is rho.  Zero corresponds to uniform
% distances (all one) and large positive values correspond to 
% separation-sensitive distances.  100 (20..200) is a typical value 
% for images of radius 500.  To gain more insight, the raw distances
% (theta below) are returned when rho=NaN is given.  A histogram
% of theta values relative to exp(-theta*rho) will show whether most
% theta's land in the linear part of the exponential.
% * If rho < 0 is given, the scale factor is taken to be abs(rho), but
% the distances are rescaled so that the overall max equals 1.  Otherwise,
% the distances diminish across the entire disk as rho increases.
% * More explicitly, consider two points (s, s') within the 2d disk.
% They map to points on the unit sphere, say
%   s = (x y z), s' = (x' y' z'), 
% The angle theta (>0) between them is found, and the distance 
% is returned as:
%   dist = exp(-theta*rho)
% Because theta is very close to zero, the code in fact uses the chord
% distance rather than the arc distance between the two points.
% * Because the distance is symmetric, dist only includes distances
% between sites s,s' where s' is less than s.  This is the case if 
% the pixel corresponding to s' comes before s in the memory footprint
% of an image.
% * Elliptical discs are allowed for.  In this case, ctr(3) is the 
% major axis semidiameter, and ell(1) is that of the minor axis.
% ell(2) is the counterclockwise rotation, in degrees, of the major
% axis off of the "m" or "y" axis.  If not given, ell(1)=ctr(3) and 
% ell(2)=0.  The implied condition ell(1) <= ctr(3) is not required 
% by the code.
% * If either s or s' is off-disk, zero is put in the corresponding
% distance value.  This is required by mrf_segment_wts.
% * The required mode string switches between sesw (mode = 'sesw')
% or transposed (mode = 'sene') pixel ordering, which works as follows.
% * The normal HMI (and normal MDI) pixel ordering starts in the
% southeast corner, and the first scan line of pixels runs toward the
% southwest corner.  This is `sesw' ordering.  The transposed ordering
% is `sene'; this ordering is what we used for the JPL MDI processing. 
% This is implemented via internal stride parameters.
% * This is implemented as a MEX file.  It agrees with makemrfdiscwts2.m
% 
% Inputs:
%   int n(1) or (2); -- if n is scalar, assume a square image
%   int del;
%   real ctr(3) or ctr(4) or ctr(5);
%   real rho;
%   string mode;
%   opt real ell(2) = [ctr(3) 0];
% 
% Outputs:
%   real dist ((del*del-1)/2,n(1),n(2));
% 
% See Also:  makemrfdiscwts2, mrf_segment_wts

% turmon oct 2006

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

/* standard boilerplate */
                   
#define NARGIN_MIN  5      /* min number of inputs */
#define NARGIN_MAX  6      /* max number of inputs */
#define NARGOUT_MIN 0      /* min number of output args */
#define NARGOUT_MAX 1      /* max number of output args */

#define ARG_n     0
#define ARG_del   1
#define ARG_ctr   2
#define ARG_rho   3
#define ARG_mode  4
#define ARG_ell   5

#define ARG_dist  0   // returned distance array

static const char *progname = "makemrfdiscwts";
#define PROGNAME makemrfdiscwts
static const char *in_specs[NARGIN_MAX] = {
  "RS|RV(2)",
  "RS(1)",
  "RV(3)|RV(4)|RV(5)",
  "RS(1)",
  "SV",
  "RV(2)" };
static const char *in_names[NARGIN_MAX] = {
  "n",
  "del",
  "ctr",
  "rho",
  "mode",
  "ell"};
static const char *out_names[NARGOUT_MAX] = {"dist"};

// defined just for brevity in a generated #include file
#define SHORTNAME mdw

/*
 * 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 ISLEGAL(l,L) ((l>=0)&&(l<L)) /* for indexing: 0 <= l < L */
#define ActiveLabel(y) ((!isnan(y)) && ((y) > 0)) /* not NaN and not 0 */
#define DIST_BLANK 0  /* distance value when s or s' is off-disk */

/*
 * makedist: compose distance array
 *
 * Regarding rotation, the viewpoint taken here is that a true disc
 * has been stretched and rotated by the imaging process.  Thus,
 * to find the distance of sites s and s', first their "true"
 * positions are recovered and then the distance is taken between
 * the two sites on an the sphere in which they are "truly" embedded.
 * 
 * In particular, we do not have the viewpoint that the two sites are
 * embedded in a 3-d ellipsoid which has been projected into a 2-d ellipse.
 * So, we do not take the view that the chord distance that we are finding
 * is a chord on an ellipsoid.  It is a chord on a sphere.  These two
 * viewpoints may be related in some easy way, or may be the same, but
 * it's not obvious to me and I have not checked.
 */
static
int
makedist(
     double *d,      /* distance output */
     double rho,
     int raw,        /* output raw distance only? */
     int rescale,    /* rescale d so max(*d) == 1 ? */
     double xcen,
     double ycen,
     double Rm,      /* radii */
     double Rn,
     double a0,      /* major-axis rotation -- NOT p-angle */
     int del,
     // switch from m,n to x,y here
     // strides have been in x,y -- but MRF code uses m,n
     // trying to be compatible with both
     int maxx,       // dims of eventual image (y/m)
     int maxy,       // x/n
     int strx,       // stride along x/n, units of doubles
     int stry        // stride along y/m
     )

{
  double *d0;         /* position in d list */
  int m, n, mp, np;   /* loop counters */
  int nu;             /* neighbor loop counter  */
  double x,  y,  z;   /* coords of s  */
  double xp, yp, zp;  /* coords of s' */
  double x0, y0;      /* temp variables */
  double theta;       /* angle from s -> s' */
  double filler;      /* blank fill value */
  double dmax;        /* max *d so far */
  const int delhalf = (del-1)/2;     /* neighbors on either side */
  const int del2half=(del*del-1)/2;  /* number of neighbors (~del^2/2) */
  const double RmInv = 1/Rm;
  const double RnInv = 1/Rn;
  const double Ca0 = cos(a0);
  const double Sa0 = sin(a0);
  
  /* setup */
  dmax = 0;  /* max so far is zero */
  if (rescale) 
    filler = mxt_getnand();  /* out of range d filled with NaN for now */
  else
    filler = DIST_BLANK;     /* out of range d get final fill value */
  /* 
   * loop over s 
   * note, m is like y, n is like x
   */
  for (m = 0; m < maxy; m++)
    for (n = 0; n < maxx; n++) {
      d0 = d + (n*strx+m*stry)*del2half; /* head of d list for (m,n) site */
      /* find s = x,y,z in unit coordinates */
      x0 = n - xcen; /* offsets from center */
      y0 = m - ycen;
      x = ( x0 * Ca0 + y0 * Sa0) * RnInv; /* rotate and scale */
      y = (-x0 * Sa0 + y0 * Ca0) * RmInv;
      z = x*x + y*y; /* intermediate result */
      /* off-disc check for s */
      if (z > 1) {
    for (nu = 0; nu < del2half; nu++)
      d0[nu] = DIST_BLANK;  /* set all nbrs to off-disk */
    continue;        /* i.e., next s = (m,n) */
      }
      z = sqrt(1 - z); /* this is the real z */
      /* (x,y,z) is now a position on the sphere where s lives */

      /* loop over s', keeping track of nbr count nu also */
      /* note: this is where the neighbor indexing is critical
       * eg, for del=3, we have 3x3 neighborhoods around the central
       * point [s], and del2half = 4 prior neighbors that are ordered 
       * as shown below.  The subsequent neighbors s' marked by [-] are 
       * not stored with s because s is a prior neighbor of each such
       * s'.  Storage is done this way so that every patch of distances
       * (ie, the numbered squares below) is still in the same indexing
       * order as the original images.  So, for example, if we have all 
       * 9 distances, a simple reshape(dists, [3 3]) will put them in
       * an immediately recognizable format.
       *     n->
       *  m  0 3 -
       *  |  1 s -
       *  V  2 - - 
       * Neighbors are indexed by nu = 0...del2half-1
       * This layout must match the layout expected by mrf_segment_wts.
       */
      for (np = -delhalf, nu = 0; np <= 0; np++)
    for (mp = -delhalf; mp <= delhalf; mp++, nu++) {
      if (nu >= del2half) break; /* no more neighbors */
      /* find s' = x',y',z' in unit coordinates */
      x0 = n + np - xcen;
      y0 = m + mp - ycen;
      xp = ( x0 * Ca0 + y0 * Sa0) * RnInv; /* rotate and scale */
      yp = (-x0 * Sa0 + y0 * Ca0) * RmInv;
      zp = xp*xp + yp*yp; /* intermediate result */
      /* off-disc check for s' */
      if (zp > 1) {
        d0[nu] = DIST_BLANK; /* set only the s' nbr to off-disk */
        continue;     /* i.e., next s' = (m+mp,n+np) */
      } 
      zp = sqrt(1 - zp); /* this is the real z */
      /* (xp,yp,zp) is now a position on the sphere where s' lives */
      /* find angle theta between s and s' */
      /* note, this approximates the arc length with the chord length */
      theta = sqrt((x-xp)*(x-xp) + (y-yp)*(y-yp) + (z-zp)*(z-zp));
      if (raw)
        d0[nu] = theta; /* note, ignore rho */
      else
        d0[nu] = exp(-rho*theta); /* usual case */
      if (d0[nu] > dmax) dmax = d0[nu]; /* update the max */
    }
    }
  /* 
   *  Rescale if that was desired
   */
  if (rescale) {
    dmax = 1.0/dmax; /* invert so the inner loop is a multiply */
    for (nu = 0; nu < maxx*maxy*del2half; nu++)
      if (isnan(d[nu]))
    d[nu] = DIST_BLANK; /* replace NaN sentinel with blank value */
      else
    d[nu] *= dmax; /* a valid distance: rescale it */
  }
  return 1;
}
   


/*
 * 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 disc 
  int maxx, maxy;      // sizes, transposed if needed
  int strx, stry;      // strides in images along x(n) and y(m) dims
  int del;             // nbhd full width 
  int dims[3];         // space for 3 dimension lengths 
  double rho;          // rho value which is used for flags too 
  double Rm, Rn;       // radii 
  double a0;           // major axis rotation angle
  int mode_sesw;       // sesw = 1, transposed = 0
  char *mode, *word;   // label for converted mode string
  int ok;
  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);

  /* n, 1x1 or 1x2 */
  M = (int) (mxGetPr(prhs[ARG_n])[0]);
  if (mxGetNumberOfElements(prhs[ARG_n]) < 2)
    N = M;
  else 
    N = (int) (mxGetPr(prhs[ARG_n])[1]);
  /* neighborhood size del */
  del = mxGetScalar(prhs[ARG_del]);

  /* 
   * `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
   */
  dims[0] = (del*del-1)/2;
  dims[1] = M;
  dims[2] = N;
  /* plhs[ARG_dist] = mxCreateDoubleMatrix(M * (del*del-1)/2, N, mxREAL); */
  plhs[ARG_dist] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
  if (!plhs[ARG_dist])
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: mxCreate call: failed for dist",
               progname), errstr));
  rho = mxGetPr(prhs[ARG_rho])[0]; /* output strength, scaling, type */

  /*
   * ellipticity params
   */
  // major axis length
  Rm = mxGetPr(prhs[ARG_ctr])[2];  
  // ellipticity, if any
  if (nrhs <= ARG_ell) {
    // no ellipticity given: make it a circle
    Rn = Rm;  // n or x scale same as m or y scale
    a0 = 0.0; // rotation angle is arbitrary in this case
  } else {
    // elliptical
    Rn = mxGetPr(prhs[ARG_ell])[0]; // length given explicitly
    a0 = mxGetPr(prhs[ARG_ell])[1]*(-1*M_PI/180.0); // deg->radians
  }

  /*
   * 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 geom input, 
   * 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
  }

  /*
   * do the computation
   */
  ok = makedist(
        mxGetPr(plhs[ARG_dist]),
        fabs(rho),       // allow for rho < 0, don't care if NaN
        isnan(rho),      // NaN means output raw chord lengths
        (int) (rho < 0), // rho < 0: rescale d so max(*d) = 1
        // disc center 
        mxGetPr(prhs[ARG_ctr])[0]-1.0, // ctr n/x, C origin 
        mxGetPr(prhs[ARG_ctr])[1]-1.0, // ctr m/y, C origin 
        // disc radii
        Rm, Rn,
        // major axis orientation
        a0,
        // sizes as ints
        del, maxx, maxy, strx, stry);
  if (!ok) {
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
             "%s: makedist() call: failed",
              progname), errstr));
  }
}



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

---