/*
 *  $Id: dune_synth.c 25200 2023-01-10 13:33:03Z yeti-dn $
 *  Copyright (C) 2022 David Necas (Yeti).
 *  E-mail: yeti@gwyddion.net.
 *
 *  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 2 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, write to the
 *  Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include "config.h"
#include <string.h>
#include <gtk/gtk.h>
#include <libgwyddion/gwymacros.h>
#include <libgwyddion/gwymath.h>
#include <libgwyddion/gwyrandgenset.h>
#include <libprocess/stats.h>
#include <libprocess/filters.h>
#include <libgwydgets/gwystock.h>
#include <libgwymodule/gwymodule-process.h>
#include <app/gwyapp.h>
#include <app/gwymoduleutils-synth.h>
#include "preview.h"
#include "libgwyddion/gwyomp.h"

#define RUN_MODES (GWY_RUN_IMMEDIATE | GWY_RUN_INTERACTIVE)

/* Cannot change this without losing reproducibility again! */
enum {
    NRANDOM_GENERATORS = 24
};

enum {
    PARAM_COVERAGE,
    PARAM_NITERS,
    PARAM_PHI,
    PARAM_PHI_SPREAD,
    PARAM_STEP,
    PARAM_STEP_RANGE,
    PARAM_P0,
    PARAM_P1,
    PARAM_MAX_SLOPE,
    PARAM_HEIGHT,
    PARAM_SEED,
    PARAM_RANDOMIZE,
    PARAM_ANIMATED,
    PARAM_ACTIVE_PAGE,
    BUTTON_LIKE_CURRENT_IMAGE,

    PARAM_DIMS0
};

enum {
    DENOM_SHIFT = 64
};

typedef enum {
    RIGHT,
    UP,
    LEFT,
    DOWN
} MainDirection;

typedef struct {
    GwyParams *params;
    GwyDataField *field;
    GwyDataField *result;
    /* Cached input image parameters. */
    gdouble zscale;  /* Negative value means there is no input image. */
} ModuleArgs;

typedef struct {
    ModuleArgs *args;
    GtkWidget *dialog;
    GwyParamTable *table_dimensions;
    GwyParamTable *table_generator;
    GwyContainer *data;
    GwyDataField *template_;
} ModuleGUI;

static gboolean         module_register     (void);
static GwyParamDef*     define_module_params(void);
static void             dune_synth          (GwyContainer *data,
                                             GwyRunType runtype);
static gboolean         execute             (ModuleArgs *args,
                                             GtkWindow *wait_window);
static GwyDialogOutcome run_gui             (ModuleArgs *args,
                                             GwyContainer *data,
                                             gint id);
static GtkWidget*       dimensions_tab_new  (ModuleGUI *gui);
static GtkWidget*       generator_tab_new   (ModuleGUI *gui);
static void             param_changed       (ModuleGUI *gui,
                                             gint id);
static void             dialog_response     (ModuleGUI *gui,
                                             gint response);
static void             preview             (gpointer user_data);

static GwyModuleInfo module_info = {
    GWY_MODULE_ABI_VERSION,
    &module_register,
    N_("Generates images by simulating formation of sand dunes."),
    "Yeti <yeti@gwyddion.net>",
    "1.0",
    "David Nečas (Yeti)",
    "2022",
};

GWY_MODULE_QUERY2(module_info, dune_synth)

static gboolean
module_register(void)
{
    gwy_process_func_register("dune_synth",
                              (GwyProcessFunc)&dune_synth,
                              N_("/S_ynthetic/_Deposition/_Dunes..."),
                              GWY_STOCK_SYNTHETIC_DUNES,
                              RUN_MODES,
                              0,
                              N_("Generate image by transporting sand"));

    return TRUE;
}

