[GRASS-SVN] r43303 - in grass-addons/raster: . r.seg r.seg/include
r.seg/include/Make r.seg/include/varseg r.seg/rseg r.seg/varseglib
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Aug 27 10:29:35 EDT 2010
Author: avitti
Date: 2010-08-27 14:29:35 +0000 (Fri, 27 Aug 2010)
New Revision: 43303
Added:
grass-addons/raster/r.seg/
grass-addons/raster/r.seg/Makefile
grass-addons/raster/r.seg/README
grass-addons/raster/r.seg/include/
grass-addons/raster/r.seg/include/Make/
grass-addons/raster/r.seg/include/Make/varseg.make
grass-addons/raster/r.seg/include/varseg/
grass-addons/raster/r.seg/include/varseg/varseg.h
grass-addons/raster/r.seg/rseg/
grass-addons/raster/r.seg/rseg/Makefile
grass-addons/raster/r.seg/rseg/description.html
grass-addons/raster/r.seg/rseg/main.c
grass-addons/raster/r.seg/varseglib/
grass-addons/raster/r.seg/varseglib/Makefile
grass-addons/raster/r.seg/varseglib/varseg.c
Log:
r.seg - module for image segmentation and discontinuity detection
Property changes on: grass-addons/raster/r.seg
___________________________________________________________________
Added: svn:ignore
+ *OBJ*
Added: grass-addons/raster/r.seg/Makefile
===================================================================
--- grass-addons/raster/r.seg/Makefile (rev 0)
+++ grass-addons/raster/r.seg/Makefile 2010-08-27 14:29:35 UTC (rev 43303)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+SUBDIRS = varseglib \
+ rseg
+
+include $(MODULE_TOPDIR)/include/Make/Dir.make
+
+default: subdirs
+
+clean: cleansubdirs
Property changes on: grass-addons/raster/r.seg/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:keywords
+ Author Date Id
Added: grass-addons/raster/r.seg/README
===================================================================
--- grass-addons/raster/r.seg/README (rev 0)
+++ grass-addons/raster/r.seg/README 2010-08-27 14:29:35 UTC (rev 43303)
@@ -0,0 +1,13 @@
+r.seg GRASS GIS module for image segmentation and edge detection
+ Alfonso Vitti <alfonso.vitti [at] ing.unitn.it>
+ see www.ing.unitn.it/~vittia/sw
+
+copy the "r.seg" directory in the "raster" directory of the GRASS source code directory
+
+from within the "r.seg" directory, as normal user run:
+make
+
+from the GRASS main source code directory, as root run:
+make install
+or, as normal user:
+sudo make install
Property changes on: grass-addons/raster/r.seg/README
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Author Date Id
Property changes on: grass-addons/raster/r.seg/include
___________________________________________________________________
Added: svn:ignore
+ *.tmp.html
Property changes on: grass-addons/raster/r.seg/include/Make
___________________________________________________________________
Added: svn:ignore
+ *OBJ*
Added: grass-addons/raster/r.seg/include/Make/varseg.make
===================================================================
--- grass-addons/raster/r.seg/include/Make/varseg.make (rev 0)
+++ grass-addons/raster/r.seg/include/Make/varseg.make 2010-08-27 14:29:35 UTC (rev 43303)
@@ -0,0 +1,3 @@
+SEG_LIBNAME = varseg
+SEGLIB = -l$(SEG_LIBNAME)
+SEGDEP = $(ARCH_LIBDIR)/$(LIB_PREFIX)$(SEG_LIBNAME)$(LIB_SUFFIX)
Property changes on: grass-addons/raster/r.seg/include/Make/varseg.make
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:keywords
+ Author Date Id
Property changes on: grass-addons/raster/r.seg/include/varseg
___________________________________________________________________
Added: svn:ignore
+ *OBJ*
Added: grass-addons/raster/r.seg/include/varseg/varseg.h
===================================================================
--- grass-addons/raster/r.seg/include/varseg/varseg.h (rev 0)
+++ grass-addons/raster/r.seg/include/varseg/varseg.h 2010-08-27 14:29:35 UTC (rev 43303)
@@ -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);
+
Property changes on: grass-addons/raster/r.seg/include/varseg/varseg.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:keywords
+ Author Date Id
Property changes on: grass-addons/raster/r.seg/rseg
___________________________________________________________________
Added: svn:ignore
+ *.tmp.html
Added: grass-addons/raster/r.seg/rseg/Makefile
===================================================================
--- grass-addons/raster/r.seg/rseg/Makefile (rev 0)
+++ grass-addons/raster/r.seg/rseg/Makefile 2010-08-27 14:29:35 UTC (rev 43303)
@@ -0,0 +1,19 @@
+### r.seg MAKEFILE ###
+
+MODULE_TOPDIR = ../../..
+
+PGM = r.seg
+
+LIBES = $(SEGLIB) $(GISLIB)
+DEPENDENCIES = $(SEGDEP) $(GISDEP)
+
+include ../include/Make/varseg.make
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+include $(MODULE_TOPDIR)/include/Make/Doxygen.make
+
+default: cmd
+
+#doxygen:
+DOXNAME = r.seg
Property changes on: grass-addons/raster/r.seg/rseg/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:keywords
+ Author Date Id
Added: grass-addons/raster/r.seg/rseg/description.html
===================================================================
--- grass-addons/raster/r.seg/rseg/description.html (rev 0)
+++ grass-addons/raster/r.seg/rseg/description.html 2010-08-27 14:29:35 UTC (rev 43303)
@@ -0,0 +1,120 @@
+<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 (r.null) any null data with the map average value (r.univar).
+<br><br>
+
+The segmentation depends on the parameters alpha and lambda:<br>
+- alpha controls how many discontinuities are allowed to exist.<br>
+- lambda controls the smoothness of the solution.<br>
+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:<br>
+- run the module with both alpha and lambda set to 1.0<br>
+- run the module with alpha set to 1.0 and different values for lambda<br>
+  e.g., 0.01, 0.1, 1, 10, 100<br>
+- 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<br>
+- see how the segmentations change and select the values that produce the result that best fits your requirements.
+<br><br>
+
+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:<br>
+- the resolution of the region matches the resolution of the input raster map;<br>
+- the boundaries of the region are lined up along the edges of the nearest cells in the input raster map.
+<br><br>
+
+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.itc.it/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 '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 c=differences
+#
+# for a better visualization of the output raster map <em>u_OF</em>, set its color table to:
+  r.colors z_OF c=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.zc.html">i.zc</a></em>,
+<em><a href="r.mfilter.html">r.mfilter</a></em>
+
+
+<h2>AUTHORS</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>
Property changes on: grass-addons/raster/r.seg/rseg/description.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: grass-addons/raster/r.seg/rseg/main.c
===================================================================
--- grass-addons/raster/r.seg/rseg/main.c (rev 0)
+++ grass-addons/raster/r.seg/rseg/main.c 2010-08-27 14:29:35 UTC (rev 43303)
@@ -0,0 +1,327 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.seg
+ *
+ * AUTHOR: Alfonso Vitti <alfonso.vitti ing.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-2010
+ *
+ * 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 <grass/gis.h>
+#include <grass/config.h>
+#include <grass/glocale.h>
+#include <grass/varseg.h>
+
+int main(int argc, char *argv[])
+{
+ struct Cell_head cellhd; /* region information */
+ 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 */
+ char *mapset; /* current mapset */
+ void *g_row; /* input row buffers */
+ void *out_u_row, *out_z_row; /* output row buffer */
+ 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) */
+ RASTER_MAP_TYPE dcell_data_type; /* GRASS raster data type (for DCELL raster) */
+ double *g, *u, *z; /* the variables for the actual computation */
+
+ 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_s; /* flag, k = use MSK instead of MS */
+
+
+ /* initialize GRASS environment */
+ G_gisinit(argv[0]); /*reads GRASS env */
+
+
+ /* initialize module */
+ module = G_define_module();
+ module->keywords = _("image segmentation, edge detection, 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 = _("lambda ___ 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 = _("alpha ___ 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 =
+ _("mxi ___ max 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 = _("tol ___ 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 =
+ _("kepsilon ___ 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 = _("beta ___ 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");
+
+
+ /* 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_cell2(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, GR_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, GR_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 = G_window_rows();
+ nc = G_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 = G_open_cell_old(in_g, mapset)) < 0)
+ G_fatal_error(_("cannot open raster file [%s]"), in_g);
+ if (G_get_cellhd(in_g, mapset, &cellhd) < 0)
+ G_fatal_error(_("cannot read [%s] header"), in_g);
+
+
+ /* allocate the buffer for storing the values of the input raster map */
+ g_row = G_allocate_raster_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;
+
+ if (G_get_raster_row(g_fd, g_row, j, dcell_data_type) < 0)
+ G_fatal_error(_("Cannot read from [%s] raster,"), in_g);
+
+ 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 */
+ G_close_cell(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 = G_open_raster_new(out_u, dcell_data_type)) < 0)
+ G_fatal_error(_("cannot open raster file [%s]"), out_u);
+ if ((out_z_fd = G_open_raster_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 = G_allocate_raster_buf(dcell_data_type);
+ out_z_row = G_allocate_raster_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);
+ }
+ if (G_put_raster_row(out_u_fd, out_u_row, dcell_data_type) < 0)
+ G_fatal_error(_("cannot write [%s] raster"), out_u);
+ if (G_put_raster_row(out_z_fd, out_z_row, dcell_data_type) < 0)
+ G_fatal_error(_("cannot write [%s] raster"), out_z);
+ }
+
+
+ /* close the input raster map and free memory */
+ G_close_cell(out_u_fd);
+ G_free(out_u_row);
+ G_close_cell(out_z_fd);
+ G_free(out_z_row);
+ G_free(g);
+ G_free(u);
+ G_free(z);
+
+
+ /* write history file */
+ G_short_history(out_u, "raster", &history);
+ G_command_history(&history);
+ sprintf(history.edhist[3], "iterations = %i", iter);
+ history.edlinecnt = 4;
+ G_write_history(out_u, &history);
+ G_write_history(out_z, &history);
+
+
+ /* exit */
+ exit(EXIT_SUCCESS);
+}
+
Property changes on: grass-addons/raster/r.seg/rseg/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Property changes on: grass-addons/raster/r.seg/varseglib
___________________________________________________________________
Added: svn:ignore
+ *.tmp.html
Added: grass-addons/raster/r.seg/varseglib/Makefile
===================================================================
--- grass-addons/raster/r.seg/varseglib/Makefile (rev 0)
+++ grass-addons/raster/r.seg/varseglib/Makefile 2010-08-27 14:29:35 UTC (rev 43303)
@@ -0,0 +1,21 @@
+### varseglib MAKEFILE ###
+
+MODULE_TOPDIR = ../../..
+
+LIB_NAME = $(SEG_LIBNAME)
+
+LIB_OBJS = varseg.o
+
+include ../include/Make/varseg.make
+
+include $(MODULE_TOPDIR)/include/Make/Lib.make
+
+include $(MODULE_TOPDIR)/include/Make/Doxygen.make
+
+default: $(ARCH_INCDIR)/varseg.h lib
+
+$(ARCH_INCDIR)/varseg.h: ../include/varseg/varseg.h
+ $(INSTALL) -m 644 ../include/varseg/varseg.h $(ARCH_INCDIR)/varseg.h
+
+#doxygen:
+DOXNAME = varseg
Property changes on: grass-addons/raster/r.seg/varseglib/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:keywords
+ Author Date Id
Added: grass-addons/raster/r.seg/varseglib/varseg.c
===================================================================
--- grass-addons/raster/r.seg/varseglib/varseg.c (rev 0)
+++ grass-addons/raster/r.seg/varseglib/varseg.c 2010-08-27 14:29:35 UTC (rev 43303)
@@ -0,0 +1,1124 @@
+
+/****************************************************************************
+ *
+ * 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 <grass/varseg.h> /* change this line according to $(ARCH_INCDIR) when using this library outside GRASS */
+
+#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;
+}
+
+/* ========================================================================== */
+/* ========================================================================== */
+
Property changes on: grass-addons/raster/r.seg/varseglib/varseg.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
More information about the grass-commit
mailing list