[GRASS-SVN] r65425 - in grass-addons/grass7/raster: . r.seg
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jun 9 05:51:48 PDT 2015
Author: avitti
Date: 2015-06-09 05:51:48 -0700 (Tue, 09 Jun 2015)
New Revision: 65425
Added:
grass-addons/grass7/raster/r.seg/
grass-addons/grass7/raster/r.seg/Makefile
grass-addons/grass7/raster/r.seg/main.c
grass-addons/grass7/raster/r.seg/r.seg.html
grass-addons/grass7/raster/r.seg/varseg.c
grass-addons/grass7/raster/r.seg/varseg.h
Log:
Piece-wise smooth approximation and discontinuity detection
Added: grass-addons/grass7/raster/r.seg/Makefile
===================================================================
--- grass-addons/grass7/raster/r.seg/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.seg/Makefile 2015-06-09 12:51:48 UTC (rev 65425)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.seg
+
+LIBES = $(RASTERLIB) $(GISLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/grass7/raster/r.seg/main.c
===================================================================
--- grass-addons/grass7/raster/r.seg/main.c (rev 0)
+++ grass-addons/grass7/raster/r.seg/main.c 2015-06-09 12:51:48 UTC (rev 65425)
@@ -0,0 +1,321 @@
+/****************************************************************************
+ *
+ * MODULE: r.seg
+ *
+ * AUTHOR: Alfonso Vitti <alfonso.vitti unitn.it>
+ *
+ * PURPOSE: generates a piece-wise smooth approximation of the input
+ * raster map and a raster map of the discontinuities (edges) of
+ * the output approximation. The discontinuities of the output
+ * approximation are preserved from being smoothed.
+ *
+ * REFERENCE: http://www.ing.unitn.it/~vittia/phd/vitti_phd.pdf
+ *
+ * COPYRIGHT: (C) 2007-2015
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ *****************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+#include "varseg.h"
+
+
+int main(int argc, char *argv[])
+{
+ char *in_g; /* input, raster map to be segmented */
+ char *out_z; /* output, raster map with detected discontinuities (edges) */
+ char *out_u; /* output, segmented raster map */
+ double lambda; /* scale coefficient */
+ double alpha; /* elasticity coefficient */
+ double kepsilon; /* discontinuities thickness */
+ double beta; /* rigidity coefficient */
+ double tol; /* convergence tolerance */
+ double *mxdf; /* maximum difference betweed two iteration steps */
+ int max_iter; /* max number of numerical iterations */
+ int iter; /* iteration index */
+ const char *mapset; /* current mapset */
+ void *g_row; /* input row buffer */
+ void *out_u_row, *out_z_row;/* output row buffers */
+ int nr, nc, nrc; /* number of rows and colums */
+ int i, j; /* row and column indexes: i=colum=x, j=row=y */
+ int jnc; /* row sequential position, for pointers */
+ int g_fd, out_u_fd, out_z_fd;/* file descriptors */
+ int usek; /* use MSK (MS with the curvature term) */
+ double *g, *u, *z; /* the variables for the actual computation */
+
+ struct Cell_head cellhd; /* GRASS region information */
+ struct History history; /* for map history */
+ struct GModule *module; /* GRASS module for parsing */
+ struct {
+ struct Option *in_g, *out_z, *out_u; /* parameters */
+ } parm;
+ struct {
+ struct Option *lambda, *kepsilon, *alpha, *beta, *tol, *max_iter; /* other parameters */
+ } opts;
+ struct Flag *flag_k; /* flag, k = use MSK instead of MS */
+ RASTER_MAP_TYPE dcell_data_type;/* GRASS raster data type (for DCELL raster) */
+
+
+ /* initialize GRASS environment */
+ G_gisinit(argv[0]);
+
+ /* initialize module */
+ module = G_define_module();
+ G_add_keyword(_("image segmentation"));
+ G_add_keyword(_("edge detection"));
+ G_add_keyword(_("smooth"));
+ module->description =
+ _("Generates a smooth approximation of the input raster and a discontinuity map");
+
+ parm.in_g = G_define_option();
+ parm.in_g->key = "in_g";
+ parm.in_g->type = TYPE_STRING;
+ parm.in_g->required = YES;
+ parm.in_g->description = _("Input raster map to segment");
+ parm.in_g->gisprompt = "old,cell,raster";
+
+ parm.out_u = G_define_option();
+ parm.out_u->key = "out_u";
+ parm.out_u->type = TYPE_STRING;
+ parm.out_u->required = YES;
+ parm.out_u->description = _("Output segmented raster map");
+ parm.out_u->gisprompt = "new,cell,raster";
+
+ parm.out_z = G_define_option();
+ parm.out_z->key = "out_z";
+ parm.out_z->type = TYPE_STRING;
+ parm.out_z->required = YES;
+ parm.out_z->description =
+ _("Output raster map with detected discontinuities");
+ parm.out_z->gisprompt = "new,cell,raster";
+
+ opts.lambda = G_define_option();
+ opts.lambda->key = "lambda";
+ opts.lambda->type = TYPE_DOUBLE;
+ opts.lambda->required = NO;
+ opts.lambda->answer = "1.0";
+ opts.lambda->description = _("Smoothness coefficient [>0]");
+
+ opts.alpha = G_define_option();
+ opts.alpha->key = "alpha";
+ opts.alpha->type = TYPE_DOUBLE;
+ opts.alpha->required = NO;
+ opts.alpha->answer = "1.0";
+ opts.alpha->description = _("Discontinuity coefficient [>0]");
+
+ opts.max_iter = G_define_option();
+ opts.max_iter->key = "mxi";
+ opts.max_iter->type = TYPE_INTEGER;
+ opts.max_iter->required = NO;
+ opts.max_iter->answer = "100";
+ opts.max_iter->description =
+ _("Maximal number of numerical iterations");
+
+ opts.tol = G_define_option();
+ opts.tol->key = "tol";
+ opts.tol->type = TYPE_DOUBLE;
+ opts.tol->required = NO;
+ opts.tol->answer = "0.001";
+ opts.tol->description = _("Convergence tolerance [>0]");
+
+ opts.kepsilon = G_define_option();
+ opts.kepsilon->key = "kepsilon";
+ opts.kepsilon->type = TYPE_DOUBLE;
+ opts.kepsilon->required = NO;
+ opts.kepsilon->answer = "1.0";
+ opts.kepsilon->description =
+ _("Discontinuity thickness [>0]");
+
+ opts.beta = G_define_option();
+ opts.beta->key = "beta";
+ opts.beta->type = TYPE_DOUBLE;
+ opts.beta->required = NO;
+ opts.beta->answer = "0.0";
+ opts.beta->description = _("Curvature coefficient [>=0]");
+ /*
+ * beta = 0 leads to MS
+ * beta > 0 leads to MSK
+ * Due to a different implementation of MSK wrt MS,
+ * the values of the parameters lambda and alpha in MSK
+ * have to be set independently from the values used in MS
+ */
+
+ flag_k = G_define_flag();
+ flag_k->key = 'k';
+ flag_k->description = _("Activate MSK model (Mumford-Shah with curvature term)");
+
+
+ /* parameters and flags parser */
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /* store parameters and flags to variables */
+ in_g = parm.in_g->answer;
+ out_u = parm.out_u->answer;
+ out_z = parm.out_z->answer;
+
+ lambda = atof(opts.lambda->answer);
+ kepsilon = atof(opts.kepsilon->answer);
+ alpha = atof(opts.alpha->answer);
+ beta = atof(opts.beta->answer);
+ tol = atof(opts.tol->answer);
+ max_iter = atoi(opts.max_iter->answer);
+
+ if (((usek = (flag_k->answer) == 0) && (beta != 0.0)))
+ G_warning(_
+ ("Beta is not zero and you have not activated the MSK formulation: \n \
+ beta will be ignored and MS (default) will be used."));
+
+ if (((usek = (flag_k->answer) == 1) && (beta == 0.0)))
+ G_warning(_
+ ("You have activated the MSK formulation, but beta is zero:\n \
+ beta should be greater than zero in MSK."));
+
+ /* check existence and names of raster maps */
+ mapset = G_find_raster2(in_g, "");
+ if (!mapset)
+ G_fatal_error(_("raster file [%s] not found"), in_g);
+ if (G_legal_filename(out_u) < 0)
+ G_fatal_error(_("[%s] is an illegal file name"), out_u);
+ G_check_input_output_name(in_g, out_u, G_FATAL_EXIT);
+ if (G_legal_filename(out_z) < 0)
+ G_fatal_error(_("[%s] is an illegal file name"), out_z);
+ G_check_input_output_name(in_g, out_z, G_FATAL_EXIT);
+ if (strcmp(out_u, out_z) == 0)
+ G_fatal_error(_("Output raster maps have the same name [%s]"), out_u);
+
+
+ /* -------------------------------------------------------------------- */
+ /* Do the work */
+ /* -------------------------------------------------------------------- */
+
+ /* data type for the computation */
+ dcell_data_type = DCELL_TYPE;
+
+
+ /* get the window dimention */
+ nr = Rast_window_rows();
+ nc = Rast_window_cols();
+ nrc = nr * nc;
+
+
+ /* allocate the memory for the varialbes used in the computation */
+ g = (DCELL *) G_malloc(sizeof(DCELL) * nrc);
+ u = (DCELL *) G_malloc(sizeof(DCELL) * nrc);
+ z = (DCELL *) G_malloc(sizeof(DCELL) * nrc);
+
+ mxdf = (double *)malloc(sizeof(double));
+
+
+ /* open the input raster map for reading */
+ if ((g_fd = Rast_open_old(in_g, mapset)) < 0)
+ G_fatal_error(_("cannot open raster file [%s]"), in_g);
+ Rast_get_cellhd(in_g, mapset, &cellhd);
+
+
+ /* allocate the buffer for storing the values of the input raster map */
+ g_row = Rast_allocate_buf(dcell_data_type);
+
+
+ /* read the input raster map + fill up the variable pointers with pixel values */
+ for (j = 0; j < nr; j++) {
+ jnc = j * nc;
+
+ Rast_get_row(g_fd, g_row, j, dcell_data_type);
+
+ for (i = 0; i < nc; i++) {
+ *(g + jnc + i) = ((DCELL *) g_row)[i];
+ *(u + jnc + i) = ((DCELL *) g_row)[i];
+ *(z + jnc + i) = 1.0;
+ }
+ }
+
+
+ /* close the input raster map and free memory */
+ Rast_close(g_fd);
+ G_free(g_row);
+
+
+ /* the first iteration is always performed */
+ iter = 1;
+
+
+ /* call the library function to perform the segmentation */
+ if (usek == 0) {
+ *mxdf = 0;
+ ms_n(g, u, z, lambda, kepsilon, alpha, mxdf, nr, nc);
+ while ((*mxdf > tol) && (iter <= max_iter)) {
+ *mxdf = 0;
+ ms_n(g, u, z, lambda, kepsilon, alpha, mxdf, nr, nc);
+ iter += 1;
+ }
+ }
+ else {
+ *mxdf = 0;
+ msk_n(g, u, z, lambda, kepsilon, alpha, beta, mxdf, nr, nc);
+ while ((*mxdf > tol) && (iter <= max_iter)) {
+ *mxdf = 0;
+ msk_n(g, u, z, lambda, kepsilon, alpha, beta, mxdf, nr, nc);
+ iter += 1;
+ }
+ }
+
+
+ /* print the total number of iteration performed */
+ G_message("\nr.seg iterations: %i\n", iter);
+ G_message("\n");
+
+ /* open the output raster maps for writing */
+ if ((out_u_fd = Rast_open_new(out_u, dcell_data_type)) < 0)
+ G_fatal_error(_("cannot open raster file [%s]"), out_u);
+ if ((out_z_fd = Rast_open_new(out_z, dcell_data_type)) < 0)
+ G_fatal_error(_("cannot open raster file [%s]"), out_z);
+
+
+ /* allocate the buffer for storing the values of the output raster maps */
+ out_u_row = Rast_allocate_buf(dcell_data_type);
+ out_z_row = Rast_allocate_buf(dcell_data_type);
+
+
+ /* fill up the output buffers with result values + write the output raster maps */
+ for (j = 0; j < nr; j++) {
+ jnc = j * nc;
+ for (i = 0; i < nc; i++) {
+ ((DCELL *) out_u_row)[i] = *(u + jnc + i);
+ ((DCELL *) out_z_row)[i] = *(z + jnc + i);
+ }
+ Rast_put_row(out_u_fd, out_u_row, dcell_data_type);
+ Rast_put_row(out_z_fd, out_z_row, dcell_data_type);
+ }
+
+
+ /* close the input raster map and free memory */
+ Rast_close(out_u_fd);
+ G_free(out_u_row);
+ Rast_close(out_z_fd);
+ G_free(out_z_row);
+ G_free(g);
+ G_free(u);
+ G_free(z);
+
+
+ /* write history file */
+ Rast_short_history(out_u, "raster", &history);
+ Rast_command_history(&history);
+ Rast_append_format_history(&history, "\nIterations = %i", iter);
+ Rast_write_history(out_u, &history);
+ Rast_write_history(out_z, &history);
+
+
+ /* exit */
+ exit(EXIT_SUCCESS);
+}
+
Added: grass-addons/grass7/raster/r.seg/r.seg.html
===================================================================
--- grass-addons/grass7/raster/r.seg/r.seg.html (rev 0)
+++ grass-addons/grass7/raster/r.seg/r.seg.html 2015-06-09 12:51:48 UTC (rev 65425)
@@ -0,0 +1,161 @@
+<h2>DESCRIPTION</h2>
+
+<em><b>r.seg</b></em> generates a piece-wise smooth approximation of the
+input raster map and a raster map of the discontinuities of the output
+approximation. <br>
+
+The discontinuities of the output approximation are preserved from being
+smoothed. The values of the discontinuity map are close to one where the
+output approximation is "homogeneous", where the output approximation
+has discontinuities (edges) the values are close to zero. <br>
+
+The module makes use of the <em>varseg</em> library which implements
+the Mumford-Shah [1] variational model for image segmentation. The
+Mumford-Shah variational model with curvature term [2] is also implemented
+in the library. The curvature term prevents the discontinuities from being
+shortened too much when the parameter alpha is set to very high values,
+(this happens very rarely). <br>
+
+Some examples of use of the module can be found <a
+href="http://www.ing.unitn.it/~vittia/sw/sw_index.html">here</a> and in <a
+href="http://download.OSGeo.org/OSGeo/foss4g/2009/SPREP/2Thu/Parkside%20GO4/1500/Thu
+G04 1545 Zatelli.pdf">this presentation [FOSS4G 2009 - pdf]</a>. <br>
+For details on the numerical implementation see [3].
+
+<h2>NOTES</h2>
+Remove any MASK before the execution of the module. If
+a MASK is present, the module stops after just one iteration.
+<br><br>
+
+Replace (<em>r.null</em>) any null data with the map average value (get with <em>r.univar</em>).
+<br><br>
+
+The segmentation depends on the parameters alpha and lambda:
+<ul>
+<li> alpha controls how many discontinuities are allowed to exist.
+<li> lambda controls the smoothness of the solution.
+<li> It is not possible to select the values of the parameters in an
+ automatic way. Test some different values to understand their
+ influence on the results. Try the following procedure:
+<ul>
+ <li> run the module with both alpha and lambda set to 1.0
+ <li> run the module with alpha set to 1.0 and different values for lambda
+ <br> e.g., 0.01, 0.1, 1, 10, 100
+ <li> run the module with lambda set to 1.0 and different values for alpha
+ <br> e.g., let's say 0.01, 0.1, 1, 10, 100
+ <li> see how the segmentations change and select the values that
+ produce the result that best fits your requirements.
+</ul>
+</ul>
+
+The module computes the segmentation by means of an iterative
+procedure.<br>
+The module stops either when the number of iterations
+reaches the maximum number of iterations [mxi] or when the maximum
+difference between the solutions of two successive iterations is less than
+the convergence tolerance [tol].<br>
+To stop the iteration procedure,
+it is easier to act on the maximum number of iterations parameter [mxi]
+than on the convergence tolerance parameter [tol].<br>
+The number of
+iterations needed to reach the convergence tolerance increases for high
+values of the parameter lambda. The larger the total number of pixels
+of the input raster map the larger the number of iterations will be.
+<br><br>
+
+The data type of the output raster maps is DOUBLE PRECISION. <br><br>
+
+The module works on one raster map at a time, imagery groups are not
+supported. <br><br>
+
+To avoid to inappropriately re-sampled the input raster map, the settings
+for the current region should be set so that:
+<ul>
+<li> the resolution of the region matches the resolution of the
+ input raster map;
+<li>the boundaries of the region are lined up along the edges of the nearest
+cells in the input raster map.
+</ul>
+
+The discontinuity thickness should be changed for test purposes only.
+<br><br>
+
+The actual need to use the MSK model should be very rare, see [3].
+Due to a different implementation of the MSK model with respect to MS
+one, the values of the parameters lambda and alpha in MSK have to be
+set independently from the values used in MS.
+
+<h2>EXAMPLE</h2>
+
+This example is based the <a
+href="http://grass.OSGeo.org/download/data.php">North Carolina GRASS sample
+data set</a>, [complete GRASS location].
+
+
+<div><pre class="code">
+# set the region to match the <em>ortho_2001_t792_1m</em> raster map:
+g.region rast=ortho_2001_t792_1m
+
+# select a smaller region:
+g.region n=221725 s=220225 w=638350 e=639550
+
+# run r.seg:
+r.seg in_g=ortho_2001_t792_1m at PERMANENT out_u=u_OF out_z=z_OF lambda=10 alpha=200 mxi=250
+
+# for a better visualization of the output raster map <em>u_OF</em>, set its color table to:
+r.colors u_OF rast=ortho_2001_t792_1m
+
+# compute the difference between the input raster map and the output raster map <em>u_OF</em>:
+r.mapcalc "diff = abs(ortho_2001_t792_1m at PERMANENT - u_OF)"
+
+# for a better visualization of the differences, compute the natural logarithm of the <em>diff</em> map:
+r.mapcalc "log_diff = log(1 + diff)"
+
+# and set its color table to the "differences" style:
+r.colors log_diff color=differences
+
+# for a better visualization of the output raster map <em>u_OF</em>, set its color table to:
+r.colors z_OF color=bgyr
+
+# run r.seg with different parameter values:
+r.seg in_g=ortho_2001_t792_1m at PERMANENT out_u=u1_OF out_z=z1_OF lambda=10 alpha=65 mxi=250
+r.seg in_g=ortho_2001_t792_1m at PERMANENT out_u=u2_OF out_z=z2_OF lambda=10 alpha=600 mxi=250
+r.seg in_g=ortho_2001_t792_1m at PERMANENT out_u=u3_OF out_z=z3_OF lambda=0.1 alpha=200 mxi=250
+r.seg in_g=ortho_2001_t792_1m at PERMANENT out_u=u4_OF out_z=z4_OF lambda=1 alpha=200 mxi=250
+
+# visualize and compare the different results
+</pre></div>
+
+<h2>REFERENCE</h2>
+
+<ul> <li> <b>[1]</b> D. Mumford and J. Shah. <em>Optimal Approximation by
+Piecewise Smooth Functions and Associated Variational Problems</em>. <br>
+Communications on Pure Applied Mathematics, 42:577-685, 1989.
+
+<li> <b>[2]</b> R. March and M. Dozio. <em>A variational method for the
+recovery of smooth boundaries</em>. <br> Image and Vision Computing,
+15:705-712, 1997.
+
+<li> <b>[3]</b> A. Vitti. <em>Free discontinuity
+problems in image and signal segmentatiion</em>. <br>
+Ph.D. Thesis - University of Trento (Italy), 2008. <br> <a
+href="http://www.ing.unitn.it/~vittia/misc/vitti_phd.pdf">http://www.ing.unitn.it/~vittia/misc/vitti_phd.pdf</a>
+</ul>
+
+
+<h2>SEE ALSO</h2>
+
+<em><a href="i.smap.html">i.smap</a></em>,
+<em><a href="i.zc.html">i.zc</a></em>,
+<em><a href="r.mfilter.html">r.mfilter</a></em>
+
+
+<h2>AUTHOR</h2>
+
+Alfonso Vitti <br>
+ Dept. Civil and
+Environmental Engineering <br>
+ University of Trento - Italy<br>
+ alfonso.vitti [at] ing.unitn.it
+
+<p><i>Last changed: $Date: 2010-08-10 12:00:00 +0200 (Tue, 10 Aug 2010)$</i>
Added: grass-addons/grass7/raster/r.seg/varseg.c
===================================================================
--- grass-addons/grass7/raster/r.seg/varseg.c (rev 0)
+++ grass-addons/grass7/raster/r.seg/varseg.c 2015-06-09 12:51:48 UTC (rev 65425)
@@ -0,0 +1,1123 @@
+/****************************************************************************
+ *
+ * LIBRARY: varseg
+ *
+ * AUTHOR: Alfonso Vitti <alfonso.vitti ing.unitn.it>
+ *
+ * PURPOSE: Implements the Mumford-Shah (MS) variational model for image
+ * segmentation.
+ * The Mumford-Shah variational model with cirvature term (MSK)
+ * is also implemented.
+ *
+ * The models generate a piece-wise smooth approximation of the
+ * input image and an image of the discontinuities (edges)
+ * of the output approximation.
+ * The discontinuities of the output approximation are preserved
+ * from being smoothed.
+ *
+ * REFERENCE: http://www.ing.unitn.it/~vittia/phd/vitti_phd.pdf
+ *
+ * COPYRIGHT: (C) 2007-2010
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2).
+ *
+ *****************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "varseg.h"
+
+#define max(a,b) ((a)>(b) ? (a) : (b))
+
+
+/* ========================================================================== */
+/*!
+ * \brief MS_N --- MUMFORD-SHAH (Gauss-Seidel method)
+ *
+ * Implements the Mumford-Shah variational model for image segmentation.
+ *
+ * A smooth approximation of the input image is produced.
+ * The model explicitely detects the discontinuities of the output approximation
+ * while preserving the discontinuities from being smoothed.
+ *
+ * A discontinuity image is also produced.
+ * The values of the discontinuity image are close to 1 where the output approximation
+ * is "homogeneous", where the output approximation has discontinuities (edges) the
+ * values are close to 0.
+ *
+ *
+ * \param[in] *g (double; pointer) - The input image
+ *
+ * \param[in] *u (double; pointer) - The output approximation
+ *
+ * \param[in] *z (double; pointer) - The output discontinuity map
+ *
+ * \param[in] lambda (double) - The scale factor parameter [positive]
+ *
+ * \param[in] kepsilon (double) - The discontinuities thickness [positive]
+ *
+ * \param[in] alpha (double) - The elasticity parameter [positive]
+ *
+ * \param[in] *mxdf (double; pointer) - The max difference in the approximation between two iteration steps
+ *
+ * \param[in] nr (int) - The number of rows in the input image
+ *
+ * \param[in] nc (int) - The number of colums in the input image
+ *
+ * \return void [shold be changed to int for error handling]
+ */
+
+/** MS_N --- MUMFORD-SHAH (Gauss-Seidel method) **/
+void ms_n(double *g, double *u, double *z, double lambda, double kepsilon,
+ double alpha, double *mxdf, int nr, int nc)
+{
+ int i, j, jnc, nrc;
+ double old_u; /* ms_n */
+ double u_hat, z_star, z_tilde, ux, uy; /* ms */
+ double num, den;
+ double k2, k, h2, h = 1;
+
+ nrc = nr * nc;
+ h2 = h * h;
+ k = kepsilon;
+ k2 = kepsilon * kepsilon;
+
+ /* Boundary conditions (Neumann) */
+
+ /*** image corners ***/
+ *u = *(u + nc + 1); /* ul corner */
+ *(u + nc - 1) = *(u + 2 * nc - 2); /* ur corner */
+ *(u + nrc - nc) = *(u + nrc - 2 * nc + 1); /* ll corner */
+ *(u + nrc - 1) = *(u + nrc - nc - 2); /* lr corner */
+ *z = *(z + nc + 1); /* ul corner */
+ *(z + nc - 1) = *(z + 2 * nc - 2); /* ur corner */
+ *(z + nrc - nc) = *(z + nrc - 2 * nc + 1); /* ll corner */
+ *(z + nrc - 1) = *(z + nrc - nc - 2); /* lr corner */
+
+ /* this is an alternative for the corners using the mean value of the x and y closest points
+ * *u = 0.5 * ( *(u+1) + *(u+nc) );
+ * *(u+nc-1) = 0.5 * ( *(u+nc-2) + *(u+2*nc-1) );
+ * *(u+nrc-nc) = 0.5 * ( *(u+nrc-nc+1) + *(u+nrc-2*nc) );
+ * *(u+nrc-1) = 0.5 * ( *(u+nrc-2) +*(u+nrc-nc-1) );
+ * *z = 0.5 * ( *(z+1) + *(z+nc) );
+ * *(z+nc-1) = 0.5 * ( *(z+nc-2) + *(z+2*nc-1) );
+ * *(z+nrc-nc) = 0.5 * ( *(z+nrc-nc+1) + *(z+nrc-2*nc) );
+ * *(z+nrc-1) = 0.5 * ( *(z+nrc-2) +*(z+nrc-nc-1) );
+ */
+
+ /*** image edges ***/
+ for (i = 1; i < nc - 1; i++) { /* top and bottom edges */
+ *(u + i) = *(u + nc + i); /* t edge */
+ *(u + nrc - nc + i) = *(u + nrc - 2 * nc + i); /* b edge */
+ *(z + i) = *(z + nc + i); /* t edge */
+ *(z + nrc - nc + i) = *(z + nrc - 2 * nc + i); /* b edge */
+ }
+ for (j = 1; j < nr - 1; j++) { /* left and right edges */
+ jnc = j * nc;
+ *(u + jnc) = *(u + jnc + 1); /* r edge */
+ *(u + jnc + nc - 1) = *(u + jnc + nc - 2); /* l edge */
+ *(z + jnc) = *(z + jnc + 1); /* r edge */
+ *(z + jnc + nc - 1) = *(z + jnc + nc - 2); /* l edge */
+ }
+
+ for (j = 1; j < nr - 1; j++) {
+ jnc = j * nc;
+ for (i = 1; i < nc - 1; i++) {
+ old_u = *(u + jnc + i);
+ /*
+ * x_i+1,j --> *(x+jnc+i+1)
+ * x_i-1,j --> *(x+jnc+i-1)
+ * x_i,j+1 --> *(x+jnc+nc+i)
+ * x_i,j-1 --> *(x+jnc-nc+i)
+ */
+ u_hat =
+ *(z + jnc + i + 1) * *(z + jnc + i + 1) * *(u + jnc + i + 1) +
+ *(z + jnc + i - 1) * *(z + jnc + i - 1) * *(u + jnc + i - 1) +
+ *(z + jnc + nc + i) * *(z + jnc + nc + i) * *(u + jnc + nc +
+ i) + *(z + jnc -
+ nc +
+ i) *
+ *(z + jnc - nc + i) * *(u + jnc - nc + i);
+ z_star =
+ *(z + jnc + i + 1) * *(z + jnc + i + 1) + *(z + jnc + i -
+ 1) * *(z + jnc +
+ i - 1) +
+ *(z + jnc + nc + i) * *(z + jnc + nc + i) + *(z + jnc - nc +
+ i) * *(z + jnc -
+ nc + i);
+
+ num = lambda * u_hat + h2 * *(g + jnc + i);
+ den = lambda * z_star + h2;
+
+ *(u + jnc + i) = num / den;
+
+ z_tilde =
+ *(z + jnc + i + 1) + *(z + jnc + i - 1) + *(z + jnc + nc +
+ i) + *(z + jnc -
+ nc + i);
+ ux = *(u + jnc + i + 1) - *(u + jnc + i - 1);
+ uy = *(u + jnc + nc + i) - *(u + jnc - nc + i);
+
+ num = 4 * z_tilde + h2 * k2;
+ den = 16 + h2 * k2 + lambda * k * (ux * ux + uy * uy) / alpha;
+
+ *(z + jnc + i) = num / den;
+
+ *mxdf = max(*mxdf, fabs(old_u - *(u + jnc + i)));
+ }
+ }
+ return;
+}
+
+/* ========================================================================== */
+/* ========================================================================== */
+
+
+/* ========================================================================== */
+/*!
+ * \brief MSK_N --- MUMFORD-SHAH with CURVATURE term (Gauss-Seidel method)
+ *
+ * Implements the Mumford-Shah variational model with curvature term for image segmentation.
+ *
+ * A smooth approximation of the input image is produced.
+ * The model explicitely detects the discontinuities of the output approximation
+ * while preserving the discontinuities from being smoothed.
+ *
+ * A discontinuity image is also produced.
+ * The values of the discontinuity image are close to 1 where the output approximation
+ * is "homogeneous", where the output approximation has discontinuities (edges) the
+ * values are close to 0.
+ *
+ * The curvature term prevents the discontinuities from being shortened too much when the parameter
+ * alpha is set to very high values.
+ *
+ *
+ * \param[in] *g (double; pointer) - The input image
+ *
+ * \param[in] *u (double; pointer) - The output approximation
+ *
+ * \param[in] *z (double; pointer) - The output discontinuity map
+ *
+ * \param[in] lambda (double) - The scale factor parameter [positive]
+ *
+ * \param[in] kepsilon (double) - The discontinuities thickness [positive]
+ *
+ * \param[in] alpha (double) - The elasticity parameter [positive]
+ *
+ * \param[in] beta (double) - The rigidity parameter [positive]
+ *
+ * \param[in] *mxdf (double; pointer) - The max difference in the approximation between two iteration steps
+ *
+ * \param[in] nr (int) - The number of rows in the input image
+ *
+ * \param[in] nc (int) - The number of colums in the input image
+ *
+ * \return void [shold be changed to int for error handling]
+ */
+
+/** MSK_N --- MUMFORD-SHAH WITH CURVATURE (Gauss-Seidel method) **/
+void msk_n(double *g, double *u, double *z, double lambda, double kepsilon,
+ double alpha, double beta, double *mxdf, int nr, int nc)
+{
+ int i, j, jnc, nrc;
+ double old_u; /* msk_n */
+ double u_tilde, ux, uy, zx, zy, z_tilde, z_hat, z_check; /* msk */
+ double a, b, c, d, e; /* msk */
+ double num, den;
+ double k2, k, h2, h = 1;
+ double ta;
+
+ nrc = nr * nc;
+ h2 = h * h;
+ k = kepsilon;
+ k2 = kepsilon * kepsilon;
+
+ /* Remove first image boundary using BC values */
+ *u = *(u + nc + 1); /* ul corner */
+ *(u + nc - 1) = *(u + 2 * nc - 2); /* ur corner */
+ *(u + nrc - nc) = *(u + nrc - 2 * nc + 1); /* ll corner */
+ *(u + nrc - 1) = *(u + nrc - nc - 2); /* lr corner */
+ *z = *(z + nc + 1); /* ul corner */
+ *(z + nc - 1) = *(z + 2 * nc - 2); /* ur corner */
+ *(z + nrc - nc) = *(z + nrc - 2 * nc + 1); /* ll corner */
+ *(z + nrc - 1) = *(z + nrc - nc - 2); /* lr corner */
+ for (i = 1; i < nc - 1; i++) { /* top and bottom edges */
+ *(u + i) = *(u + nc + i); /* t edge */
+ *(u + nrc - nc + i) = *(u + nrc - 2 * nc + i); /* b edge */
+ *(z + i) = *(z + nc + i); /* t edge */
+ *(z + nrc - nc + i) = *(z + nrc - 2 * nc + i); /* b edge */
+ }
+ for (j = 1; j < nr - 1; j++) { /* left and right edges */
+ jnc = j * nc;
+ *(u + jnc) = *(u + jnc + 1); /* r edge */
+ *(u + jnc + nc - 1) = *(u + jnc + nc - 2); /* l edge */
+ *(z + jnc) = *(z + jnc + 1); /* r edge */
+ *(z + jnc + nc - 1) = *(z + jnc + nc - 2); /* l edge */
+ }
+
+ /* Boundary conditions (Neumann) */
+
+ /*** image corners ***/
+ *(u + nc + 1) = *(u + 2 * nc + 2); /* ul corner */
+ *(u + 2 * nc - 2) = *(u + 3 * nc - 3); /* ur corner */
+ *(u + nrc - 2 * nc + 1) = *(u + nrc - 3 * nc + 2); /* ll corner */
+ *(u + nrc - nc - 2) = *(u + nrc - 2 * nc - 3); /* lr corner */
+ *(z + nc + 1) = *(z + 2 * nc + 2); /* ul corner */
+ *(z + 2 * nc - 2) = *(z + 3 * nc - 3); /* ur corner */
+ *(z + nrc - 2 * nc + 1) = *(z + nrc - 3 * nc + 2); /* ll corner */
+ *(z + nrc - nc - 2) = *(z + nrc - 2 * nc - 3); /* lr corner */
+
+ /* this is an alternative for the corners using the mean value of the x and y closest points
+ * *(u+nc+1) = 0.5 * ( *(u+nc+1) + *(u+2*nc+1) );
+ * *(u+2*nc-2) = 0.5 * ( *(u+2*nc-3) + *(u+3*nc-2) );
+ * *(u+nrc-2*nc+1) = 0.5 * ( *(u+nrc-2*nc+2) + *(u+nrc-3*nc+1) );
+ * *(u+nrc-nc-2) = 0.5 * ( *(u+nrc-nc-3) +*(u+nrc-2*nc-2) );
+ * *z = 0.5 * ( *(z+nc+1) + *(z+2*nc+1) );
+ * *(z+nc-1) = 0.5 * ( *(z+2*nc-3) + *(z+3*nc-2) );
+ * *(z+nrc-nc) = 0.5 * ( *(z+nrc-2*nc+2) + *(z+nrc-3*nc+1) );
+ * *(z+nrc-1) = 0.5 * ( *(z+nrc-nc-3) +*(z+nrc-2*nc-2) );
+ */
+
+ /*** image edges ***/
+ for (i = 2; i < nc - 2; i++) { /* top and bottom edges */
+ *(u + nc + i) = *(u + 2 * nc + i); /* t edge */
+ *(u + nrc - 2 * nc + i) = *(u + nrc - 3 * nc + i); /* b edge */
+ *(z + nc + i) = *(z + 2 * nc + i); /* t edge */
+ *(z + nrc - 2 * nc + i) = *(z + nrc - 3 * nc + i); /* b edge */
+ }
+ for (j = 2; j < nr - 2; j++) { /* left and right edges */
+ jnc = j * nc;
+ *(u + jnc + 1) = *(u + jnc + 2); /* r edge */
+ *(u + jnc + nc - 2) = *(u + jnc + nc - 3); /* l edge */
+ *(z + jnc + 1) = *(z + jnc + 2); /* r edge */
+ *(z + jnc + nc - 2) = *(z + jnc + nc - 3); /* l edge */
+ }
+
+ for (j = 2; j < nr - 2; j++) {
+ jnc = j * nc;
+ for (i = 2; i < nc - 2; i++) {
+ old_u = *(u + jnc + i);
+
+ u_tilde =
+ *(u + jnc + i + 1) + *(u + jnc + i - 1) + *(u + jnc + nc +
+ i) + *(u + jnc -
+ nc + i);
+ ux = *(u + jnc + i + 1) - *(u + jnc + i - 1);
+ uy = *(u + jnc + nc + i) - *(u + jnc - nc + i);
+ zx = *(z + jnc + i + 1) - *(z + jnc + i - 1);
+ zy = *(z + jnc + nc + i) - *(z + jnc - nc + i);
+ num =
+ 0.5 * lambda * *(z + jnc + i) * (zx * ux + zy * uy) +
+ lambda * *(z + jnc + i) * *(z + jnc + i) * u_tilde +
+ h2 * *(g + jnc + i);
+ den = 4 * lambda * *(z + jnc + i) * *(z + jnc + i) + h2;
+
+ *(u + jnc + i) = num / den;
+
+ z_tilde =
+ *(z + jnc + i + 1) + *(z + jnc + i - 1) + *(z + jnc + nc +
+ i) + *(z + jnc -
+ nc + i);
+ z_hat =
+ *(z + jnc + nc + i + 1) + *(z + jnc - nc + i + 1) + *(z +
+ jnc +
+ nc + i -
+ 1) +
+ *(z + jnc - nc + i - 1);
+ z_check =
+ *(z + jnc + i + 2) + *(z + jnc + i - 2) + *(z + jnc + 2 * nc +
+ i) + *(z + jnc -
+ 2 * nc +
+ i);
+ a = alpha * z_tilde + 4 * alpha * h2 * k2 * *(z + jnc + i) * *(z +
+ jnc
+ +
+ i)
+ * *(z + jnc + i);
+ b = 4 * beta * (8 * z_tilde - 2 * z_hat - z_check) / h2;
+ c = -16 * beta * k2 * (3 * *(z + jnc + i) * *(z + jnc + i) +
+ 1) * (z_tilde - 4 * *(z + jnc + i));
+ d = 64 * beta * k2 * *(z + jnc +
+ i) * (3 * *(z + jnc + i) * *(z + jnc + i) -
+ 1);
+ e = 64 * beta * h2 * k2 * k2 * *(z + jnc + i) * *(z + jnc +
+ i) * *(z + jnc +
+ i) * (3 *
+ *(z
+ +
+ jnc
+ +
+ i)
+ *
+ *(z
+ +
+ jnc
+ +
+ i)
+ -
+ 2);
+ num = a + b + c + d + e;
+ ta = a;
+
+ a = 4 * alpha +
+ 2 * alpha * h2 * k2 * (3 * *(z + jnc + i) * *(z + jnc + i) -
+ 1) + 0.25 * lambda * k * (ux * ux +
+ uy * uy);
+ b = 16 * beta * k2 * (h2 * k2 + 5 / (h2 * k2) - 4);
+ c = -12 * beta * k2 * (zx * zx + zy * zy);
+ d = -96 * beta * k2 * *(z + jnc + i) * (z_tilde -
+ 4 * *(z + jnc + i));
+ e = 192 * beta * k2 * *(z + jnc + i) * (1 - h2 * k2) +
+ 240 * beta * h2 * k2 * k2 * *(z + jnc + i) * *(z + jnc +
+ i) * *(z +
+ jnc +
+ i) *
+ *(z + jnc + i);
+ den = a + b + c + d + e;
+
+ *(z + jnc + i) = num / den;
+
+ *mxdf = max(*mxdf, fabs(old_u - *(u + jnc + i)));
+ }
+ }
+ return;
+}
+
+/* ========================================================================== */
+/* ========================================================================== */
+
+
+/* ========================================================================== */
+/*!
+ * \brief MS_O --- MUMFORD-SHAH (Jacobi method)
+ *
+ * Implements the Mumford-Shah variational model for image segmentation.
+ *
+ * A smooth approximation of the input image is produced.
+ * The model explicitely detects the discontinuities of the output approximation
+ * while preserving the discontinuities from being smoothed.
+ *
+ * A discontinuity image is also produced.
+ * The values of the discontinuity image are close to 1 where the output approximation
+ * is "homogeneous", where the output approximation has discontinuities (edges) the
+ * values are close to 0.
+ *
+ *
+ * \param[in] *g (double; pointer) - The input image
+ *
+ * \param[in] *u (double; pointer) - The output approximation
+ *
+ * \param[in] *z (double; pointer) - The output discontinuity map
+ *
+ * \param[in] lambda (double) - The scale factor parameter [positive]
+ *
+ * \param[in] kepsilon (double) - The discontinuities thickness [positive]
+ *
+ * \param[in] alpha (double) - The elasticity parameter [positive]
+ *
+ * \param[in] *mxdf (double; pointer) - The max difference in the approximation between two iteration steps
+ *
+ * \param[in] nr (int) - The number of rows in the input image
+ *
+ * \param[in] nc (int) - The number of colums in the input image
+ *
+ * \return void [shold be changed to int for error handling]
+ */
+
+/** MS_O --- MUMFORD-SHAH (Jacobi method) **/
+void ms_o(double *g, double *u, double *z, double lambda, double kepsilon,
+ double alpha, double *mxdf, int nr, int nc)
+{
+ int i, j, jnc, nrc;
+ double *old_u, *old_z; /* ms_o */
+ double u_hat, z_star, z_tilde, ux, uy; /* ms */
+ double num, den;
+ double k2, k, h2, h = 1;
+
+ nrc = nr * nc;
+ h2 = h * h;
+ k = kepsilon;
+ k2 = kepsilon * kepsilon;
+
+ /* Boundary conditions (Neumann) */
+
+ /*** image corners ***/
+ *u = *(u + nc + 1); /* ul corner */
+ *(u + nc - 1) = *(u + 2 * nc - 2); /* ur corner */
+ *(u + nrc - nc) = *(u + nrc - 2 * nc + 1); /* ll corner */
+ *(u + nrc - 1) = *(u + nrc - nc - 2); /* lr corner */
+ *z = *(z + nc + 1); /* ul corner */
+ *(z + nc - 1) = *(z + 2 * nc - 2); /* ur corner */
+ *(z + nrc - nc) = *(z + nrc - 2 * nc + 1); /* ll corner */
+ *(z + nrc - 1) = *(z + nrc - nc - 2); /* lr corner */
+
+ /* this is an alternative for the corners using the mean value of the x and y closest points
+ * *u = 0.5 * ( *(u+1) + *(u+nc) );
+ * *(u+nc-1) = 0.5 * ( *(u+nc-2) + *(u+2*nc-1) );
+ * *(u+nrc-nc) = 0.5 * ( *(u+nrc-nc+1) + *(u+nrc-2*nc) );
+ * *(u+nrc-1) = 0.5 * ( *(u+nrc-2) +*(u+nrc-nc-1) );
+ * *z = 0.5 * ( *(z+1) + *(z+nc) );
+ * *(z+nc-1) = 0.5 * ( *(z+nc-2) + *(z+2*nc-1) );
+ * *(z+nrc-nc) = 0.5 * ( *(z+nrc-nc+1) + *(z+nrc-2*nc) );
+ * *(z+nrc-1) = 0.5 * ( *(z+nrc-2) +*(z+nrc-nc-1) );
+ */
+
+ /*** image edges ***/
+ for (i = 1; i < nc - 1; i++) { /* top and bottom edges */
+ *(u + i) = *(u + nc + i); /* t edge */
+ *(u + nrc - nc + i) = *(u + nrc - 2 * nc + i); /* b edge */
+ *(z + i) = *(z + nc + i); /* t edge */
+ *(z + nrc - nc + i) = *(z + nrc - 2 * nc + i); /* b edge */
+ }
+ for (j = 1; j < nr - 1; j++) { /* left and right edges */
+ jnc = j * nc;
+ *(u + jnc) = *(u + jnc + 1); /* r edge */
+ *(u + jnc + nc - 1) = *(u + jnc + nc - 2); /* l edge */
+ *(z + jnc) = *(z + jnc + 1); /* r edge */
+ *(z + jnc + nc - 1) = *(z + jnc + nc - 2); /* l edge */
+ }
+
+ old_u = (double *)malloc(sizeof(double) * nrc);
+ old_z = (double *)malloc(sizeof(double) * nrc);
+
+ for (j = 0; j < nr; j++) {
+ jnc = j * nc;
+ for (i = 0; i < nc; i++) {
+ *(old_u + jnc + i) = *(u + jnc + i);
+ *(old_z + jnc + i) = *(z + jnc + i);
+ }
+ }
+
+ for (j = 1; j < nr - 1; j++) {
+ jnc = j * nc;
+ for (i = 1; i < nc - 1; i++) {
+ u_hat =
+ *(old_z + jnc + i + 1) * *(old_z + jnc + i + 1) * *(old_u +
+ jnc + i +
+ 1) +
+ *(old_z + jnc + i - 1) * *(old_z + jnc + i - 1) * *(old_u +
+ jnc + i -
+ 1) +
+ *(old_z + jnc + nc + i) * *(old_z + jnc + nc + i) * *(old_u +
+ jnc +
+ nc +
+ i) +
+ *(old_z + jnc - nc + i) * *(old_z + jnc - nc + i) * *(old_u +
+ jnc -
+ nc + i);
+ z_star =
+ *(old_z + jnc + i + 1) * *(old_z + jnc + i + 1) + *(old_z +
+ jnc + i -
+ 1) *
+ *(old_z + jnc + i - 1) + *(old_z + jnc + nc + i) * *(old_z +
+ jnc +
+ nc + i) +
+ *(old_z + jnc - nc + i) * *(old_z + jnc - nc + i);
+
+ num = lambda * u_hat + h2 * *(g + jnc + i);
+ den = lambda * z_star + h2;
+
+ *(u + jnc + i) = num / den;
+
+ z_tilde =
+ *(old_z + jnc + i + 1) + *(old_z + jnc + i - 1) + *(old_z +
+ jnc + nc +
+ i) +
+ *(old_z + jnc - nc + i);
+ ux = *(old_u + jnc + i + 1) - *(old_u + jnc + i - 1);
+ uy = *(old_u + jnc + nc + i) - *(old_u + jnc - nc + i);
+
+ num = 4 * z_tilde + h2 * k2;
+ den = 16 + h2 * k2 + lambda * k * (ux * ux + uy * uy) / alpha;
+
+ *(z + jnc + i) = num / den;
+
+ *mxdf = max(*mxdf, fabs(*(old_u + jnc + i) - *(u + jnc + i)));
+ }
+ }
+ free(old_u);
+ free(old_z);
+ return;
+}
+
+/* ========================================================================== */
+/* ========================================================================== */
+
+
+/* ========================================================================== */
+/*!
+ * \brief MSK_O --- MUMFORD-SHAH with CURVATURE term (Jacobi method)
+ *
+ * Implements the Mumford-Shah variational model with curvature term for image segmentation.
+ *
+ * A smooth approximation of the input image is produced.
+ * The model explicitely detects the discontinuities of the output approximation
+ * while preserving the discontinuities from being smoothed.
+ *
+ * A discontinuity image is also produced.
+ * The values of the discontinuity image are close to 1 where the output approximation
+ * is "homogeneous", where the output approximation has discontinuities (edges) the
+ * values are close to 0.
+ *
+ * The curvature term prevents the discontinuities from being shortened too much when the parameter
+ * alpha is set to very high values.
+ *
+ *
+ * \param[in] *g (double; pointer) - The input image
+ *
+ * \param[in] *u (double; pointer) - The output approximation
+ *
+ * \param[in] *z (double; pointer) - The output discontinuity map
+ *
+ * \param[in] lambda (double) - The scale factor parameter [positive]
+ *
+ * \param[in] kepsilon (double) - The discontinuities thickness [positive]
+ *
+ * \param[in] alpha (double) - The elasticity parameter [positive]
+ *
+ * \param[in] beta (double) - The rigidity parameter [positive]
+ *
+ * \param[in] *mxdf (double; pointer) - The max difference in the approximation between two iteration steps
+ *
+ * \param[in] nr (int) - The number of rows in the input image
+ *
+ * \param[in] nc (int) - The number of colums in the input image
+ *
+ * \return void [shold be changed to int for error handling]
+ */
+
+/** MSK_O --- MUMFORD-SHAH WITH CURVATURE (Jacobi method) **/
+void msk_o(double *g, double *u, double *z, double lambda, double kepsilon,
+ double alpha, double beta, double *mxdf, int nr, int nc)
+{
+ int i, j, jnc, nrc;
+ double *old_u, *old_z; /* msk_o */
+ double u_tilde, ux, uy, zx, zy, z_tilde, z_hat, z_check; /* msk */
+ double a, b, c, d, e; /* msk */
+ double num, den;
+ double k2, k, h2, h = 1;
+ double ta;
+
+ nrc = nr * nc;
+ h2 = h * h;
+ k = kepsilon;
+ k2 = kepsilon * kepsilon;
+
+ /* Remove first image boundary using BC values */
+ *u = *(u + nc + 1); /* ul corner */
+ *(u + nc - 1) = *(u + 2 * nc - 2); /* ur corner */
+ *(u + nrc - nc) = *(u + nrc - 2 * nc + 1); /* ll corner */
+ *(u + nrc - 1) = *(u + nrc - nc - 2); /* lr corner */
+ *z = *(z + nc + 1); /* ul corner */
+ *(z + nc - 1) = *(z + 2 * nc - 2); /* ur corner */
+ *(z + nrc - nc) = *(z + nrc - 2 * nc + 1); /* ll corner */
+ *(z + nrc - 1) = *(z + nrc - nc - 2); /* lr corner */
+ for (i = 1; i < nc - 1; i++) { /* top and bottom edges */
+ *(u + i) = *(u + nc + i); /* t edge */
+ *(u + nrc - nc + i) = *(u + nrc - 2 * nc + i); /* b edge */
+ *(z + i) = *(z + nc + i); /* t edge */
+ *(z + nrc - nc + i) = *(z + nrc - 2 * nc + i); /* b edge */
+ }
+ for (j = 1; j < nr - 1; j++) { /* left and right edges */
+ jnc = j * nc;
+ *(u + jnc) = *(u + jnc + 1); /* r edge */
+ *(u + jnc + nc - 1) = *(u + jnc + nc - 2); /* l edge */
+ *(z + jnc) = *(z + jnc + 1); /* r edge */
+ *(z + jnc + nc - 1) = *(z + jnc + nc - 2); /* l edge */
+ }
+
+ /* Boundary conditions (Neumann) */
+
+ /*** image corners ***/
+ *(u + nc + 1) = *(u + 2 * nc + 2); /* ul corner */
+ *(u + 2 * nc - 2) = *(u + 3 * nc - 3); /* ur corner */
+ *(u + nrc - 2 * nc + 1) = *(u + nrc - 3 * nc + 2); /* ll corner */
+ *(u + nrc - nc - 2) = *(u + nrc - 2 * nc - 3); /* lr corner */
+ *(z + nc + 1) = *(z + 2 * nc + 2); /* ul corner */
+ *(z + 2 * nc - 2) = *(z + 3 * nc - 3); /* ur corner */
+ *(z + nrc - 2 * nc + 1) = *(z + nrc - 3 * nc + 2); /* ll corner */
+ *(z + nrc - nc - 2) = *(z + nrc - 2 * nc - 3); /* lr corner */
+
+ /* this is an alternative for the corners using the mean value of the x and y closest points
+ * *(u+nc+1) = 0.5 * ( *(u+nc+1) + *(u+2*nc+1) );
+ * *(u+2*nc-2) = 0.5 * ( *(u+2*nc-3) + *(u+3*nc-2) );
+ * *(u+nrc-2*nc+1) = 0.5 * ( *(u+nrc-2*nc+2) + *(u+nrc-3*nc+1) );
+ * *(u+nrc-nc-2) = 0.5 * ( *(u+nrc-nc-3) +*(u+nrc-2*nc-2) );
+ * *z = 0.5 * ( *(z+nc+1) + *(z+2*nc+1) );
+ * *(z+nc-1) = 0.5 * ( *(z+2*nc-3) + *(z+3*nc-2) );
+ * *(z+nrc-nc) = 0.5 * ( *(z+nrc-2*nc+2) + *(z+nrc-3*nc+1) );
+ * *(z+nrc-1) = 0.5 * ( *(z+nrc-nc-3) +*(z+nrc-2*nc-2) );
+ */
+
+ /*** image edges ***/
+ for (i = 2; i < nc - 2; i++) { /* top and bottom edges */
+ *(u + nc + i) = *(u + 2 * nc + i); /* t edge */
+ *(u + nrc - 2 * nc + i) = *(u + nrc - 3 * nc + i); /* b edge */
+ *(z + nc + i) = *(z + 2 * nc + i); /* t edge */
+ *(z + nrc - 2 * nc + i) = *(z + nrc - 3 * nc + i); /* b edge */
+ }
+ for (j = 2; j < nr - 2; j++) { /* left and right edges */
+ jnc = j * nc;
+ *(u + jnc + 1) = *(u + jnc + 2); /* r edge */
+ *(u + jnc + nc - 2) = *(u + jnc + nc - 3); /* l edge */
+ *(z + jnc + 1) = *(z + jnc + 2); /* r edge */
+ *(z + jnc + nc - 2) = *(z + jnc + nc - 3); /* l edge */
+ }
+
+ old_u = (double *)malloc(sizeof(double) * nrc);
+ old_z = (double *)malloc(sizeof(double) * nrc);
+
+ for (j = 0; j < nr; j++) {
+ jnc = j * nc;
+ for (i = 0; i < nc; i++) {
+ *(old_u + jnc + i) = *(u + jnc + i);
+ *(old_z + jnc + i) = *(z + jnc + i);
+ }
+ }
+
+ for (j = 2; j < nr - 2; j++) {
+ jnc = j * nc;
+ for (i = 2; i < nc - 2; i++) {
+
+ u_tilde =
+ *(old_u + jnc + i + 1) + *(old_u + jnc + i - 1) + *(old_u +
+ jnc + nc +
+ i) +
+ *(old_u + jnc - nc + i);
+ ux = *(old_u + jnc + i + 1) - *(old_u + jnc + i - 1);
+ uy = *(old_u + jnc + nc + i) - *(old_u + jnc - nc + i);
+ zx = *(old_z + jnc + i + 1) - *(old_z + jnc + i - 1);
+ zy = *(old_z + jnc + nc + i) - *(old_z + jnc - nc + i);
+ num =
+ 0.5 * lambda * *(old_z + jnc + i) * (zx * ux + zy * uy) +
+ lambda * *(old_z + jnc + i) * *(old_z + jnc + i) * u_tilde +
+ h2 * *(g + jnc + i);
+ den = 4 * lambda * *(old_z + jnc + i) * *(old_z + jnc + i) + h2;
+
+ *(u + jnc + i) = num / den;
+
+ z_tilde =
+ *(old_z + jnc + i + 1) + *(old_z + jnc + i - 1) + *(old_z +
+ jnc + nc +
+ i) +
+ *(old_z + jnc - nc + i);
+ z_hat =
+ *(old_z + jnc + nc + i + 1) + *(old_z + jnc - nc + i + 1) +
+ *(old_z + jnc + nc + i - 1) + *(old_z + jnc - nc + i - 1);
+ z_check =
+ *(old_z + jnc + i + 2) + *(old_z + jnc + i - 2) + *(old_z +
+ jnc +
+ 2 * nc +
+ i) +
+ *(old_z + jnc - 2 * nc + i);
+ a = alpha * z_tilde + 4 * alpha * h2 * k2 * *(old_z + jnc +
+ i) * *(old_z + jnc +
+ i) *
+ *(old_z + jnc + i);
+ b = 4 * beta * (8 * z_tilde - 2 * z_hat - z_check) / h2;
+ c = -16 * beta * k2 * (3 * *(old_z + jnc + i) *
+ *(old_z + jnc + i) + 1) * (z_tilde -
+ 4 * *(old_z +
+ jnc + i));
+ d = 64 * beta * k2 * *(old_z + jnc +
+ i) * (3 * *(old_z + jnc + i) * *(old_z +
+ jnc + i) -
+ 1);
+ e = 64 * beta * h2 * k2 * k2 * *(old_z + jnc + i) * *(old_z +
+ jnc +
+ i) *
+ *(old_z + jnc +
+ i) * (3 * *(old_z + jnc + i) * *(old_z + jnc + i) - 2);
+ num = a + b + c + d + e;
+ ta = a;
+
+ a = 4 * alpha +
+ 2 * alpha * h2 * k2 * (3 * *(old_z + jnc + i) *
+ *(old_z + jnc + i) - 1) +
+ 0.25 * lambda * k * (ux * ux + uy * uy);
+ b = 16 * beta * k2 * (h2 * k2 + 5 / (h2 * k2) - 4);
+ c = -12 * beta * k2 * (zx * zx + zy * zy);
+ d = -96 * beta * k2 * *(old_z + jnc + i) * (z_tilde -
+ 4 * *(z + jnc + i));
+ e = 192 * beta * k2 * *(old_z + jnc + i) * (1 - h2 * k2) +
+ 240 * beta * h2 * k2 * k2 * *(old_z + jnc + i) * *(old_z +
+ jnc +
+ i) *
+ *(old_z + jnc + i) * *(old_z + jnc + i);
+ den = a + b + c + d + e;
+
+ *(z + jnc + i) = num / den;
+
+ *mxdf = max(*mxdf, fabs(*(old_u + jnc + i) - *(u + jnc + i)));
+ }
+ }
+ free(old_u);
+ free(old_z);
+ return;
+}
+
+/* ========================================================================== */
+/* ========================================================================== */
+
+
+/* ========================================================================== */
+/*!
+ * \brief MS_T --- MUMFORD-SHAH (Gauss-Seidel method) - (Different approximation of "u" wrt MS_N)
+ *
+ * Implements the Mumford-Shah variational model for image segmentation.
+ *
+ * A smooth approximation of the input image is produced.
+ * The model explicitely detects the discontinuities of the output approximation
+ * while preserving the discontinuities from being smoothed.
+ *
+ * A discontinuity image is also produced.
+ * The values of the discontinuity image are close to 1 where the output approximation
+ * is "homogeneous", where the output approximation has discontinuities (edges) the
+ * values are close to 0.
+ *
+ * The approximation of "u" is different wrt the one used in MS_N.
+ *
+ *
+ * \param[in] *g (double; pointer) - The input image
+ *
+ * \param[in] *u (double; pointer) - The output approximation
+ *
+ * \param[in] *z (double; pointer) - The output discontinuity map
+ *
+ * \param[in] lambda (double) - The scale factor parameter [positive]
+ *
+ * \param[in] kepsilon (double) - The discontinuities thickness [positive]
+ *
+ * \param[in] alpha (double) - The elasticity parameter [positive]
+ *
+ * \param[in] *mxdf (double; pointer) - The max difference in the approximation between two iteration steps
+ *
+ * \param[in] nr (int) - The number of rows in the input image
+ *
+ * \param[in] nc (int) - The number of colums in the input image
+ *
+ * \return void [shold be changed to int for error handling]
+ */
+
+/** MS_T --- MUMFORD-SHAH (Gauss-Seidel method) - (Different approximation of "u" wrt MS_N) **/
+void ms_t(double *g, double *u, double *z, double lambda, double kepsilon,
+ double alpha, double *mxdf, int nr, int nc)
+{
+ int i, j, jnc, nrc;
+ double old_u; /* ms_n */
+ double old_t; /* ms_t */
+ double u_hat, z_star, z_tilde, ux, uy; /* ms */
+ double u_tilde, zx, zy; /* msk */
+ double num, den;
+ double k2, k, h2, h = 1;
+
+ nrc = nr * nc;
+ h2 = h * h;
+ k = kepsilon;
+ k2 = kepsilon * kepsilon;
+
+ /* Boundary conditions (Neumann) */
+ *u = *(u + nc + 1); /* ul corner */
+ *(u + nc - 1) = *(u + 2 * nc - 2); /* ur corner */
+ *(u + nrc - nc) = *(u + nrc - 2 * nc + 1); /* ll corner */
+ *(u + nrc - 1) = *(u + nrc - nc - 2); /* lr corner */
+ *z = *(z + nc + 1); /* ul corner */
+ *(z + nc - 1) = *(z + 2 * nc - 2); /* ur corner */
+ *(z + nrc - nc) = *(z + nrc - 2 * nc + 1); /* ll corner */
+ *(z + nrc - 1) = *(z + nrc - nc - 2); /* lr corner */
+ for (i = 1; i < nc - 1; i++) { /* top and bottom edges */
+ *(u + i) = *(u + nc + i); /* t edge */
+ *(u + nrc - nc + i) = *(u + nrc - 2 * nc + i); /* b edge */
+ *(z + i) = *(z + nc + i); /* t edge */
+ *(z + nrc - nc + i) = *(z + nrc - 2 * nc + i); /* b edge */
+ }
+ for (j = 1; j < nr - 1; j++) { /* left and right edges */
+ jnc = j * nc;
+ *(u + jnc) = *(u + jnc + 1); /* r edge */
+ *(u + jnc + nc - 1) = *(u + jnc + nc - 2); /* l edge */
+ *(z + jnc) = *(z + jnc + 1); /* r edge */
+ *(z + jnc + nc - 1) = *(z + jnc + nc - 2); /* l edge */
+ }
+ for (j = 1; j < nr - 1; j++) {
+ jnc = j * nc;
+ for (i = 1; i < nc - 1; i++) {
+ old_u = *(u + jnc + i);
+
+ /* March version and MS */
+ u_hat =
+ *(z + jnc + i + 1) * *(z + jnc + i + 1) * *(u + jnc + i + 1) +
+ *(z + jnc + i - 1) * *(z + jnc + i - 1) * *(u + jnc + i - 1) +
+ *(z + jnc + nc + i) * *(z + jnc + nc + i) * *(u + jnc + nc +
+ i) + *(z + jnc -
+ nc +
+ i) *
+ *(z + jnc - nc + i) * *(u + jnc - nc + i);
+ z_star =
+ *(z + jnc + i + 1) * *(z + jnc + i + 1) + *(z + jnc + i -
+ 1) * *(z + jnc +
+ i - 1) +
+ *(z + jnc + nc + i) * *(z + jnc + nc + i) + *(z + jnc - nc +
+ i) * *(z + jnc -
+ nc + i);
+ num = lambda * u_hat + h2 * *(g + jnc + i);
+ den = lambda * z_star + h2;
+ *(u + jnc + i) = num / den;
+
+ old_t = *(u + jnc + i);
+
+ /* math thesis version and MSK */
+ u_tilde =
+ *(u + jnc + i + 1) + *(u + jnc + i - 1) + *(u + jnc + nc +
+ i) + *(u + jnc -
+ nc + i);
+ ux = *(u + jnc + i + 1) - *(u + jnc + i - 1);
+ uy = *(u + jnc + nc + i) - *(u + jnc - nc + i);
+ zx = *(z + jnc + i + 1) - *(z + jnc + i - 1);
+ zy = *(z + jnc + nc + i) - *(z + jnc - nc + i);
+ num =
+ 0.5 * lambda * *(z + jnc + i) * (zx * ux + zy * uy) +
+ lambda * *(z + jnc + i) * *(z + jnc + i) * u_tilde +
+ h2 * *(g + jnc + i);
+ den = 4 * lambda * *(z + jnc + i) * *(z + jnc + i) + h2;
+ *(u + jnc + i) = num / den;
+
+ z_tilde =
+ *(z + jnc + i + 1) + *(z + jnc + i - 1) + *(z + jnc + nc +
+ i) + *(z + jnc -
+ nc + i);
+ num = 4 * z_tilde + h2 * k2;
+ den = 16 + h2 * k2 + lambda * k * (ux * ux + uy * uy) / alpha;
+ *(z + jnc + i) = num / den;
+
+ *mxdf = max(*mxdf, fabs(old_u - *(u + jnc + i)));
+ }
+ }
+ return;
+}
+
+/* ========================================================================== */
+/* ========================================================================== */
+
+
+/* ========================================================================== */
+/*!
+ * \brief MSK_T --- MUMFORD-SHAH with CURVATURE term (Gauss-Seidel method) - (Different approximation of "u" wrt MSK_N)
+ *
+ * Implements the Mumford-Shah variational model with curvature term for image segmentation.
+ *
+ * A smooth approximation of the input image is produced.
+ * The model explicitely detects the discontinuities of the output approximation
+ * while preserving the discontinuities from being smoothed.
+ *
+ * A discontinuity image is also produced.
+ * The values of the discontinuity image are close to 1 where the output approximation
+ * is "homogeneous", where the output approximation has discontinuities (edges) the
+ * values are close to 0.
+ *
+ * The curvature term prevents the discontinuities from being shortened too much when the parameter
+ * alpha is set to very high values.
+ *
+ * The approximation of "u" is different wrt the one used in MSK_N.
+ *
+ *
+ * \param[in] *g (double; pointer) - The input image
+ *
+ * \param[in] *u (double; pointer) - The output approximation
+ *
+ * \param[in] *z (double; pointer) - The output discontinuity map
+ *
+ * \param[in] lambda (double) - The scale factor parameter [positive]
+ *
+ * \param[in] kepsilon (double) - The discontinuities thickness [positive]
+ *
+ * \param[in] alpha (double) - The elasticity parameter [positive]
+ *
+ * \param[in] beta (double) - The rigidity parameter [positive]
+ *
+ * \param[in] *mxdf (double; pointer) - The max difference in the approximation between two iteration steps
+ *
+ * \param[in] nr (int) - The number of rows in the input image
+ *
+ * \param[in] nc (int) - The number of colums in the input image
+ *
+ * \return void [shold be changed to int for error handling]
+ */
+
+/** MSK_T --- TEST --- MUMFORD-SHAH with CURVATURE term (Gauss-Seidel method) - (Different approximation of "u" wrt MSK_N) **/
+void msk_t(double *g, double *u, double *z, double lambda, double kepsilon,
+ double alpha, double beta, double *mxdf, int nr, int nc)
+{
+ int i, j, jnc, nrc;
+ double old_u; /* ms_n */
+ double u_hat, z_star, z_tilde, ux, uy; /* ms */
+ double u_tilde, zx, zy, z_hat, z_check; /* msk */
+ double a, b, c, d, e; /* msk */
+ double num, den;
+ double k2, k, h2, h = 1;
+ double ta;
+
+ nrc = nr * nc;
+ h2 = h * h;
+ k = kepsilon;
+ k2 = kepsilon * kepsilon;
+
+ /* Remove first image boundary using BC values */
+ *u = *(u + nc + 1); /* ul corner */
+ *(u + nc - 1) = *(u + 2 * nc - 2); /* ur corner */
+ *(u + nrc - nc) = *(u + nrc - 2 * nc + 1); /* ll corner */
+ *(u + nrc - 1) = *(u + nrc - nc - 2); /* lr corner */
+ *z = *(z + nc + 1); /* ul corner */
+ *(z + nc - 1) = *(z + 2 * nc - 2); /* ur corner */
+ *(z + nrc - nc) = *(z + nrc - 2 * nc + 1); /* ll corner */
+ *(z + nrc - 1) = *(z + nrc - nc - 2); /* lr corner */
+ for (i = 1; i < nc - 1; i++) { /* top and bottom edges */
+ *(u + i) = *(u + nc + i); /* t edge */
+ *(u + nrc - nc + i) = *(u + nrc - 2 * nc + i); /* b edge */
+ *(z + i) = *(z + nc + i); /* t edge */
+ *(z + nrc - nc + i) = *(z + nrc - 2 * nc + i); /* b edge */
+ }
+ for (j = 1; j < nr - 1; j++) { /* left and right edges */
+ jnc = j * nc;
+ *(u + jnc) = *(u + jnc + 1); /* r edge */
+ *(u + jnc + nc - 1) = *(u + jnc + nc - 2); /* l edge */
+ *(z + jnc) = *(z + jnc + 1); /* r edge */
+ *(z + jnc + nc - 1) = *(z + jnc + nc - 2); /* l edge */
+ }
+
+ /* Boundary conditions (Neumann) */
+ *(u + nc + 1) = *(u + 2 * nc + 2); /* ul corner */
+ *(u + 2 * nc - 2) = *(u + 3 * nc - 3); /* ur corner */
+ *(u + nrc - 2 * nc + 1) = *(u + nrc - 3 * nc + 2); /* ll corner */
+ *(u + nrc - nc - 2) = *(u + nrc - 2 * nc - 3); /* lr corner */
+ *(z + nc + 1) = *(z + 2 * nc + 2); /* ul corner */
+ *(z + 2 * nc - 2) = *(z + 3 * nc - 3); /* ur corner */
+ *(z + nrc - 2 * nc + 1) = *(z + nrc - 3 * nc + 2); /* ll corner */
+ *(z + nrc - nc - 2) = *(z + nrc - 2 * nc - 3); /* lr corner */
+ for (i = 2; i < nc - 2; i++) { /* top and bottom edges */
+ *(u + nc + i) = *(u + 2 * nc + i); /* t edge */
+ *(u + nrc - 2 * nc + i) = *(u + nrc - 3 * nc + i); /* b edge */
+ *(z + nc + i) = *(z + 2 * nc + i); /* t edge */
+ *(z + nrc - 2 * nc + i) = *(z + nrc - 3 * nc + i); /* b edge */
+ }
+ for (j = 2; j < nr - 2; j++) { /* left and right edges */
+ jnc = j * nc;
+ *(u + jnc + 1) = *(u + jnc + 2); /* r edge */
+ *(u + jnc + nc - 2) = *(u + jnc + nc - 3); /* l edge */
+ *(z + jnc + 1) = *(z + jnc + 2); /* r edge */
+ *(z + jnc + nc - 2) = *(z + jnc + nc - 3); /* l edge */
+ }
+
+ for (j = 2; j < nr - 2; j++) {
+ jnc = j * nc;
+ for (i = 2; i < nc - 2; i++) {
+ old_u = *(u + jnc + i);
+
+ /* math thesis version and MKK */
+ u_tilde =
+ *(u + jnc + i + 1) + *(u + jnc + i - 1) + *(u + jnc + nc +
+ i) + *(u + jnc -
+ nc + i);
+ ux = *(u + jnc + i + 1) - *(u + jnc + i - 1);
+ uy = *(u + jnc + nc + i) - *(u + jnc - nc + i);
+ zx = *(z + jnc + i + 1) - *(z + jnc + i - 1);
+ zy = *(z + jnc + nc + i) - *(z + jnc - nc + i);
+ num =
+ 0.5 * lambda * *(z + jnc + i) * (zx * ux + zy * uy) +
+ lambda * *(z + jnc + i) * *(z + jnc + i) * u_tilde +
+ h2 * *(g + jnc + i);
+ den = 4 * lambda * *(z + jnc + i) * *(z + jnc + i) + h2;
+ *(u + jnc + i) = num / den;
+
+ /* March version and MS */
+ u_hat =
+ *(z + jnc + i + 1) * *(z + jnc + i + 1) * *(u + jnc + i + 1) +
+ *(z + jnc + i - 1) * *(z + jnc + i - 1) * *(u + jnc + i - 1) +
+ *(z + jnc + nc + i) * *(z + jnc + nc + i) * *(u + jnc + nc +
+ i) + *(z + jnc -
+ nc +
+ i) *
+ *(z + jnc - nc + i) * *(u + jnc - nc + i);
+ z_star =
+ *(z + jnc + i + 1) * *(z + jnc + i + 1) + *(z + jnc + i -
+ 1) * *(z + jnc +
+ i - 1) +
+ *(z + jnc + nc + i) * *(z + jnc + nc + i) + *(z + jnc - nc +
+ i) * *(z + jnc -
+ nc + i);
+ num = lambda * u_hat + h2 * *(g + jnc + i);
+ den = lambda * z_star + h2;
+ *(u + jnc + i) = num / den;
+
+ z_tilde =
+ *(z + jnc + i + 1) + *(z + jnc + i - 1) + *(z + jnc + nc +
+ i) + *(z + jnc -
+ nc + i);
+ z_hat =
+ *(z + jnc + nc + i + 1) + *(z + jnc - nc + i + 1) + *(z +
+ jnc +
+ nc + i -
+ 1) +
+ *(z + jnc - nc + i - 1);
+ z_check =
+ *(z + jnc + i + 2) + *(z + jnc + i - 2) + *(z + jnc + 2 * nc +
+ i) + *(z + jnc -
+ 2 * nc +
+ i);
+ a = alpha * z_tilde + 4 * alpha * h2 * k2 * *(z + jnc + i) * *(z +
+ jnc
+ +
+ i)
+ * *(z + jnc + i);
+ b = 4 * beta * (8 * z_tilde - 2 * z_hat - z_check) / h2;
+ c = -16 * beta * k2 * (3 * *(z + jnc + i) * *(z + jnc + i) +
+ 1) * (z_tilde - 4 * *(z + jnc + i));
+ d = 64 * beta * k2 * *(z + jnc +
+ i) * (3 * *(z + jnc + i) * *(z + jnc + i) -
+ 1);
+ e = 64 * beta * h2 * k2 * k2 * *(z + jnc + i) * *(z + jnc +
+ i) * *(z + jnc +
+ i) * (3 *
+ *(z
+ +
+ jnc
+ +
+ i)
+ *
+ *(z
+ +
+ jnc
+ +
+ i)
+ -
+ 2);
+ num = a + b + c + d + e;
+ ta = a;
+
+ a = 4 * alpha +
+ 2 * alpha * h2 * k2 * (3 * *(z + jnc + i) * *(z + jnc + i) -
+ 1) + 0.25 * lambda * k * (ux * ux +
+ uy * uy);
+ b = 16 * beta * k2 * (h2 * k2 + 5 / (h2 * k2) - 4);
+ c = -12 * beta * k2 * (zx * zx + zy * zy);
+ d = -96 * beta * k2 * *(z + jnc + i) * (z_tilde -
+ 4 * *(z + jnc + i));
+ e = 192 * beta * k2 * *(z + jnc + i) * (1 - h2 * k2) +
+ 240 * beta * h2 * k2 * k2 * *(z + jnc + i) * *(z + jnc +
+ i) * *(z +
+ jnc +
+ i) *
+ *(z + jnc + i);
+ den = a + b + c + d + e;
+ *(z + jnc + i) = num / den;
+
+ *mxdf = max(*mxdf, fabs(old_u - *(u + jnc + i)));
+ }
+ }
+ return;
+}
+
+/* ========================================================================== */
+/* ========================================================================== */
+
Added: grass-addons/grass7/raster/r.seg/varseg.h
===================================================================
--- grass-addons/grass7/raster/r.seg/varseg.h (rev 0)
+++ grass-addons/grass7/raster/r.seg/varseg.h 2015-06-09 12:51:48 UTC (rev 65425)
@@ -0,0 +1,18 @@
+/** MS_N --- MUMFORD-SHAH (Gauss-Seidel method) **/
+void ms_n(double *, double *, double *, double, double, double, double *, int, int);
+
+/** MSK_N --- MUMFORD-SHAH with CURVATURE term (Gauss-Seidel method) **/
+void msk_n(double *, double *, double *, double, double, double, double, double *, int, int);
+
+/** MS_O --- MUMFORD-SHAH (Jacobi method) **/
+void ms_o(double *, double *, double *, double, double, double, double *, int, int);
+
+/** MSK_O --- MUMFORD-SHAH with CURVATURE term (Jacobi method) **/
+void msk_o(double *, double *, double *, double, double, double, double, double *, int, int);
+
+/** MS_T --- TEST --- MUMFORD-SHAH (Gauss-Seidel method) - (Different approximation of "u" wrt MS_N) **/
+void ms_t(double *, double *, double *, double, double, double, double *, int, int);
+
+/** MSK_T --- TEST --- MUMFORD-SHAH with CURVATURE term (Gauss-Seidel method) - (Different approximation of "u" wrt MSK_N) **/
+void msk_t(double *, double *, double *, double, double, double, double, double *, int, int);
+
More information about the grass-commit
mailing list