static GwyParamDef*
define_module_params(void)
{
    static GwyParamDef *paramdef = NULL;

    if (paramdef)
        return paramdef;

    paramdef = gwy_param_def_new();
    gwy_param_def_set_function_name(paramdef, gwy_process_func_current());
    gwy_param_def_add_double(paramdef, PARAM_COVERAGE, "coverage", _("Co_verage"), 0.1, 1000.0, 5.0);
    gwy_param_def_add_int(paramdef, PARAM_NITERS, "niters", _("_Number of iterations"), 1, 100000, 1000);
    gwy_param_def_add_angle(paramdef, PARAM_PHI, "phi", _("_Direction"), FALSE, 1, 0.0);
    gwy_param_def_add_double(paramdef, PARAM_PHI_SPREAD, "phi_spread", _("Spread"), 0.0, 1.0, 0.2);
    gwy_param_def_add_double(paramdef, PARAM_STEP, "step", _("_Minimum step"), 0.25, 100.0, 2.0);
    gwy_param_def_add_double(paramdef, PARAM_STEP_RANGE, "step_range", _("Step range"), 0.0, 100.0, 2.0);
    gwy_param_def_add_double(paramdef, PARAM_P0, "p0", _("Sticking probability on _rock"), 0.0, 1.0, 0.2);
    gwy_param_def_add_double(paramdef, PARAM_P1, "p1", _("Sticking _probability on sand"), 0.0, 1.0, 0.5);
    gwy_param_def_add_int(paramdef, PARAM_MAX_SLOPE, "max_slope", _("Maximum s_lope"), 1, 10, 5);
    gwy_param_def_add_double(paramdef, PARAM_HEIGHT, "height", _("_Height scale"), 1e-4, 1000.0, 1.0);
    gwy_param_def_add_seed(paramdef, PARAM_SEED, "seed", NULL);
    gwy_param_def_add_randomize(paramdef, PARAM_RANDOMIZE, PARAM_SEED, "randomize", NULL, TRUE);
    gwy_param_def_add_boolean(paramdef, PARAM_ANIMATED, "animated", _("Progressive preview"), TRUE);
    gwy_param_def_add_active_page(paramdef, PARAM_ACTIVE_PAGE, "active_page", NULL);

    gwy_synth_define_dimensions_params(paramdef, PARAM_DIMS0);
    return paramdef;
}

static void
dune_synth(GwyContainer *data, GwyRunType runtype)
{
    GwyDialogOutcome outcome = GWY_DIALOG_PROCEED;
    ModuleArgs args;
    GwyDataField *field;
    gint id;

    g_return_if_fail(runtype & RUN_MODES);
    gwy_clear(&args, 1);
    gwy_app_data_browser_get_current(GWY_APP_DATA_FIELD, &field,
                                     GWY_APP_DATA_FIELD_ID, &id,
                                     0);
    args.field = field;
    args.zscale = field ? gwy_data_field_get_rms(field) : -1.0;

    args.params = gwy_params_new_from_settings(define_module_params());
    gwy_synth_sanitise_params(args.params, PARAM_DIMS0, field);
    if (runtype == GWY_RUN_INTERACTIVE) {
        outcome = run_gui(&args, data, id);
        gwy_params_save_to_settings(args.params);
        if (outcome == GWY_DIALOG_CANCEL)
            goto end;
    }
    args.result = gwy_synth_make_result_data_field((args.field = field), args.params, FALSE);
    if (gwy_params_get_boolean(args.params, PARAM_ANIMATED))
        gwy_app_wait_preview_data_field(args.result, data, id);
    if (!execute(&args, gwy_app_find_window_for_channel(data, id)))
        goto end;
    gwy_synth_add_result_to_file(args.result, data, id, args.params);

end:
    GWY_OBJECT_UNREF(args.result);
    g_object_unref(args.params);
}

