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

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

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

% clean_edge_label: clean up labels at extreme edge of a disk
% 
% [res,nclean]=clean_edge_label(img,center,delta,fill,mode);
% * Remove possibly-tainted labels at edge of disk by cleaning
% pixels in an annulus of width delta pixels.  Also, sets all 
% off-disk values to the value of the (1,1) pixel of img.
% * There are two cleaning modes.  Easiest is, set values in
% the edge annulus to `fill' (NaN is OK).  But, if `mode' 
% contains the string "adaptive", we instead propagate the value
% just inside the annulus radially outward to the limb, and the
% `fill' value is unused.
% * The number of on-disk pixels altered is in nclean.
% * Primary disk parameters (x-center, y-center, radius) are in 
% `center'.  The sun is the disk defined by these three numbers.  
% As a convenience, 'center' can be a "geom" vector, a 5-tuple
% with beta and p-angle at the end (these are not needed by
% this routine).
% * If delta < 0, img is propagated to res without further ado.
% (The off-disk clearing is not done either.)
% In this case, the center parameter need not be correct.
% * The "mode" string also switches between sesw (mode = 'sesw')
% or transposed (mode = 'sene') pixel ordering.  You must specify
% one or the other.
% * 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.
% 
% 
% Inputs:
%   real img(m,n);
%   real center(3) or (5);
%   real delta;
%   real fill;
%   string mode;
% 
% Outputs:
%   real res(m,n);
%   int nclean;
% 

% implemented as a mex file

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

/* constants and globals used for error checking */

#define NARGIN_MAX  5        /* number of inputs */
#define NARGIN_MIN  5        /* number of inputs */
#define NARGOUT_MIN 1    /* number of output args */
#define NARGOUT_MAX 2    /* number of output args */

#define ARG_img      0  // input image
#define ARG_center   1  // apparent center and radius of the sun
#define ARG_delta    2  // width of annulus
#define ARG_fill     3  // infill value
#define ARG_mode     4  // mode of operation

#define ARG_out      0  /* output image */
#define ARG_nclean   1  /* pixel change count */

static const char *progname = "clean_edge_label";
#define PROGNAME clean_edge_label
static const char *in_specs[NARGIN_MAX] = {
  "RM",
  "RV(3)|RV(5)",
  "RS",
  "RS",
  "SV" };
static const char *in_names[NARGIN_MAX] = {
  "image",
  "center",
  "delta",
  "fill",
  "mode"};
static const char *out_names[NARGOUT_MAX] = {
  "res",
  "nclean"};

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

/*
 * 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, int *Madapt) 
{
  char *word, *word2, *endptr;
  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;
    /* adaptive? */
    else if (strcasecmp(word, "adaptive") == 0) *Madapt = 1;
    /* unrecognized word: not OK */
    else return word;
  }
  return 0; /* OK */
}


/* computational routine: 
 * identify and fix edge labels
 * Formerly: no non-edge labels are not propagated from image -> imagep
 * New Oct 2010: off-disk pixels are all set to image[0], and all
 * pixels are propagated here.
 */

