diff options
Diffstat (limited to 'plug-ins/common/contrast-retinex.c')
| -rw-r--r-- | plug-ins/common/contrast-retinex.c | 840 |
1 files changed, 840 insertions, 0 deletions
diff --git a/plug-ins/common/contrast-retinex.c b/plug-ins/common/contrast-retinex.c new file mode 100644 index 0000000..f48e5c9 --- /dev/null +++ b/plug-ins/common/contrast-retinex.c @@ -0,0 +1,840 @@ +/* GIMP - The GNU Image Manipulation Program + * Copyright (C) 1995 Spencer Kimball and Peter Mattis + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <https://www.gnu.org/licenses/>. + */ + +#include "config.h" + +#include <string.h> + +#include "libgimp/gimp.h" +#include "libgimp/gimpui.h" + +#include "libgimp/stdplugins-intl.h" + + +#define PLUG_IN_PROC "plug-in-retinex" +#define PLUG_IN_BINARY "contrast-retinex" +#define PLUG_IN_ROLE "gimp-contrast-retinex" +#define MAX_RETINEX_SCALES 8 +#define MIN_GAUSSIAN_SCALE 16 +#define MAX_GAUSSIAN_SCALE 250 +#define SCALE_WIDTH 150 +#define ENTRY_WIDTH 4 + + +typedef struct +{ + gint scale; + gint nscales; + gint scales_mode; + gfloat cvar; +} RetinexParams; + +typedef enum +{ + filter_uniform, + filter_low, + filter_high +} FilterMode; + +/* + Definit comment sont repartis les + differents filtres en fonction de + l'echelle (~= ecart type de la gaussienne) + */ +#define RETINEX_UNIFORM 0 +#define RETINEX_LOW 1 +#define RETINEX_HIGH 2 + +static gfloat RetinexScales[MAX_RETINEX_SCALES]; + +typedef struct +{ + gint N; + gfloat sigma; + gdouble B; + gdouble b[4]; +} gauss3_coefs; + + +/* + * Declare local functions. + */ +static void query (void); +static void run (const gchar *name, + gint nparams, + const GimpParam *param, + gint *nreturn_vals, + GimpParam **return_vals); + +/* Gimp */ +static gboolean retinex_dialog (gint32 drawable_ID); +static void retinex (gint32 drawable_ID, + GimpPreview *preview); +static void retinex_preview (gpointer drawable_ID, + GimpPreview *preview); + +static void retinex_scales_distribution (gfloat *scales, + gint nscales, + gint mode, + gint s); + +static void compute_mean_var (gfloat *src, + gfloat *mean, + gfloat *var, + gint size, + gint bytes); +/* + * Gauss + */ +static void compute_coefs3 (gauss3_coefs *c, + gfloat sigma); + +static void gausssmooth (gfloat *in, + gfloat *out, + gint size, + gint rowtride, + gauss3_coefs *c); + +/* + * MSRCR = MultiScale Retinex with Color Restoration + */ +static void MSRCR (guchar *src, + gint width, + gint height, + gint bytes, + gboolean preview_mode); + + +/* + * Private variables. + */ +static RetinexParams rvals = +{ + 240, /* Scale */ + 3, /* Scales */ + RETINEX_UNIFORM, /* Echelles reparties uniformement */ + 1.2 /* A voir */ +}; + +static GimpPlugInInfo PLUG_IN_INFO = +{ + NULL, /* init_proc */ + NULL, /* quit_proc */ + query, /* query_proc */ + run, /* run_proc */ +}; + +MAIN () + +static void +query (void) +{ + static const GimpParamDef args[] = + { + { GIMP_PDB_INT32, "run-mode", "The run mode { RUN-INTERACTIVE (0), RUN-NONINTERACTIVE (1) }" }, + { GIMP_PDB_IMAGE, "image", "Input image (unused)" }, + { GIMP_PDB_DRAWABLE, "drawable", "Input drawable" }, + { GIMP_PDB_INT32, "scale", "Biggest scale value" }, + { GIMP_PDB_INT32, "nscales", "Number of scales" }, + { GIMP_PDB_INT32, "scales-mode", "Retinex distribution through scales" }, + { GIMP_PDB_FLOAT, "cvar", "Variance value" } + }; + + gimp_install_procedure (PLUG_IN_PROC, + N_("Enhance contrast using the Retinex method"), + "The Retinex Image Enhancement Algorithm is an " + "automatic image enhancement method that enhances " + "a digital image in terms of dynamic range " + "compression, color independence from the spectral " + "distribution of the scene illuminant, and " + "color/lightness rendition.", + "Fabien Pelisson", + "Fabien Pelisson", + "2003", + N_("Retine_x..."), + "RGB*", + GIMP_PLUGIN, + G_N_ELEMENTS (args), 0, + args, NULL); + + gimp_plugin_menu_register (PLUG_IN_PROC, "<Image>/Colors/Tone Mapping"); +} + +static void +run (const gchar *name, + gint nparams, + const GimpParam *param, + gint *nreturn_vals, + GimpParam **return_vals) +{ + static GimpParam values[1]; + GimpRunMode run_mode; + gint32 drawable_ID; + GimpPDBStatusType status = GIMP_PDB_SUCCESS; + gint x, y, width, height; + + INIT_I18N (); + gegl_init (NULL, NULL); + + *nreturn_vals = 1; + *return_vals = values; + + values[0].type = GIMP_PDB_STATUS; + values[0].data.d_status = status; + + run_mode = param[0].data.d_int32; + drawable_ID = param[2].data.d_drawable; + + if (! gimp_drawable_mask_intersect (drawable_ID, + &x, &y, &width, &height) || + width < MIN_GAUSSIAN_SCALE || + height < MIN_GAUSSIAN_SCALE) + { + status = GIMP_PDB_EXECUTION_ERROR; + values[0].data.d_status = status; + return; + } + + switch (run_mode) + { + case GIMP_RUN_INTERACTIVE: + /* Possibly retrieve data */ + gimp_get_data (PLUG_IN_PROC, &rvals); + + /* First acquire information with a dialog */ + if (! retinex_dialog (drawable_ID)) + return; + break; + + case GIMP_RUN_NONINTERACTIVE: + /* Make sure all the arguments are there! */ + if (nparams != 7) + { + status = GIMP_PDB_CALLING_ERROR; + } + else + { + rvals.scale = (param[3].data.d_int32); + rvals.nscales = (param[4].data.d_int32); + rvals.scales_mode = (param[5].data.d_int32); + rvals.cvar = (param[6].data.d_float); + } + break; + + case GIMP_RUN_WITH_LAST_VALS: + gimp_get_data (PLUG_IN_PROC, &rvals); + break; + + default: + break; + } + + if (status == GIMP_PDB_SUCCESS && + (gimp_drawable_is_rgb (drawable_ID))) + { + gimp_progress_init (_("Retinex")); + + retinex (drawable_ID, NULL); + + if (run_mode != GIMP_RUN_NONINTERACTIVE) + gimp_displays_flush (); + + /* Store data */ + if (run_mode == GIMP_RUN_INTERACTIVE) + gimp_set_data (PLUG_IN_PROC, &rvals, sizeof (RetinexParams)); + } + else + { + status = GIMP_PDB_EXECUTION_ERROR; + } + + values[0].data.d_status = status; +} + + +static gboolean +retinex_dialog (gint32 drawable_ID) +{ + GtkWidget *dialog; + GtkWidget *main_vbox; + GtkWidget *preview; + GtkWidget *table; + GtkWidget *combo; + GtkObject *adj; + gboolean run; + + gimp_ui_init (PLUG_IN_BINARY, FALSE); + + dialog = gimp_dialog_new (_("Retinex Image Enhancement"), PLUG_IN_ROLE, + NULL, 0, + gimp_standard_help_func, PLUG_IN_PROC, + + _("_Cancel"), GTK_RESPONSE_CANCEL, + _("_OK"), GTK_RESPONSE_OK, + + NULL); + + gtk_dialog_set_alternative_button_order (GTK_DIALOG (dialog), + GTK_RESPONSE_OK, + GTK_RESPONSE_CANCEL, + -1); + + gimp_window_set_transient (GTK_WINDOW (dialog)); + + main_vbox = gtk_box_new (GTK_ORIENTATION_VERTICAL, 12); + gtk_container_set_border_width (GTK_CONTAINER (main_vbox), 12); + gtk_box_pack_start (GTK_BOX (gtk_dialog_get_content_area (GTK_DIALOG (dialog))), + main_vbox, TRUE, TRUE, 0); + gtk_widget_show (main_vbox); + + preview = gimp_zoom_preview_new_from_drawable_id (drawable_ID); + gtk_box_pack_start (GTK_BOX (main_vbox), preview, TRUE, TRUE, 0); + gtk_widget_show (preview); + + g_signal_connect_swapped (preview, "invalidated", + G_CALLBACK (retinex_preview), + GINT_TO_POINTER (drawable_ID)); + + table = gtk_table_new (4, 3, FALSE); + gtk_table_set_col_spacings (GTK_TABLE (table), 6); + gtk_table_set_row_spacings (GTK_TABLE (table), 6); + gtk_box_pack_start (GTK_BOX (main_vbox), table, FALSE, FALSE, 0); + gtk_widget_show (table); + + combo = gimp_int_combo_box_new (_("Uniform"), filter_uniform, + _("Low"), filter_low, + _("High"), filter_high, + NULL); + + gimp_int_combo_box_connect (GIMP_INT_COMBO_BOX (combo), rvals.scales_mode, + G_CALLBACK (gimp_int_combo_box_get_active), + &rvals.scales_mode); + g_signal_connect_swapped (combo, "changed", + G_CALLBACK (gimp_preview_invalidate), + preview); + + gimp_table_attach_aligned (GTK_TABLE (table), 0, 0, + _("_Level:"), 0.0, 0.5, + combo, 2, FALSE); + gtk_widget_show (combo); + + adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 1, + _("_Scale:"), SCALE_WIDTH, ENTRY_WIDTH, + rvals.scale, + MIN_GAUSSIAN_SCALE, MAX_GAUSSIAN_SCALE, 1, 1, 0, + TRUE, 0, 0, NULL, NULL); + + g_signal_connect (adj, "value-changed", + G_CALLBACK (gimp_int_adjustment_update), + &rvals.scale); + g_signal_connect_swapped (adj, "value-changed", + G_CALLBACK (gimp_preview_invalidate), + preview); + + adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 2, + _("Scale _division:"), SCALE_WIDTH, ENTRY_WIDTH, + rvals.nscales, + 0, MAX_RETINEX_SCALES, 1, 1, 0, + TRUE, 0, 0, NULL, NULL); + + g_signal_connect (adj, "value-changed", + G_CALLBACK (gimp_int_adjustment_update), + &rvals.nscales); + g_signal_connect_swapped (adj, "value-changed", + G_CALLBACK (gimp_preview_invalidate), + preview); + + adj = gimp_scale_entry_new (GTK_TABLE (table), 0, 3, + _("Dy_namic:"), SCALE_WIDTH, ENTRY_WIDTH, + rvals.cvar, 0, 4, 0.1, 0.1, 1, + TRUE, 0, 0, NULL, NULL); + + g_signal_connect (adj, "value-changed", + G_CALLBACK (gimp_float_adjustment_update), + &rvals.cvar); + g_signal_connect_swapped (adj, "value-changed", + G_CALLBACK (gimp_preview_invalidate), + preview); + + gtk_widget_show (dialog); + + run = (gimp_dialog_run (GIMP_DIALOG (dialog)) == GTK_RESPONSE_OK); + + gtk_widget_destroy (dialog); + + return run; +} + +/* + * Applies the algorithm + */ +static void +retinex (gint32 drawable_ID, + GimpPreview *preview) +{ + GeglBuffer *src_buffer; + GeglBuffer *dest_buffer; + const Babl *format; + guchar *src = NULL; + guchar *psrc = NULL; + gint x, y, width, height; + gint size, bytes; + + /* + * Get the size of the current image or its selection. + */ + if (preview) + { + src = gimp_zoom_preview_get_source (GIMP_ZOOM_PREVIEW (preview), + &width, &height, &bytes); + } + else + { + if (! gimp_drawable_mask_intersect (drawable_ID, + &x, &y, &width, &height)) + return; + + if (gimp_drawable_has_alpha (drawable_ID)) + format = babl_format ("R'G'B'A u8"); + else + format = babl_format ("R'G'B' u8"); + + bytes = babl_format_get_bytes_per_pixel (format); + + /* Allocate memory */ + size = width * height * bytes; + src = g_try_malloc (sizeof (guchar) * size); + + if (src == NULL) + { + g_warning ("Failed to allocate memory"); + return; + } + + memset (src, 0, sizeof (guchar) * size); + + /* Fill allocated memory with pixel data */ + src_buffer = gimp_drawable_get_buffer (drawable_ID); + + gegl_buffer_get (src_buffer, GEGL_RECTANGLE (x, y, width, height), 1.0, + format, src, + GEGL_AUTO_ROWSTRIDE, GEGL_ABYSS_NONE); + } + + /* + Algorithm for Multi-scale Retinex with color Restoration (MSRCR). + */ + psrc = src; + MSRCR (psrc, width, height, bytes, preview != NULL); + + if (preview) + { + gimp_preview_draw_buffer (preview, psrc, width * bytes); + } + else + { + dest_buffer = gimp_drawable_get_shadow_buffer (drawable_ID); + + gegl_buffer_set (dest_buffer, GEGL_RECTANGLE (x, y, width, height), 0, + format, psrc, + GEGL_AUTO_ROWSTRIDE); + + g_object_unref (src_buffer); + g_object_unref (dest_buffer); + + gimp_progress_update (1.0); + + gimp_drawable_merge_shadow (drawable_ID, TRUE); + gimp_drawable_update (drawable_ID, x, y, width, height); + } + + g_free (src); +} + +static void +retinex_preview (gpointer drawable_ID, + GimpPreview *preview) +{ + retinex (GPOINTER_TO_INT (drawable_ID), preview); +} + +/* + * calculate scale values for desired distribution. + */ +static void +retinex_scales_distribution (gfloat *scales, + gint nscales, + gint mode, + gint s) +{ + if (nscales == 1) + { /* For one filter we choose the median scale */ + scales[0] = (gint) s / 2; + } + else if (nscales == 2) + { /* For two filters we choose the median and maximum scale */ + scales[0] = (gint) s / 2; + scales[1] = (gint) s; + } + else + { + gfloat size_step = (gfloat) s / (gfloat) nscales; + gint i; + + switch(mode) + { + case RETINEX_UNIFORM: + for(i = 0; i < nscales; ++i) + scales[i] = 2. + (gfloat) i * size_step; + break; + + case RETINEX_LOW: + size_step = (gfloat) log(s - 2.0) / (gfloat) nscales; + for (i = 0; i < nscales; ++i) + scales[i] = 2. + pow (10, (i * size_step) / log (10)); + break; + + case RETINEX_HIGH: + size_step = (gfloat) log(s - 2.0) / (gfloat) nscales; + for (i = 0; i < nscales; ++i) + scales[i] = s - pow (10, (i * size_step) / log (10)); + break; + + default: + break; + } + } +} + +/* + * Calculate the coefficients for the recursive filter algorithm + * Fast Computation of gaussian blurring. + */ +static void +compute_coefs3 (gauss3_coefs *c, gfloat sigma) +{ + /* + * Papers: "Recursive Implementation of the gaussian filter.", + * Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995. + * formula: 11b computation of q + * 8c computation of b0..b1 + * 10 alpha is normalization constant B + */ + gfloat q, q2, q3; + + if (sigma >= 2.5) + { + q = 0.98711 * sigma - 0.96330; + } + else if ((sigma >= 0.5) && (sigma < 2.5)) + { + q = 3.97156 - 4.14554 * (gfloat) sqrt ((double) 1 - 0.26891 * sigma); + } + else + { + q = 0.1147705018520355224609375; + } + + q2 = q * q; + q3 = q * q2; + c->b[0] = (1.57825+(2.44413*q)+(1.4281 *q2)+(0.422205*q3)); + c->b[1] = ( (2.44413*q)+(2.85619*q2)+(1.26661 *q3)); + c->b[2] = ( -((1.4281*q2)+(1.26661 *q3))); + c->b[3] = ( (0.422205*q3)); + c->B = 1.0-((c->b[1]+c->b[2]+c->b[3])/c->b[0]); + c->sigma = sigma; + c->N = 3; + +/* + g_printerr ("q %f\n", q); + g_printerr ("q2 %f\n", q2); + g_printerr ("q3 %f\n", q3); + g_printerr ("c->b[0] %f\n", c->b[0]); + g_printerr ("c->b[1] %f\n", c->b[1]); + g_printerr ("c->b[2] %f\n", c->b[2]); + g_printerr ("c->b[3] %f\n", c->b[3]); + g_printerr ("c->B %f\n", c->B); + g_printerr ("c->sigma %f\n", c->sigma); + g_printerr ("c->N %d\n", c->N); +*/ +} + +static void +gausssmooth (gfloat *in, gfloat *out, gint size, gint rowstride, gauss3_coefs *c) +{ + /* + * Papers: "Recursive Implementation of the gaussian filter.", + * Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995. + * formula: 9a forward filter + * 9b backward filter + * fig7 algorithm + */ + gint i,n, bufsize; + gfloat *w1,*w2; + + /* forward pass */ + bufsize = size+3; + size -= 1; + w1 = (gfloat *) g_try_malloc (bufsize * sizeof (gfloat)); + w2 = (gfloat *) g_try_malloc (bufsize * sizeof (gfloat)); + w1[0] = in[0]; + w1[1] = in[0]; + w1[2] = in[0]; + for ( i = 0 , n=3; i <= size ; i++, n++) + { + w1[n] = (gfloat)(c->B*in[i*rowstride] + + ((c->b[1]*w1[n-1] + + c->b[2]*w1[n-2] + + c->b[3]*w1[n-3] ) / c->b[0])); + } + + /* backward pass */ + w2[size+1]= w1[size+3]; + w2[size+2]= w1[size+3]; + w2[size+3]= w1[size+3]; + for (i = size, n = i; i >= 0; i--, n--) + { + w2[n]= out[i * rowstride] = (gfloat)(c->B*w1[n+3] + + ((c->b[1]*w2[n+1] + + c->b[2]*w2[n+2] + + c->b[3]*w2[n+3] ) / c->b[0])); + } + + g_free (w1); + g_free (w2); +} + +/* + * This function is the heart of the algo. + * (a) Filterings at several scales and sumarize the results. + * (b) Calculation of the final values. + */ +static void +MSRCR (guchar *src, gint width, gint height, gint bytes, gboolean preview_mode) +{ + + gint scale,row,col; + gint i,j; + gint size; + gint channel; + guchar *psrc = NULL; /* backup pointer for src buffer */ + gfloat *dst = NULL; /* float buffer for algorithm */ + gfloat *pdst = NULL; /* backup pointer for float buffer */ + gfloat *in, *out; + gint channelsize; /* Float memory cache for one channel */ + gfloat weight; + gauss3_coefs coef; + gfloat mean, var; + gfloat mini, range, maxi; + gfloat alpha; + gfloat gain; + gfloat offset; + gdouble max_preview = 0.0; + + if (!preview_mode) + { + gimp_progress_init (_("Retinex: filtering")); + max_preview = 3 * rvals.nscales; + } + + /* Allocate all the memory needed for algorithm*/ + size = width * height * bytes; + dst = g_try_malloc (size * sizeof (gfloat)); + if (dst == NULL) + { + g_warning ("Failed to allocate memory"); + return; + } + memset (dst, 0, size * sizeof (gfloat)); + + channelsize = (width * height); + in = (gfloat *) g_try_malloc (channelsize * sizeof (gfloat)); + if (in == NULL) + { + g_free (dst); + g_warning ("Failed to allocate memory"); + return; /* do some clever stuff */ + } + + out = (gfloat *) g_try_malloc (channelsize * sizeof (gfloat)); + if (out == NULL) + { + g_free (in); + g_free (dst); + g_warning ("Failed to allocate memory"); + return; /* do some clever stuff */ + } + + + /* + Calculate the scales of filtering according to the + number of filter and their distribution. + */ + + retinex_scales_distribution (RetinexScales, + rvals.nscales, rvals.scales_mode, rvals.scale); + + /* + Filtering according to the various scales. + Summerize the results of the various filters according to a + specific weight(here equivalent for all). + */ + weight = 1./ (gfloat) rvals.nscales; + + /* + The recursive filtering algorithm needs different coefficients according + to the selected scale (~ = standard deviation of Gaussian). + */ + for (channel = 0; channel < 3; channel++) + { + gint pos; + + for (i = 0, pos = channel; i < channelsize ; i++, pos += bytes) + { + /* 0-255 => 1-256 */ + in[i] = (gfloat)(src[pos] + 1.0); + } + for (scale = 0; scale < rvals.nscales; scale++) + { + compute_coefs3 (&coef, RetinexScales[scale]); + /* + * Filtering (smoothing) Gaussian recursive. + * + * Filter rows first + */ + for (row=0 ;row < height; row++) + { + pos = row * width; + gausssmooth (in + pos, out + pos, width, 1, &coef); + } + + memcpy(in, out, channelsize * sizeof(gfloat)); + memset(out, 0 , channelsize * sizeof(gfloat)); + + /* + * Filtering (smoothing) Gaussian recursive. + * + * Second columns + */ + for (col=0; col < width; col++) + { + pos = col; + gausssmooth(in + pos, out + pos, height, width, &coef); + } + + /* + Summarize the filtered values. + In fact one calculates a ratio between the original values and the filtered values. + */ + for (i = 0, pos = channel; i < channelsize; i++, pos += bytes) + { + dst[pos] += weight * (log (src[pos] + 1.) - log (out[i])); + } + + if (!preview_mode) + gimp_progress_update ((channel * rvals.nscales + scale) / + max_preview); + } + } + g_free(in); + g_free(out); + + /* + Final calculation with original value and cumulated filter values. + The parameters gain, alpha and offset are constants. + */ + /* Ci(x,y)=log[a Ii(x,y)]-log[ Ei=1-s Ii(x,y)] */ + + alpha = 128.; + gain = 1.; + offset = 0.; + + for (i = 0; i < size; i += bytes) + { + gfloat logl; + + psrc = src+i; + pdst = dst+i; + + logl = log((gfloat)psrc[0] + (gfloat)psrc[1] + (gfloat)psrc[2] + 3.); + + pdst[0] = gain * ((log(alpha * (psrc[0]+1.)) - logl) * pdst[0]) + offset; + pdst[1] = gain * ((log(alpha * (psrc[1]+1.)) - logl) * pdst[1]) + offset; + pdst[2] = gain * ((log(alpha * (psrc[2]+1.)) - logl) * pdst[2]) + offset; + } + +/* if (!preview_mode) + gimp_progress_update ((2.0 + (rvals.nscales * 3)) / + ((rvals.nscales * 3) + 3));*/ + + /* + Adapt the dynamics of the colors according to the statistics of the first and second order. + The use of the variance makes it possible to control the degree of saturation of the colors. + */ + pdst = dst; + + compute_mean_var (pdst, &mean, &var, size, bytes); + mini = mean - rvals.cvar*var; + maxi = mean + rvals.cvar*var; + range = maxi - mini; + + if (!range) + range = 1.0; + + for (i = 0; i < size; i+= bytes) + { + psrc = src + i; + pdst = dst + i; + + for (j = 0 ; j < 3 ; j++) + { + gfloat c = 255 * ( pdst[j] - mini ) / range; + + psrc[j] = (guchar) CLAMP (c, 0, 255); + } + } + + g_free (dst); +} + +/* + * Calculate the average and variance in one go. + */ +static void +compute_mean_var (gfloat *src, gfloat *mean, gfloat *var, gint size, gint bytes) +{ + gfloat vsquared; + gint i,j; + gfloat *psrc; + + vsquared = 0; + *mean = 0; + for (i = 0; i < size; i+= bytes) + { + psrc = src+i; + for (j = 0 ; j < 3 ; j++) + { + *mean += psrc[j]; + vsquared += psrc[j] * psrc[j]; + } + } + + *mean /= (gfloat) size; /* mean */ + vsquared /= (gfloat) size; /* mean (x^2) */ + *var = ( vsquared - (*mean * *mean) ); + *var = sqrt(*var); /* var */ +} |