static GwyDialogOutcome
run_gui(ModuleArgs *args, GwyContainer *data, gint id)
{
    GwyDialogOutcome outcome;
    GwyDialog *dialog;
    GtkWidget *hbox, *dataview;
    GtkNotebook *notebook;
    ModuleGUI gui;

    gwy_clear(&gui, 1);
    gui.args = args;
    gui.template_ = args->field;

    if (gui.template_)
        args->field = gwy_synth_make_preview_data_field(gui.template_, PREVIEW_SIZE);
    else
        args->field = gwy_data_field_new(PREVIEW_SIZE, PREVIEW_SIZE, PREVIEW_SIZE, PREVIEW_SIZE, TRUE);
    args->result = gwy_synth_make_result_data_field(args->field, args->params, TRUE);

    gui.data = gwy_container_new();
    gwy_container_set_object(gui.data, gwy_app_get_data_key_for_id(0), args->result);
    if (gui.template_)
        gwy_app_sync_data_items(data, gui.data, id, 0, FALSE, GWY_DATA_ITEM_GRADIENT, 0);

    gui.dialog = gwy_dialog_new(_("Dunes"));
    dialog = GWY_DIALOG(gui.dialog);
    gwy_dialog_add_buttons(dialog, GWY_RESPONSE_UPDATE, GWY_RESPONSE_RESET, GTK_RESPONSE_CANCEL, GTK_RESPONSE_OK, 0);

    dataview = gwy_create_preview(gui.data, 0, PREVIEW_SIZE, FALSE);
    hbox = gwy_create_dialog_preview_hbox(GTK_DIALOG(dialog), GWY_DATA_VIEW(dataview), FALSE);

    notebook = GTK_NOTEBOOK(gtk_notebook_new());
    gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(notebook), TRUE, TRUE, 0);

    gtk_notebook_append_page(notebook, dimensions_tab_new(&gui), gtk_label_new(_("Dimensions")));
    gtk_notebook_append_page(notebook, generator_tab_new(&gui), gtk_label_new(_("Generator")));
    gwy_param_active_page_link_to_notebook(args->params, PARAM_ACTIVE_PAGE, notebook);

    g_signal_connect_swapped(gui.table_dimensions, "param-changed", G_CALLBACK(param_changed), &gui);
    g_signal_connect_swapped(gui.table_generator, "param-changed", G_CALLBACK(param_changed), &gui);
    g_signal_connect_swapped(dialog, "response", G_CALLBACK(dialog_response), &gui);
    gwy_dialog_set_preview_func(dialog, GWY_PREVIEW_UPON_REQUEST, preview, &gui, NULL);

    outcome = gwy_dialog_run(dialog);

    g_object_unref(gui.data);
    GWY_OBJECT_UNREF(args->field);
    GWY_OBJECT_UNREF(args->result);

    return outcome;
}

static GtkWidget*
dimensions_tab_new(ModuleGUI *gui)
{
    gui->table_dimensions = gwy_param_table_new(gui->args->params);
    gwy_synth_append_dimensions_to_param_table(gui->table_dimensions, 0);
    gwy_dialog_add_param_table(GWY_DIALOG(gui->dialog), gui->table_dimensions);

    return gwy_param_table_widget(gui->table_dimensions);
}

static GtkWidget*
generator_tab_new(ModuleGUI *gui)
{
    GwyParamTable *table;

    table = gui->table_generator = gwy_param_table_new(gui->args->params);

    gwy_param_table_append_header(table, -1, _("Simulation Parameters"));
    gwy_param_table_append_slider(table, PARAM_COVERAGE);
    gwy_param_table_slider_set_mapping(table, PARAM_COVERAGE, GWY_SCALE_MAPPING_LOG);
    gwy_param_table_append_slider(table, PARAM_NITERS);
    gwy_param_table_slider_set_mapping(table, PARAM_NITERS, GWY_SCALE_MAPPING_LOG);
    gwy_param_table_append_separator(table);
    gwy_param_table_append_slider(table, PARAM_PHI);
    gwy_param_table_append_slider(table, PARAM_PHI_SPREAD);
    gwy_param_table_append_separator(table);
    gwy_param_table_append_slider(table, PARAM_STEP);
    gwy_param_table_slider_set_digits(table, PARAM_STEP, 2);
    gwy_param_table_append_slider(table, PARAM_STEP_RANGE);
    gwy_param_table_slider_set_digits(table, PARAM_STEP, 2);
    gwy_param_table_append_separator(table);
    gwy_param_table_append_slider(table, PARAM_P0);
    gwy_param_table_append_slider(table, PARAM_P1);
    gwy_param_table_append_slider(table, PARAM_MAX_SLOPE);

    gwy_param_table_append_header(table, -1, _("Output"));
    gwy_param_table_append_slider(table, PARAM_HEIGHT);
    gwy_param_table_slider_set_mapping(table, PARAM_HEIGHT, GWY_SCALE_MAPPING_LOG);
    if (gui->template_) {
        gwy_param_table_append_button(table, BUTTON_LIKE_CURRENT_IMAGE, -1, GWY_RESPONSE_SYNTH_INIT_Z,
                                      _("_Like Current Image"));
    }

    gwy_param_table_append_header(table, -1, _("Options"));
    gwy_param_table_append_seed(table, PARAM_SEED);
    gwy_param_table_append_checkbox(table, PARAM_RANDOMIZE);
    gwy_param_table_append_separator(table);
    gwy_param_table_append_checkbox(table, PARAM_ANIMATED);

    gwy_dialog_add_param_table(GWY_DIALOG(gui->dialog), table);

    return gwy_param_table_widget(table);
}