static
void
propagate(int *count,         // count of changed pixels 
      double *image,      // original image  
      double *imagep,     // updated image 
      int maxx, int maxy, // image, imagep size 
      int strx,           // stride along x, units of doubles
      int stry,           // stride along y
      double cenx,        // image center 
      double ceny, 
      double radius,      // image radius 
      double delta,       // annulus width 
      double fill_val,    // infill value
      int mode_adapt
      )
{
  int x, y;          // loop vars
  int offset, inset; // offsets into image arrays
  double r2;         // radius^2 at (x,y)
  double dx, dy;     // offset from center to (x,y)
  int x_off, y_off;  // integer indexes for fill-in value
  int ct = 0;        // local count for return value
  double scale;
  double offdisk;
  const double radius2 = radius*radius;
  const double radius_inside2 = (radius - delta)*(radius - delta);

  // printf("%s: xcen, ycen = (%.3f, %.3f)\n", progname, cenx, ceny);
  if (maxx * maxy > 0)
    offdisk = image[0]; // typically, NaN
  else
    offdisk = 0; // will not be used
  for (x = 0; x < maxx; x++)
    for (y = 0; y < maxy; y++) {
      offset = x*strx + y*stry;  // still in units of doubles
      // distances from center
      dx = x - cenx;
      dy = y - ceny;
      // current point's radius-squared
      r2 = dx*dx + dy*dy;
      // case on radius
      if (r2 > radius2) {
    // set offdisk pixels to the special value
    imagep[offset] = offdisk;
    // can be helpful to count #times we *needed* to fill with nan
        // if (!isnan(image[offset])) ct_offdisk++; 
      } else if (r2 < radius_inside2) {
    // pass inside pixels thru
    imagep[offset] = image[offset];
      } else {
    // Within annulus: check mode
    if (!mode_adapt) {
      // set, unconditionally, to fill_val
      imagep[offset] = fill_val;
    } else {
      // must fill in: generate new index
      scale = sqrt(radius_inside2/r2); // < 1, pull into annulus interior
      // offset to fill-in value
      dx *= scale;  dy *= scale;
      // round inward toward center
      if (dx > 0) 
        x_off = floor(cenx + dx);
      else
        x_off = ceil(cenx + dx);
      if (dy > 0) 
        y_off = floor(ceny + dy);
      else
        y_off = ceil(ceny + dy);
      // inner offset
      inset = x_off*strx + y_off*stry;
      // set up new value
      imagep[offset] = image[inset];
    } // end if (mode)
    // in both "mode" cases: increment count if the pixel changed
    // (note, this count does not include off-disk pixels)
    if (image[offset] != imagep[offset])
      ct++;
      }
    } // end for (x,y)
  *count = ct;
  /*
  printf("%s: Altered %d pixels\n", progname, ct);
  */
}

/* gateway routine */

#ifdef StaticP  /* undefined under mex */
StaticP
#endif
void
mexFunction(int nlhs, 
        mxArray *plhs[], 
        int nrhs,
        const mxArray *prhs[])
{
  int m, n, maxx, maxy;
  int strx, stry;            // strides in images along x and y dims
  double *center;            // pointers to small vectors of parameters
  int mode_sesw;             // sesw = 1, transposed = 0
  int mode_adapt;            // adaptive edge-setting mode?
  char *mode, *word;         // label for converted mode string
  char errstr[200];
  double datamin, datamax;   // hold range values in output image
  int count;                 // number of pixels changed

  /* 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;
  }
  /* argument checking */
  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);

  /* 
   * `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));
  mode_sesw = -1; 
  mode_adapt = 0; 
  if ((word = extract_mode(mode, &mode_sesw, &mode_adapt)) != 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

  /* center parameters */
  center  = mxGetPr(prhs[ARG_center]);

  /* get size of image */
  m = mxGetM(prhs[ARG_img]);
  n = mxGetN(prhs[ARG_img]);

  /*
   * 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
  }

  /* allocate output image */
  plhs[ARG_out] = mxCreateDoubleMatrix(m, n, mxREAL);
  if (!plhs[ARG_out])
    mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
               "%s: Failed calloc for output image",
               progname), errstr));
  if (nlhs > 1) {
    plhs[ARG_nclean] = mxCreateDoubleMatrix(1, 1, mxREAL);
    if (!plhs[ARG_nclean])
      mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
                 "%s: Failed calloc for output count",
                 progname), errstr));
  }

  // getrange(prhs[ARG_img], &datamin, &datamax);
  // setrange(plhs[ARG_out],  datamin,  datamax);
  setrange(plhs[ARG_out],  0.0,  254.0);   /* set range */

  /* fix the edges */
  if (mxGetScalar(prhs[ARG_delta]) >= 0.0) {
    propagate(&count,
          mxGetPr(prhs[ARG_img]), 
          mxGetPr(plhs[ARG_out]), 
          maxx, maxy,
          strx, stry, 
          center[0]-1,   // xcenter: C origin = (0,0) not (1,1) 
          center[1]-1,   // ycenter: C origin = (0,0) not (1,1) 
          center[2],     // radius
          mxGetScalar(prhs[ARG_delta]),  // width parameter
          mxGetScalar(prhs[ARG_fill ]),  // fill-in value
          mode_adapt
          );
  } else {
    /* delta < 0: set out = img, disregard center */
    memcpy(mxGetPr(plhs[ARG_out]), 
       mxGetPr(prhs[ARG_img]), 
       m * n * sizeof(double));
    count = 0;
  }
  /* plug in count if needed */
  if (nlhs > 1)
    *mxGetPr(plhs[ARG_nclean]) = (double) count;
} 



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

---