JSOC_Documentation: proj/mag/ident/libs/mexfunctions/hmi

00001 #include "mex.h"  /* must appear first */
00002 #include <stdio.h>
00003 #include <math.h>
00004 #include "mexhead.h" /* my mex defines */
00005 #include "Doc/hmi_segment_docstring.h"  /* autogenerated from this file */
00006 
00007 /**************************************************************
00008 
00009 %hmi_segment    driver for HMI segmentation
00010 % 
00011 % [y,s,post,nclean]=hmi_segment(xm,xp,edge,iter,T,beta,alpha,geom,rho,m1,...)
00012 % * Integrated routine for deriving HMI segmentations.  Uses
00013 % models (m1,m2,...), plus a magnetogram-photogram pair (xm,xp),
00014 % to deduce an integer labeling.  Besides images and models,
00015 % it also requires some disk parameters and labeling smoothness 
00016 % parameters.
00017 % * The defaults for iter and T are inherited from mrf_segment_wts,
00018 % and should be looked up there.
00019 % * The posterior probability `post' of the mrf fit can be returned;
00020 % this can quickly verify repeatability.
00021 % * The extreme limb can be a problem for labeling; `edge' helps with
00022 % this.  We set the last edge(1) pixels to the value edge(2), and 
00023 % return the number of such pixels modified, optionally, in nclean.
00024 % Typically edge=[3 1]; both can be real-valued.
00025 % 
00026 % Inputs:
00027 %   real xm(m,n)
00028 %   real xp(m,n)
00029 %   real edge(2)
00030 %   int iter(1) or (2)
00031 %   real T[0] or [1] or [2] or [3] or [4]
00032 %   real beta[1] or [K,K]
00033 %   real alpha[K] or [0] = []
00034 %   real geom(5)
00035 %   real rho
00036 %   real m1(l,k1)
00037 %   ...
00038 %   real mR(l,kR)
00039 % 
00040 % Outputs:
00041 %   int y(m,n)
00042 %   opt real s(R,nS)
00043 %   opt real post
00044 %   opt int nclean
00045 % 
00046 % See Also:  makemrfdiscwts, mrf_segment_wts, mixNprobNd, 
00047 %            clean_edge_label, roi_stats_mag
00048 
00049 % turmon oct 2009, june 2010, june 2011
00050 
00051 ****************************************************************/
00052 
00053 /* standard boilerplate */
00054 #define NARGIN_MIN  11     /* min number of inputs (2 classes) */
00055 #define NARGIN_MAX  13     /* max number of inputs (4 classes) */
00056 #define NARGOUT_MIN 1      /* min number of output args */
00057 #define NARGOUT_MAX 4      /* max number of output args */
00058 
00059 #define ARG_xm    0
00060 #define ARG_xp    1
00061 #define ARG_edge  2 
00062 #define ARG_iter  3 
00063 #define ARG_T     4 
00064 #define ARG_beta  5 
00065 #define ARG_alpha 6 
00066 #define ARG_geom  7 
00067 #define ARG_rho   8 
00068 #define ARG_m1    9 
00069 #define ARG_m2    10
00070 #define ARG_m3    11
00071 #define ARG_m4    12
00072 
00073 #define ARG_y      0
00074 #define ARG_s      1
00075 #define ARG_post   2
00076 #define ARG_nclean 3
00077 
00078 static const char *progname = "hmi_segment";
00079 #define PROGNAME hmi_segment
00080 static const char *in_specs[NARGIN_MAX] = {
00081   "RM",
00082   "RM",
00083   "RV(2)",                // edge
00084   "RS|RV(2)",             // iter
00085   "RS|RV(2)|RV(3)|RV(4)", // T
00086   "RS(1)|RM",             // beta
00087   "RV",                   // alpha
00088   "RV(5)",                // geom
00089   "RS(1)",                // rho
00090   "RM",                   // m1, etc.
00091   "RM",
00092   "RM",
00093   "RM"};
00094 static const char *in_names[NARGIN_MAX] = {
00095   "xm",
00096   "xp",
00097   "edge",
00098   "iter",
00099   "T",
00100   "beta",
00101   "alpha",
00102   "geom",
00103   "rho",
00104   "m1",
00105   "m2",
00106   "m3",
00107   "m4"};
00108 static const char *out_names[NARGOUT_MAX] = {"y", "s", "post", "nclean"};
00109 
00110 // short form of PROGNAME for #include identifiers
00111 #define SHORTNAME Hseg
00112 
00113 // declarations of mex-functions we will use
00114 #include "makemrfdiscwts.h"
00115 #include "mixNprob2d.h"
00116 #include "mrf_segment_wts.h"
00117 #include "clean_edge_label.h"
00118 #include "roi_stats_mag.h"
00119 
00120 // define the argument numbers we actually use
00121 #define MXT_mdw_NARGIN_USE  MXT_mdw_NARGIN_MIN // omit last arg (ell)
00122 #define MXT_mdw_NARGOUT_USE MXT_mdw_NARGOUT_MAX
00123 #define MXT_m2d_NARGIN_USE  MXT_m2d_NARGIN_MIN // omit last arg (mode)
00124 #define MXT_m2d_NARGOUT_USE MXT_m2d_NARGOUT_MAX
00125 #define MXT_msw_NARGIN_USE  MXT_msw_NARGIN_MAX // actually variadic
00126 #define MXT_msw_NARGOUT_USE MXT_msw_NARGOUT_MAX
00127 #define MXT_cel_NARGIN_USE  MXT_cel_NARGIN_MAX
00128 #define MXT_cel_NARGOUT_USE MXT_cel_NARGOUT_MAX
00129 #define MXT_rsm_NARGIN_USE  MXT_rsm_NARGIN_MAX
00130 #define MXT_rsm_NARGOUT_USE 1 // names, combo unused
00131 
00132 
00133 /************************************************************************
00134  *
00135  * Helper routine
00136  *
00137  ************************************************************************/
00138 
00139 /*
00140  * Set off-disk image pixels to nan, in-place
00141  *
00142  * This is currently used for probability maps, but it could
00143  * also be used for the source images.
00144  *
00145  * This assumes standard HMI sesw ordering.
00146  */
00147 
00148 static
00149 void
00150 segment_offdisk_nan(double *p, int maxx, int maxy, double *ctr)
00151 {
00152   const double cenx = ctr[0] - 1.0; // from 1-based to 0-based
00153   const double ceny = ctr[1] - 1.0; // ditto
00154   const double rsun2 = ctr[2]*ctr[2];
00155   const double nan = mxt_getnand();
00156   int x, y, inx;
00157   double dx, dy, r2;
00158 
00159   for (y = 0; y < maxy; y++)
00160     for (x = 0; x < maxx; x++) {
00161       // sene (jpl-mdi): inx = x*maxy + y
00162       // inx = x*maxy + y;
00163       // sesw (native): inx = x*maxy + y
00164       inx = y*maxx + x;
00165       if (isnan(p[inx]))
00166     continue; // already nan
00167       // offsets from center
00168       dx = x - cenx;
00169       dy = y - ceny;
00170       // radius, squared
00171       r2 = dx*dx + dy*dy;
00172       // plunk the nan if off-disk
00173       if (r2 > rsun2)
00174     p[inx] = nan;
00175     } /* end for (x,y) */
00176 }
00177 
00178 /*
00179  * Gateway routine
00180  */
00181 #ifdef StaticP  /* undefined under mex */
00182 StaticP
00183 #endif
00184 void 
00185 mexFunction(int nlhs, 
00186         mxArray *plhs[], 
00187         int nrhs,
00188         const mxArray *prhsC[])
00189 {
00190   mxArray **prhs = (mxArray **) prhsC;    // squelch warnings
00191   // p=mixNprob2d(xm,xp,model)
00192   mxArray *prhs_m2d[MXT_m2d_NARGIN_USE];
00193   mxArray *plhs_m2d[MXT_m2d_NARGOUT_USE];
00194   // dist=makemrfdiscwts(n,del,ctr,rho)
00195   mxArray *prhs_mdw[MXT_mdw_NARGIN_USE];
00196   mxArray *plhs_mdw[MXT_mdw_NARGOUT_USE];
00197   // [yp,post] = mrf_segment_wts(iter,T,beta,alpha,dist,y,p1,p2)
00198   mxArray *prhs_msw[MXT_msw_NARGIN_USE]; // variadic, this is max #args
00199   mxArray *plhs_msw[MXT_msw_NARGOUT_USE];
00200   // res=clean_edge_label(img,center,delta,mode)
00201   mxArray *prhs_cel[MXT_cel_NARGIN_USE];
00202   mxArray *plhs_cel[MXT_cel_NARGOUT_USE];
00203   // roi_stats_mag
00204   mxArray *prhs_rsm[MXT_rsm_NARGIN_USE];
00205   mxArray *plhs_rsm[MXT_rsm_NARGOUT_USE];
00206   // general declarations
00207   int M, N;         // size of disc 
00208   int Nmod;         // number of models
00209   int i;
00210   const int verbose = 0;
00211   char errstr[256];
00212   char *msg;        // returned error message
00213 
00214   /* Hook for introspection (function signature, docstring) */
00215   if (nrhs < 0) { 
00216     plhs[0] = mxt_PackSignature((mxt_Signature) (-nrhs), 
00217                 NARGIN_MIN, NARGIN_MAX, 
00218                 NARGOUT_MIN, NARGOUT_MAX, 
00219                 in_names, in_specs, out_names, docstring);
00220     return;
00221   }
00222   /*
00223    * check args
00224    */
00225   if ((nrhs < NARGIN_MIN) || (nrhs > NARGIN_MAX))
00226     mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
00227                "%s: Expect %d <= input args <= %d",
00228                progname, NARGIN_MIN, NARGIN_MAX), errstr));
00229   if ((nlhs < NARGOUT_MIN) || (nlhs > NARGOUT_MAX))
00230     mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
00231                "%s: Expect %d <= output args <= %d",
00232                progname, NARGOUT_MIN, NARGOUT_MAX), errstr));
00233   if (verbose) printf("in %s\n", progname);
00234   mexargparse(nrhs, prhsC, in_names, in_specs, NULL, progname);
00235 
00236   start_sizechecking();
00237   sizeinit(prhs[ARG_xp]);
00238   sizeagree(prhs[ARG_xm]); /* size(xp) == size(xm) */
00239   sizecheck_msg(progname, in_names, ARG_xp);
00240 
00241   /*
00242    * do the computation
00243    */
00244   M = mxGetM(prhs[ARG_xp]);
00245   N = mxGetN(prhs[ARG_xp]);
00246   Nmod = nrhs - ARG_m1; // number of models
00247 
00248   // disk-weighting function call
00249   //  (ell input argument omitted)
00250   prhs_mdw[MXT_mdw_ARG_n   ] = mxCreateDoubleMatrix(1, 2, mxREAL); // size
00251   prhs_mdw[MXT_mdw_ARG_del ] = mxCreateDoubleScalar(3.0); // del
00252   prhs_mdw[MXT_mdw_ARG_ctr ] = prhs[ARG_geom]; // ctr (5-tuple is OK)
00253   prhs_mdw[MXT_mdw_ARG_rho ] = prhs[ARG_rho]; // rho
00254   prhs_mdw[MXT_mdw_ARG_mode] = mxCreateString("sesw");
00255   if (!prhs_mdw[MXT_mdw_ARG_n  ] ||
00256       !prhs_mdw[MXT_mdw_ARG_del ] ||
00257       !prhs_mdw[MXT_mdw_ARG_mode])
00258     mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
00259                "%s: Failed to allocate makemrfdiscwts inputs",
00260                progname), errstr));
00261   // we know it's non-null, so set it up
00262   mxGetPr(prhs_mdw[MXT_mdw_ARG_n])[0] = (double) M;
00263   mxGetPr(prhs_mdw[MXT_mdw_ARG_n])[1] = (double) N;
00264   if (verbose) printf("segment: calling makemrfdiscwts\n");
00265   msg = main_makemrfdiscwts(MXT_mdw_NARGOUT_USE, plhs_mdw, 
00266                 MXT_mdw_NARGIN_USE,  prhs_mdw);
00267   if (msg)
00268     mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
00269                "%s: Trouble in makemrfdiscwts (%s)",
00270                progname, msg), errstr));
00271   // free what we allocated here, if we no longer need it
00272   mxDestroyArray(prhs_mdw[MXT_mdw_ARG_n   ]);
00273   mxDestroyArray(prhs_mdw[MXT_mdw_ARG_del ]);
00274   mxDestroyArray(prhs_mdw[MXT_mdw_ARG_mode]);
00275 
00276   // set up prob-finder input args
00277   prhs_m2d[MXT_m2d_ARG_i1] = prhs[ARG_xm]; // i1 = mgram
00278   prhs_m2d[MXT_m2d_ARG_i2] = prhs[ARG_xp]; // i2 = pgram
00279   for (i = 0; i < Nmod; i++) {
00280     prhs_m2d[MXT_m2d_ARG_model] = prhs[ARG_m1+i]; // m1, etc.
00281     if (verbose) printf("segment: calling mixNprob2d\n");
00282     msg = main_mixNprob2d(MXT_m2d_NARGOUT_USE, plhs_m2d, 
00283               MXT_m2d_NARGIN_USE,  prhs_m2d);
00284     if (msg)
00285       mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
00286                  "%s: Trouble in mixNprob2d[mod=%d] (%s)",
00287                  progname, i, msg), errstr));
00288     if (verbose) printf("\tprob%d[M/2,N/2] = %g\n", i, 
00289             mxGetPr(plhs_m2d[MXT_m2d_ARG_p])[(M/2)*N+N/2]);
00290     // set off-disk to nan, in-place; only need to do one class
00291     //  (this helper routine assumes sesw ordering)
00292     //  (note that the disk center below is 1-based, not 0-based)
00293     if (i == 0)
00294       segment_offdisk_nan(mxGetPr(plhs_m2d[MXT_m2d_ARG_p]),
00295               M, N, mxGetPr(prhs[ARG_geom]));
00296     // plunk into mrf input arg
00297     prhs_msw[MXT_msw_ARG_lprob1+i] = plhs_m2d[MXT_m2d_ARG_p]; 
00298   }
00299   // free what we allocated here, if we no longer need it
00300   // (we allocated nothing)
00301 
00302   // segmentation
00303   prhs_msw[MXT_msw_ARG_iters] = prhs[ARG_iter];
00304   prhs_msw[MXT_msw_ARG_T]     = prhs[ARG_T];
00305   prhs_msw[MXT_msw_ARG_beta]  = prhs[ARG_beta];
00306   prhs_msw[MXT_msw_ARG_alpha] = prhs[ARG_alpha];
00307   prhs_msw[MXT_msw_ARG_dist]  = plhs_mdw[MXT_mdw_ARG_dist];
00308   prhs_msw[MXT_msw_ARG_y]     = mxCreateDoubleMatrix(0, 0, mxREAL); // IC
00309   if (!prhs_msw[MXT_msw_ARG_y])
00310     mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
00311                "%s: Failed to allocate mrf_segment_wts input",
00312                progname), errstr));
00313   if (verbose) printf("segment: calling mrf_segment\n");
00314   // segment is variadic, so NARGIN uses #models
00315   msg = main_mrf_segment_wts(MXT_msw_NARGOUT_USE,     plhs_msw, 
00316                  MXT_msw_ARG_lprob1+Nmod, prhs_msw);
00317   if (msg)
00318     mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
00319                "%s: Trouble in mrf_segment (%s)",
00320                progname, msg), errstr));
00321   if (verbose) printf("\tseg[M/2,N/2] = %g\n", 
00322               mxGetPr(plhs_msw[MXT_msw_ARG_yp])[(M/2)*N+N/2]);
00323   // free what we allocated here, if we no longer need it
00324   mxDestroyArray(prhs_msw[MXT_msw_ARG_dist]);
00325   mxDestroyArray(prhs_msw[MXT_msw_ARG_y]);
00326   for (i = 0; i < Nmod; i++)
00327     mxDestroyArray(prhs_msw[MXT_msw_ARG_lprob1+i]);
00328 
00329   // set up edge-clean input args
00330   // note, edge = [delta_pix edge_fill_value]
00331   prhs_cel[MXT_cel_ARG_image ] = plhs_msw[MXT_msw_ARG_yp];
00332   prhs_cel[MXT_cel_ARG_center] = prhs[ARG_geom];  // 5-tuple OK
00333   prhs_cel[MXT_cel_ARG_delta ] = mxCreateDoubleScalar(mxGetPr(prhs[ARG_edge])[0]);
00334   prhs_cel[MXT_cel_ARG_fill  ] = mxCreateDoubleScalar(mxGetPr(prhs[ARG_edge])[1]);
00335   prhs_cel[MXT_cel_ARG_mode  ] = mxCreateString("sesw");
00336   if (!prhs_cel[MXT_cel_ARG_delta] ||
00337       !prhs_cel[MXT_cel_ARG_fill ] ||
00338       !prhs_cel[MXT_cel_ARG_mode ])
00339     mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
00340                "%s: Failed to allocate clean_edge_label inputs",
00341                progname), errstr));
00342   if (verbose) printf("segment: calling clean_edge_label\n");
00343   msg = main_clean_edge_label(MXT_cel_NARGOUT_USE, plhs_cel, 
00344                   MXT_cel_NARGIN_USE,  prhs_cel);
00345   if (msg)
00346     mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
00347                "%s: Trouble in clean_edge_label (%s)",
00348                progname, msg), errstr));
00349   if (verbose) printf("\tsegp[M/2,N/2] = %g\n", 
00350               mxGetPr(plhs_cel[MXT_cel_ARG_res])[(M/2)*N+N/2]);
00351   // free what we allocated here, if we no longer need it
00352   mxDestroyArray(prhs_cel[MXT_cel_ARG_image]);
00353   mxDestroyArray(prhs_cel[MXT_cel_ARG_delta]);
00354   mxDestroyArray(prhs_cel[MXT_cel_ARG_fill ]);
00355   mxDestroyArray(prhs_cel[MXT_cel_ARG_mode ]);
00356 
00357   // Optionally find per-class summary statistics
00358   if (nlhs > ARG_s) {
00359     //   s = roi_stats_mag(x,y,mag,geom,nroi,mode)
00360     // (only the first output arg is used)
00361     // (x and y arg are both the activity mask)
00362     // (specify Nmod output roi's in case a class is missing from res)
00363     prhs_rsm[MXT_rsm_ARG_x   ] = plhs_cel[MXT_cel_ARG_res]; // 1..Nr
00364     prhs_rsm[MXT_rsm_ARG_y   ] = plhs_cel[MXT_cel_ARG_res]; // 0/1
00365     prhs_rsm[MXT_rsm_ARG_mag ] = prhs[ARG_xm];   // magnetogram
00366     prhs_rsm[MXT_rsm_ARG_geom] = prhs[ARG_geom]; // require the 5-tuple
00367     prhs_rsm[MXT_rsm_ARG_nroi] = mxCreateDoubleScalar((double) Nmod);
00368     prhs_rsm[MXT_rsm_ARG_mode] = mxCreateString("sesw");
00369     if (!prhs_rsm[MXT_rsm_ARG_mode])
00370       mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
00371                  "%s: Failed to allocate roi_stats_mag inputs",
00372                  progname), errstr));
00373 
00374     if (verbose) printf("segment: calling roi_stats_mag\n");
00375     msg = main_roi_stats_mag(MXT_rsm_NARGOUT_USE, plhs_rsm, 
00376                  MXT_rsm_NARGIN_USE,  prhs_rsm);
00377     if (msg)
00378       mexErrMsgTxt((snprintf(errstr, sizeof(errstr),
00379                  "%s: Trouble in roi_stats_mag (%s)",
00380                  progname, msg), errstr));
00381     // free arrays allocated in this block
00382     mxDestroyArray(prhs_rsm[MXT_rsm_ARG_mode]);
00383     mxDestroyArray(prhs_rsm[MXT_rsm_ARG_nroi]);
00384   }
00385 
00386   /*
00387    *  set up output args
00388    */
00389   if (verbose) printf("segment: setting up outputs...\n");
00390   // y
00391   plhs[ARG_y] = plhs_cel[MXT_cel_ARG_res];
00392   // s
00393   if (nlhs > ARG_s) {
00394     // s was found above
00395     plhs[ARG_s] = plhs_rsm[MXT_rsm_ARG_s]; 
00396   } else {
00397     // s was not found, so no need to free it
00398   }
00399   // post
00400   if (nlhs > ARG_post) {
00401     // asked for posterior output
00402     plhs[ARG_post] = plhs_msw[MXT_msw_ARG_post];
00403   } else {
00404     // did not need posterior (but found it anyway)
00405     mxDestroyArray(plhs_msw[MXT_msw_ARG_post]);
00406   }
00407   // nclean
00408   if (nlhs > ARG_nclean) {
00409     // asked for #pixels cleaned output
00410     plhs[ARG_nclean] = plhs_cel[MXT_cel_ARG_nclean];
00411   } else {
00412     // did not need posterior (but found it anyway)
00413     mxDestroyArray(plhs_cel[MXT_cel_ARG_nclean]);
00414   }
00415 }
00416 
00417 
00418 /* Hook for generic tail matter */
00419 #ifdef MEX2C_TAIL_HOOK
00420 #include "mex2c_tail.h"
00421 #endif
proj/mag/ident/libs/mexfunctions/hmi_segment.c