static void
param_changed(ModuleGUI *gui, gint id)
{
    GwyParamTable *table = gui->table_generator;

    if (gwy_synth_handle_param_changed(gui->table_dimensions, id))
        id = -1;

    if (id < 0 || id == PARAM_DIMS0 + GWY_DIMS_PARAM_ZUNIT) {
        static const gint zids[] = { PARAM_HEIGHT };

        gwy_synth_update_value_unitstrs(table, zids, G_N_ELEMENTS(zids));
        gwy_synth_update_like_current_button_sensitivity(table, BUTTON_LIKE_CURRENT_IMAGE);
    }
}

static void
dialog_response(ModuleGUI *gui, gint response)
{
    ModuleArgs *args = gui->args;

    if (response == GWY_RESPONSE_SYNTH_INIT_Z) {
        gdouble zscale = args->zscale;
        gint power10z;

        if (zscale > 0.0) {
            gwy_params_get_unit(args->params, PARAM_DIMS0 + GWY_DIMS_PARAM_ZUNIT, &power10z);
            gwy_param_table_set_double(gui->table_generator, PARAM_HEIGHT, zscale/pow10(power10z));
        }
    }
    else if (response == GWY_RESPONSE_SYNTH_TAKE_DIMS) {
        gwy_synth_use_dimensions_template(gui->table_dimensions);
    }
}

static void
preview(gpointer user_data)
{
    ModuleGUI *gui = (ModuleGUI*)user_data;

    if (execute(gui->args, GTK_WINDOW(gui->dialog)))
        gwy_data_field_data_changed(gui->args->result);
}

/* Steepest descent trickle down, starting at the bottom from newly dug hole. */
static inline void
fill_hole_if_too_steep(gint *desert, gint xres, gint yres, gint j, gint i, gint maxslope)
{
    gint ym1block = (yres - 1)*xres;
    gint k = i*xres + j;

    while (TRUE) {
        gint kup = (k >= xres ? k - xres : k + ym1block);
        gint kdown = (k < ym1block ? k + xres : k - ym1block);
        gint kleft = (k % xres ? k-1 : k + xres-1);
        gint kright = ((k + 1) % xres ? k+1 : k - (xres-1));
        gint hup = desert[kup], hdown = desert[kdown], hleft = desert[kleft], hright = desert[kright];
        gint kmax;

        if (hup > hdown) {
            if (hup > hleft)
                kmax = (hup > hright ? kup : kright);
            else
                kmax = (hleft > hright ? kleft : kright);
        }
        else {
            if (hdown > hleft)
                kmax = (hdown > hright ? kdown : kright);
            else
                kmax = (hleft > hright ? kleft : kright);
        }

        if (desert[kmax] - desert[k] <= maxslope)
            return;

        desert[k]++;
        desert[kmax]--;
        k = kmax;
    }
}

/* Steepest descent trickle down, starting at the top from a newly deposited sand slab. */
static inline void
slide_down_if_too_steep(gint *desert, gint xres, gint yres, gint j, gint i, gint maxslope)
{
    gint ym1block = (yres - 1)*xres;
    gint k = i*xres + j;

    while (TRUE) {
        gint kup = (k >= xres ? k - xres : k + ym1block);
        gint kdown = (k < ym1block ? k + xres : k - ym1block);
        gint kleft = (k % xres ? k-1 : k + xres-1);
        gint kright = ((k + 1) % xres ? k+1 : k - (xres-1));
        gint hup = desert[kup], hdown = desert[kdown], hleft = desert[kleft], hright = desert[kright];
        gint kmin;

        if (hup < hdown) {
            if (hup < hleft)
                kmin = (hup < hright ? kup : kright);
            else
                kmin = (hleft < hright ? kleft : kright);
        }
        else {
            if (hdown < hleft)
                kmin = (hdown < hright ? kdown : kright);
            else
                kmin = (hleft < hright ? kleft : kright);
        }

        if (desert[k] - desert[kmin] <= maxslope)
            return;

        desert[k]--;
        desert[kmin]++;
        k = kmin;
    }
}

static void
copy_desert_to_field(GwyDataField *field, const gint *desert)
{
    gdouble *d = gwy_data_field_get_data(field);
    gint i, n = gwy_data_field_get_xres(field)*gwy_data_field_get_yres(field);

    for (i = 0; i < n; i++)
        d[i] = desert[i];
}

static void
compute_local_properties(const gint *desert, gint xres, gint yres, gint j, gint i,
                         gdouble *h, gdouble *gx, gdouble *gy, gdouble *curv)
{
    gint sz_up, sz_down, sz_left, sz_right, sz_cross;
    gint k = i*xres + j;

    if (G_LIKELY(i && j && j+1 < xres && i+1 < yres)) {
        sz_up = desert[k-xres-1] + desert[k-xres] + desert[k-xres+1];
        sz_down = desert[k+xres-1] + desert[k+xres] + desert[k+xres+1];
        sz_left = desert[k-xres-1] + desert[k-1] + desert[k+xres-1];
        sz_right = desert[k-xres+1] + desert[k+1] + desert[k+xres+1];
        sz_cross = desert[k-xres] + desert[k-1] + desert[k+1] + desert[k+xres];
    }
    else {
        gint jplus = (j + 1) % xres, jminus = (j + xres-1) % xres;
        gint iplus = (i + 1) % yres, iminus = (i + yres-1) % yres;
        sz_up = (desert[iminus*xres + jminus] + desert[iminus*xres + j] + desert[iminus*xres + jplus]);
        sz_down = (desert[iplus*xres + jminus] + desert[iplus*xres + j] + desert[iplus*xres + jplus]);
        sz_left = (desert[iminus*xres + jminus] + desert[i*xres + jminus] + desert[iplus*xres + jminus]);
        sz_right = (desert[iminus*xres + jplus] + desert[i*xres + jplus] + desert[iplus*xres + jplus]);
        sz_cross = desert[iminus*xres + j] + desert[i*xres + jminus] + desert[i*xres + jplus] + desert[iplus*xres + j];
    }

    *h = desert[k];
    *gx = (sz_right - sz_left)/6.0;
    *gy = (sz_down - sz_up)/6.0;
    *curv = 0.25*sz_cross - desert[k];
}

static guint
take_rand_int(GwyRandGenSet *rngset, gint *navail, gint size, guint *random_values)
{
    if (G_LIKELY(*navail))
        return random_values[size - (*navail)--];
    gwy_rand_gen_set_fill_ints(rngset, random_values, size);
    *navail = size-1;
    return random_values[0];
}

static guint
take_rand_int_range(GwyRandGenSet *rngset, gint *navail, gint size, guint *random_values,
                    guint minvalue, guint maxvalue)
{
    guint r, R = maxvalue+1 - minvalue;

    g_return_val_if_fail(maxvalue >= minvalue, minvalue);
    if (G_UNLIKELY(minvalue == 0 && maxvalue == G_MAXUINT32))
        return take_rand_int(rngset, navail, size, random_values);

    do {
        r = take_rand_int(rngset, navail, size, random_values);
    } while (r/R == G_MAXUINT32/R);

    return minvalue + r % R;
}

static gdouble
take_rand_double(GwyRandGenSet *rngset, gint *navail, gint size, gdouble *random_values)
{
    if (G_LIKELY(*navail))
        return random_values[size - (*navail)--];
    gwy_rand_gen_set_fill_doubles(rngset, random_values, size);
    *navail = size-1;
    return random_values[0];
}

static gboolean
execute(ModuleArgs *args, GtkWindow *wait_window)
{
    GwyParams *params = args->params;
    gboolean do_initialise = gwy_params_get_boolean(params, PARAM_DIMS0 + GWY_DIMS_PARAM_INITIALIZE);
    gdouble height = gwy_params_get_double(params, PARAM_HEIGHT);
    gboolean animated = gwy_params_get_boolean(params, PARAM_ANIMATED);
    gsize niters_param = gwy_params_get_int(params, PARAM_NITERS);
    gdouble coverage = gwy_params_get_double(params, PARAM_COVERAGE);
    gdouble p0 = gwy_params_get_double(params, PARAM_P0);
    gdouble p1 = gwy_params_get_double(params, PARAM_P1);
    gdouble phi_mean = gwy_params_get_double(params, PARAM_PHI);
    gdouble phi_spread = gwy_params_get_double(params, PARAM_PHI_SPREAD);
    gdouble step = gwy_params_get_double(params, PARAM_STEP);
    gdouble step_range = gwy_params_get_double(params, PARAM_STEP_RANGE);
    gint maxslope = gwy_params_get_int(params, PARAM_MAX_SLOPE);
    GwyRandGenSet *rngset;
    GwySynthUpdateType update;
    GwyDataField *field = args->field, *result = args->result;
    gsize iter, niters;
    gint power10z, xres, yres, n, k, i, j, rngbufsize;
    guint r, p0_int, p1_int;
    gint n_ints = 0, n_doubles = 0;
    gint *desert;
    const gdouble *f;
    gdouble preview_time = (animated ? 1.25 : 0.0);
    guint32 *random_ints;
    gdouble *random_doubles;
    GTimer *timer;
    gboolean finished = FALSE;

    gwy_app_wait_start(wait_window, _("Initializing..."));

    gwy_params_get_unit(params, PARAM_DIMS0 + GWY_DIMS_PARAM_ZUNIT, &power10z);
    height *= pow10(power10z);

    rngset = gwy_rand_gen_set_new(NRANDOM_GENERATORS);
    gwy_rand_gen_set_init(rngset, gwy_params_get_int(params, PARAM_SEED));

    xres = gwy_data_field_get_xres(result);
    yres = gwy_data_field_get_yres(result);
    n = xres*yres;
    /* One iteration given by user is one complete pass through the image, but we count individual sand chunk
     * movements. */
    niters = niters_param*n;
    p0_int = (guint)floor(G_MAXUINT32*p0 + 0.5);
    p1_int = (guint)floor(G_MAXUINT32*p1 + 0.5);
    rngbufsize = 4*n;

    desert = g_new0(gint, n);
    random_ints = g_new(guint32, rngbufsize);
    random_doubles = g_new(gdouble, rngbufsize);
    if (field && do_initialise) {
        gdouble min, avg, q;

        min = gwy_data_field_get_min(field);
        avg = gwy_data_field_get_avg(field);
        if (avg > min) {
            q = coverage/(avg - min);
            f = gwy_data_field_get_data_const(field);
            for (k = 0; k < n; k++)
                desert[k] = GWY_ROUND(q*(f[k] - min));
        }
    }
    else {
        gint h0 = floor(coverage);
        gdouble hplus = coverage - h0;

        gwy_rand_gen_set_fill_doubles(rngset, random_doubles, n);
        for (k = 0; k < n; k++)
            desert[k] = h0 + (random_doubles[k] < hplus);
    }

    timer = g_timer_new();
    gwy_synth_update_progress(NULL, 0, 0, 0);
    if (!gwy_app_wait_set_message(_("Running computation...")))
        goto end;

    for (iter = 0; iter < niters; iter++) {
        gint x, y, dx, dy, scaled_xres = xres*DENOM_SHIFT, scaled_yres = yres*DENOM_SHIFT;
        gdouble vx, vy, h, gx, gy, curv, w, p_height, p_slope, p_curv;

        k = take_rand_int_range(rngset, &n_ints, rngbufsize, random_ints, 0, n-1);
        if (!desert[k])
            continue;

        /* XXX: There may be various direction generation methods, even dynamically changing ones. */
        if (TRUE) {
            gdouble phi = take_rand_double(rngset, &n_doubles, rngbufsize, random_doubles);
            gdouble v = step + step_range*take_rand_double(rngset, &n_doubles, rngbufsize, random_doubles);

            phi = phi_mean + 2.0*G_PI*phi_spread*(phi - 0.5);
            vx = v*cos(phi);
            vy = v*sin(phi);
        }

        /* This is generic once we have vx and vy. */
        dx = GWY_ROUND(vx*DENOM_SHIFT);
        dy = GWY_ROUND(vy*DENOM_SHIFT);
        i = k/xres;
        j = k % xres;

        /* Compute pick-up probability using
         * (a) possibly height, with higher elevation (h) corresponding to higher pick-up probability, even though
         *     there is no clear overall normalisation; we can use probability like
         *     p(h) = h/(h + c₀)
         * (b) possibly slope in the wind direction, with higher upslope (g) meaning lower pick-up probability; again,
         *     we can use an increasing function like a sigmoid
         *     p(g) = 1/2*(1 - g/√(g² + c₁))
         * (c) possibly curvature, with smaller curvature (k) meaning higher pick-up probability (i.e. negative
         *     curvature means higher probability);
         *     p(k) = 1/2*(1 - k/√(k² + c₂))
         */
        compute_local_properties(desert, xres, yres, j, i, &h, &gx, &gy, &curv);
        p_height = h/(1.0 + h);
        w = (gx*dx + gy*dy)/sqrt(dx*dx + dy*dy);
        p_slope = 0.5*(1 + w/sqrt(w*w + 1.0));
        p_curv = 0.5*(1 - curv/sqrt(curv*curv + 5.0));
        if (take_rand_double(rngset, &n_doubles, rngbufsize, random_doubles) > p_height*p_slope*p_curv)
            continue;

        x = take_rand_int_range(rngset, &n_ints, rngbufsize, random_ints, 0, DENOM_SHIFT*DENOM_SHIFT-1);
        y = i*DENOM_SHIFT + x/DENOM_SHIFT;
        x = j*DENOM_SHIFT + x % DENOM_SHIFT;

        desert[k]--;
        fill_hole_if_too_steep(desert, xres, yres, j, i, maxslope);
        while (TRUE) {
            x = (x + scaled_xres + dx) % scaled_xres;
            y = (y + scaled_yres + dy) % scaled_yres;
            j = x/DENOM_SHIFT;
            i = y/DENOM_SHIFT;
            k = i*xres + j;
            r = take_rand_int(rngset, &n_ints, rngbufsize, random_ints);
            if ((desert[k] && r < p1_int) || r < p0_int)
                break;
        }
        desert[k]++;
        slide_down_if_too_steep(desert, xres, yres, j, i, maxslope);

        if (iter % 100000 == 0) {
            update = gwy_synth_update_progress(timer, preview_time, iter, niters);
            if (update == GWY_SYNTH_UPDATE_CANCELLED)
                return FALSE;
            if (update == GWY_SYNTH_UPDATE_DO_PREVIEW) {
                copy_desert_to_field(result, desert);
                gwy_data_field_data_changed(result);
            }
        }
    }

    copy_desert_to_field(result, desert);
    gwy_data_field_multiply(result, height);
    finished = TRUE;

end:
    gwy_app_wait_finish();
    gwy_rand_gen_set_free(rngset);
    g_timer_destroy(timer);
    g_free(desert);
    g_free(random_ints);
    g_free(random_doubles);

    return finished;
}

/* vim: set cin columns=120 tw=118 et ts=4 sw=4 cino=>1s,e0,n0,f0,{0,}0,^0,\:1s,=0,g1s,h0,t0,+1s,c3,(0,u0 : */