← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to raster/r.li/r.li.mps/mps.c

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/r.li/r.li.mps/mps.c
1
 
/*
2
 
 * \brief calculates mean patch size index
3
 
 *
4
 
 *  \AUTHOR: Serena Pallecchi student of Computer Science University of Pisa (Italy)
5
 
 *                      Commission from Faunalia Pontedera (PI) www.faunalia.it
6
 
 *
7
 
 *   This program is free software under the GPL (>=v2)
8
 
 *   Read the COPYING file that comes with GRASS for details.
9
 
 *       
10
 
 *       \BUGS: please send bugs reports to pallecch@cli.di.unipi.it
11
 
 */
12
 
 
13
 
 
14
 
#include <grass/gis.h>
15
 
#include <grass/glocale.h>
 
1
/****************************************************************************
 
2
 *
 
3
 * MODULE:       r.li.patchnum
 
4
 * AUTHOR(S):    Serena Pallecchi (original contributor)
 
5
 *                student of Computer Science University of Pisa (Italy)
 
6
 *               Commission from Faunalia Pontedera (PI) www.faunalia.it
 
7
 *               Rewrite: Markus Metz
 
8
 *               Patch identification: Michael Shapiro - CERL
 
9
 *
 
10
 * PURPOSE:      calculates mean patch size index
 
11
 * COPYRIGHT:    (C) 2007-2014 by the GRASS Development Team
 
12
 *
 
13
 *               This program is free software under the GNU General Public
 
14
 *               License (>=v2). Read the file COPYING that comes with GRASS
 
15
 *               for details.
 
16
 *
 
17
 *****************************************************************************/
16
18
 
17
19
#include <stdlib.h>
18
20
#include <fcntl.h>
19
21
#include <math.h>
20
22
 
21
 
#include "../r.li.daemon/defs.h"
22
 
#include "../r.li.daemon/avlDefs.h"
23
 
#include "../r.li.daemon/avlID.h"
 
23
#include <grass/gis.h>
 
24
#include <grass/raster.h>
 
25
#include <grass/glocale.h>
 
26
 
24
27
#include "../r.li.daemon/daemon.h"
25
 
 
26
 
 
27
 
int calculate(int fd, area_des ad, struct Cell_head hd, double *result);
28
 
int calculateD(int fd, area_des ad, struct Cell_head hd, double *result);
29
 
int calculateF(int fd, area_des ad, struct Cell_head hd, double *result);
30
 
 
 
28
#include "../r.li.daemon/GenericCell.h"
 
29
 
 
30
/* template is patchnum */
 
31
 
 
32
/* cell count and type of each patch */
 
33
struct pst {
 
34
    long count;
 
35
    generic_cell type;
 
36
};
 
37
 
 
38
rli_func meanPatchSize;
 
39
int calculate(int fd, struct area_entry *ad, double *result);
 
40
int calculateD(int fd, struct area_entry *ad, double *result);
 
41
int calculateF(int fd, struct area_entry *ad, double *result);
31
42
 
32
43
int main(int argc, char *argv[])
33
44
{
38
49
    module = G_define_module();
39
50
    module->description =
40
51
        _("Calculates mean patch size index on a raster map, using a 4 neighbour algorithm");
41
 
    module->keywords = _("raster, landscape structure analysis, patch index");
 
52
    G_add_keyword(_("raster"));
 
53
    G_add_keyword(_("landscape structure analysis"));
 
54
    G_add_keyword(_("patch index"));
 
55
 
42
56
    /* define options */
43
57
 
44
 
    raster = G_define_standard_option(G_OPT_R_MAP);
 
58
    raster = G_define_standard_option(G_OPT_R_INPUT);
45
59
 
46
 
    conf = G_define_option();
47
 
    conf->key = "conf";
 
60
    conf = G_define_standard_option(G_OPT_F_INPUT);
 
61
    conf->key = "config";
48
62
    conf->description = _("Configuration file");
49
 
    conf->type = TYPE_STRING;
50
63
    conf->required = YES;
51
 
    conf->gisprompt = "old_file,file,input";
52
64
 
53
65
    output = G_define_standard_option(G_OPT_R_OUTPUT);
54
66
 
57
69
 
58
70
    return calculateIndex(conf->answer, meanPatchSize, NULL, raster->answer,
59
71
                          output->answer);
60
 
 
61
72
}
62
73
 
63
74
 
64
 
 
65
 
int meanPatchSize(int fd, char **par, area_des ad, double *result)
 
75
int meanPatchSize(int fd, char **par, struct area_entry *ad, double *result)
66
76
{
67
 
 
68
 
    char *mapset;
69
 
 
70
 
    int ris = 0;
71
 
 
 
77
    int ris = RLI_OK;
72
78
    double indice = 0;
73
79
 
74
 
    struct Cell_head hd;
75
 
 
76
 
 
77
 
 
78
 
    mapset = G_find_cell(ad->raster, "");
79
 
    if (G_get_cellhd(ad->raster, mapset, &hd) == -1)
80
 
        return RLI_ERRORE;
81
 
 
82
 
 
83
 
 
84
80
    switch (ad->data_type) {
85
81
    case CELL_TYPE:
86
82
        {
87
 
            ris = calculate(fd, ad, hd, &indice);
 
83
            ris = calculate(fd, ad, &indice);
88
84
            break;
89
85
        }
90
86
    case DCELL_TYPE:
91
87
        {
92
 
            ris = calculateD(fd, ad, hd, &indice);
 
88
            ris = calculateD(fd, ad, &indice);
93
89
            break;
94
90
        }
95
91
    case FCELL_TYPE:
96
92
        {
97
 
            ris = calculateF(fd, ad, hd, &indice);
 
93
            ris = calculateF(fd, ad, &indice);
98
94
            break;
99
95
        }
100
96
    default:
103
99
            return RLI_ERRORE;
104
100
        }
105
101
    }
 
102
 
106
103
    if (ris != RLI_OK) {
107
104
        return RLI_ERRORE;
108
105
    }
109
106
 
110
107
    *result = indice;
111
 
    return RLI_OK;
112
 
 
113
 
}
114
 
 
115
 
 
116
 
int calculate(int fd, area_des ad, struct Cell_head hd, double *result)
117
 
{
118
 
    CELL *buf;
119
 
    CELL *buf_sup;
120
 
 
121
 
    CELL corrCell;
122
 
    CELL precCell;
123
 
    CELL supCell;
124
 
 
125
 
    int i, j;
126
 
    int mask_fd = -1, *mask_buf;
127
 
    int ris = 0;
128
 
    int masked = FALSE;
129
 
 
130
 
 
131
 
    long npatch = 0;
132
 
    long tot = 0;
133
 
    long zero = 0;
134
 
    long uno = 1;
135
 
    long idCorr = 0;
136
 
    long lastId = 0;
137
 
    long doppi = 0;
138
 
    long *mask_patch_sup;
139
 
    long *mask_patch_corr;
140
 
 
141
 
    double indice = 0;
142
 
    double area = 0;            /*if all cells are null area=0 */
143
 
    double areaCorrect = 0;
144
 
    double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
145
 
 
146
 
    avlID_tree albero = NULL;
147
 
 
148
 
    avlID_table *array;
149
 
 
150
 
 
151
 
 
152
 
    /* open mask if needed */
153
 
    if (ad->mask == 1) {
154
 
        if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
155
 
            return RLI_ERRORE;
156
 
        mask_buf = G_malloc(ad->cl * sizeof(int));
157
 
        if (mask_buf == NULL) {
158
 
            G_fatal_error("malloc mask_buf failed");
159
 
            return RLI_ERRORE;
160
 
        }
161
 
        masked = TRUE;
162
 
    }
163
 
 
164
 
    mask_patch_sup = G_malloc(ad->cl * sizeof(long));
165
 
    if (mask_patch_sup == NULL) {
166
 
        G_fatal_error("malloc mask_patch_sup failed");
167
 
        return RLI_ERRORE;
168
 
    }
169
 
 
170
 
    mask_patch_corr = G_malloc(ad->cl * sizeof(long));
171
 
    if (mask_patch_corr == NULL) {
172
 
        G_fatal_error("malloc mask_patch_corr failed");
173
 
        return RLI_ERRORE;
174
 
    }
175
 
 
176
 
    buf_sup = G_allocate_cell_buf();
177
 
    if (buf_sup == NULL) {
178
 
        G_fatal_error("malloc buf_sup failed");
179
 
        return RLI_ERRORE;
180
 
    }
181
 
 
182
 
    buf = G_allocate_cell_buf();
183
 
    if (buf == NULL) {
184
 
        G_fatal_error("malloc buf failed");
185
 
        return RLI_ERRORE;
186
 
    }
187
 
 
188
 
    G_set_c_null_value(buf_sup + ad->x, ad->cl);        /*the first time buf_sup is all null */
189
 
 
190
 
    for (i = 0; i < ad->cl; i++) {
191
 
        mask_patch_sup[i] = 0;
192
 
        mask_patch_corr[i] = 0;
193
 
    }
194
 
 
195
 
    /*for each raster row */
196
 
 
197
 
    for (j = 0; j < ad->rl; j++) {
198
 
        if (j > 0) {
199
 
            buf_sup = RLI_get_cell_raster_row(fd, j - 1 + ad->y, ad);
200
 
        }
201
 
        buf = RLI_get_cell_raster_row(fd, j + ad->y, ad);
202
 
        if (masked) {
203
 
            if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0) {
204
 
                G_fatal_error("mask read failed");
205
 
                return RLI_ERRORE;
206
 
            }
207
 
        }
208
 
        G_set_c_null_value(&precCell, 1);
209
 
        for (i = 0; i < ad->cl; i++) {  /*for each cell in the row */
210
 
 
211
 
            corrCell = buf[i + ad->x];
212
 
 
213
 
            if ((masked) && (mask_buf[i + ad->x] == 0)) {
214
 
                G_set_c_null_value(&corrCell, 1);
215
 
            }
216
 
 
217
 
            /*valid cell */
218
 
            if (!(G_is_null_value(&corrCell, CELL_TYPE))) {
219
 
                area++;
220
 
                if (i > 0)
221
 
                    precCell = buf[i - 1 + ad->x];
222
 
 
223
 
                if (j == 0)
224
 
                    G_set_c_null_value(&supCell, 1);
225
 
                else
226
 
                    supCell = buf_sup[i + ad->x];
227
 
 
228
 
 
229
 
                if (corrCell != precCell) {
230
 
                    if (corrCell != supCell) {
231
 
                        /*new patch */
232
 
                        if (idCorr == 0) {      /*first patch */
233
 
                            lastId = 1;
234
 
                            idCorr = 1;
235
 
                            mask_patch_corr[i] = idCorr;
236
 
                        }
237
 
                        else {  /*not first patch */
238
 
                            /* put in the tree previous values */
239
 
                            if (albero == NULL) {
240
 
                                albero = avlID_make(idCorr, uno);
241
 
                                if (albero == NULL) {
242
 
                                    G_fatal_error("avlID_make error");
243
 
                                    return RLI_ERRORE;
244
 
                                }
245
 
                                npatch++;
246
 
 
247
 
                            }
248
 
                            else {      /* tree not null */
249
 
 
250
 
                                ris = avlID_add(&albero, idCorr, uno);
251
 
                                switch (ris) {
252
 
                                case AVL_ERR:
253
 
                                    {
254
 
                                        G_fatal_error("avlID_add error");
255
 
                                        return RLI_ERRORE;
256
 
                                    }
257
 
                                case AVL_ADD:
258
 
                                    {
259
 
                                        npatch++;
260
 
                                        break;
261
 
                                    }
262
 
                                case AVL_PRES:
263
 
                                    {
264
 
                                        break;
265
 
                                    }
266
 
                                default:
267
 
                                    {
268
 
                                        G_fatal_error
269
 
                                            ("avlID_add unknown error");
270
 
                                        return RLI_ERRORE;
271
 
                                    }
272
 
                                }
273
 
                            }
274
 
 
275
 
                            lastId++;
276
 
                            idCorr = lastId;
277
 
                            mask_patch_corr[i] = idCorr;
278
 
                        }
279
 
                    }
280
 
                    else {      /*current cell and upper cell are equal */
281
 
 
282
 
                        if ((corrCell == precCell) &&
283
 
                            (mask_patch_sup[i] != mask_patch_corr[i - 1])) {
284
 
                            long r = 0;
285
 
                            long del = mask_patch_sup[i];
286
 
 
287
 
                            r = avlID_sub(&albero, del);
288
 
                            if (r == 0) {
289
 
                                G_fatal_error("avlID_sub error");
290
 
                                return RLI_ERRORE;
291
 
                            }
292
 
                            /*Remove one patch because it makes part of a patch already found */
293
 
                            ris = avlID_add(&albero, idCorr, uno);
294
 
                            switch (ris) {
295
 
                            case AVL_ERR:
296
 
                                {
297
 
                                    G_fatal_error("avlID_add error");
298
 
                                    return RLI_ERRORE;
299
 
                                }
300
 
                            case AVL_ADD:
301
 
                                {
302
 
                                    npatch++;
303
 
                                    break;
304
 
                                }
305
 
                            case AVL_PRES:
306
 
                                {
307
 
                                    break;
308
 
                                }
309
 
                            default:
310
 
                                {
311
 
                                    G_fatal_error("avlID_add unknown error");
312
 
                                    return RLI_ERRORE;
313
 
                                }
314
 
                            }
315
 
                            r = i;
316
 
                            while (i < ad->cl) {
317
 
                                if (mask_patch_sup[r] == del) {
318
 
                                    mask_patch_sup[r] = idCorr;
319
 
                                }
320
 
                                else {
321
 
                                    r = ad->cl + 1;
322
 
                                }
323
 
                            }
324
 
                        }
325
 
 
326
 
                        if (albero == NULL) {
327
 
                            albero = avlID_make(idCorr, uno);
328
 
                            if (albero == NULL) {
329
 
                                G_fatal_error("avlID_make error");
330
 
                                return RLI_ERRORE;
331
 
                            }
332
 
                            npatch++;
333
 
                        }
334
 
                        else {  /* tree not null */
335
 
 
336
 
                            ris = avlID_add(&albero, idCorr, uno);
337
 
                            switch (ris) {
338
 
                            case AVL_ERR:
339
 
                                {
340
 
                                    G_fatal_error("avlID_add error");
341
 
                                    return RLI_ERRORE;
342
 
                                }
343
 
                            case AVL_ADD:
344
 
                                {
345
 
                                    npatch++;
346
 
                                    break;
347
 
                                }
348
 
                            case AVL_PRES:
349
 
                                {
350
 
                                    break;
351
 
                                }
352
 
                            default:
353
 
                                {
354
 
                                    G_fatal_error("avlID_add unknown error");
355
 
                                    return RLI_ERRORE;
356
 
                                }
357
 
                            }
358
 
                        }
359
 
 
360
 
                        idCorr = mask_patch_sup[i];
361
 
                        mask_patch_corr[i] = idCorr;
362
 
                    }
363
 
                }
364
 
                else {          /*current cell and previous cell are equal */
365
 
 
366
 
 
367
 
                    if ((corrCell == supCell) &&
368
 
                        (mask_patch_sup[i] != mask_patch_corr[i - 1])) {
369
 
                        int l;
370
 
 
371
 
                        mask_patch_corr[i] = mask_patch_sup[i];
372
 
                        l = i - 1;
373
 
                        while (l >= 0) {
374
 
                            if (mask_patch_corr[l] == idCorr) {
375
 
                                mask_patch_corr[l] = mask_patch_sup[i];
376
 
                                l--;
377
 
                            }
378
 
                            else {
379
 
                                l = (-1);
380
 
                            }
381
 
                        }
382
 
                        lastId--;
383
 
                        idCorr = mask_patch_sup[i];
384
 
                    }
385
 
                    else {
386
 
                        mask_patch_corr[i] = idCorr;
387
 
                    }
388
 
 
389
 
                }
390
 
            }
391
 
            else {              /*cell not to consider or cell is null */
392
 
 
393
 
                mask_patch_corr[i] = 0;
394
 
            }
395
 
        }
396
 
        mask_patch_sup = mask_patch_corr;
397
 
    }
398
 
 
399
 
 
400
 
 
401
 
    if (area != 0) {
402
 
        if (albero == NULL) {
403
 
            albero = avlID_make(idCorr, uno);
404
 
            if (albero == NULL) {
405
 
                G_fatal_error("avlID_make error");
406
 
                return RLI_ERRORE;
407
 
            }
408
 
            npatch++;
409
 
        }
410
 
        else {
411
 
            ris = avlID_add(&albero, idCorr, uno);
412
 
            switch (ris) {
413
 
            case AVL_ERR:
414
 
                {
415
 
                    G_fatal_error("avlID_add error");
416
 
                    return RLI_ERRORE;
417
 
                }
418
 
            case AVL_ADD:
419
 
                {
420
 
                    npatch++;
421
 
                    break;
422
 
                }
423
 
            case AVL_PRES:
424
 
                {
425
 
                    break;
426
 
                }
427
 
            default:
428
 
                {
429
 
                    G_fatal_error("avlID_add unknown error");
430
 
                    return RLI_ERRORE;
431
 
                }
432
 
            }
433
 
        }
434
 
 
435
 
 
436
 
        array = G_malloc(npatch * sizeof(avlID_tableRow));
437
 
        if (array == NULL) {
438
 
            G_fatal_error("malloc array failed");
439
 
            return RLI_ERRORE;
440
 
        }
441
 
        tot = avlID_to_array(albero, zero, array);
442
 
 
443
 
        if (tot != npatch) {
444
 
            G_warning
445
 
                ("avlID_to_array unaspected value. the result could be wrong");
446
 
            return RLI_ERRORE;
447
 
        }
448
 
 
449
 
        for (i = 0; i < npatch; i++) {
450
 
            if (array[i]->tot == 0)
451
 
                doppi++;
452
 
        }
453
 
        npatch = npatch - doppi;
454
 
 
455
 
        /*calculate distance */
456
 
        G_begin_distance_calculations();
457
 
        /* EW Dist at North edge */
458
 
        EW_DIST1 = G_distance(hd.east, hd.north, hd.west, hd.north);
459
 
        /* EW Dist at South Edge */
460
 
        EW_DIST2 = G_distance(hd.east, hd.south, hd.west, hd.south);
461
 
        /* NS Dist at East edge */
462
 
        NS_DIST1 = G_distance(hd.east, hd.north, hd.east, hd.south);
463
 
        /* NS Dist at West edge */
464
 
        NS_DIST2 = G_distance(hd.west, hd.north, hd.west, hd.south);
465
 
 
466
 
        areaCorrect = (((EW_DIST1 + EW_DIST2) / 2) / hd.cols) *
467
 
            (((NS_DIST1 + NS_DIST2) / 2) / hd.rows) * (area);
468
 
        indice = areaCorrect / npatch;
469
 
        G_free(array);
470
 
    }
471
 
    else
472
 
        indice = (double)(0);
473
 
 
474
 
 
475
 
    *result = indice;
476
 
 
477
 
 
478
 
    if (masked)
479
 
        G_free(mask_buf);
480
 
 
481
 
    G_free(mask_patch_corr);
482
 
 
483
 
    G_free(buf_sup);
484
 
    return RLI_OK;
485
 
}
486
 
 
487
 
 
488
 
int calculateD(int fd, area_des ad, struct Cell_head hd, double *result)
489
 
{
490
 
    DCELL *buf;
491
 
    DCELL *buf_sup;
492
 
 
493
 
    DCELL corrCell;
494
 
    DCELL precCell;
495
 
    DCELL supCell;
496
 
 
497
 
    int i, j;
498
 
    int mask_fd = -1, *mask_buf;
499
 
    int ris = 0;
500
 
    int masked = FALSE;
501
 
 
502
 
    long npatch = 0;
503
 
    long tot = 0;
504
 
    long zero = 0;
505
 
    long uno = 1;
506
 
    long idCorr = 0;
507
 
    long lastId = 0;
508
 
    long doppi = 0;
509
 
    long *mask_patch_sup;
510
 
    long *mask_patch_corr;
511
 
 
512
 
    double indice = 0;
513
 
    double area = 0;            /*if all cells are null area=0 */
514
 
    double areaCorrect = 0;
515
 
    double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
516
 
 
517
 
    avlID_tree albero = NULL;
518
 
 
519
 
    avlID_table *array;
520
 
 
521
 
 
522
 
 
523
 
    /* open mask if needed */
524
 
    if (ad->mask == 1) {
525
 
        if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
526
 
            return RLI_ERRORE;
527
 
        mask_buf = G_malloc(ad->cl * sizeof(int));
528
 
        if (mask_buf == NULL) {
529
 
            G_fatal_error("malloc mask_buf failed");
530
 
            return RLI_ERRORE;
531
 
        }
532
 
        masked = TRUE;
533
 
    }
534
 
 
535
 
    mask_patch_sup = G_malloc(ad->cl * sizeof(long));
536
 
    if (mask_patch_sup == NULL) {
537
 
        G_fatal_error("malloc mask_patch_sup failed");
538
 
        return RLI_ERRORE;
539
 
    }
540
 
 
541
 
    mask_patch_corr = G_malloc(ad->cl * sizeof(long));
542
 
    if (mask_patch_corr == NULL) {
543
 
        G_fatal_error("malloc mask_patch_corr failed");
544
 
        return RLI_ERRORE;
545
 
    }
546
 
 
547
 
    buf_sup = G_allocate_d_raster_buf();
548
 
    if (buf_sup == NULL) {
549
 
        G_fatal_error("malloc buf_sup failed");
550
 
        return RLI_ERRORE;
551
 
    }
552
 
 
553
 
    buf = G_allocate_d_raster_buf();
554
 
    if (buf == NULL) {
555
 
        G_fatal_error("malloc buf failed");
556
 
        return RLI_ERRORE;
557
 
    }
558
 
 
559
 
    G_set_d_null_value(buf_sup + ad->x, ad->cl);        /*the first time buf_sup is all null */
560
 
 
561
 
    for (i = 0; i < ad->cl; i++) {
562
 
        mask_patch_sup[i] = 0;
563
 
        mask_patch_corr[i] = 0;
564
 
    }
565
 
 
566
 
    /*read each row and put in an avlId tree cell value with the number of cells which have that value */
567
 
    for (j = 0; j < ad->rl; j++) {
568
 
        if (j > 0) {
569
 
            buf_sup = RLI_get_dcell_raster_row(fd, j - 1 + ad->y, ad);
570
 
        }
571
 
        buf = RLI_get_dcell_raster_row(fd, j + ad->y, ad);
572
 
        if (masked) {
573
 
            if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0) {
574
 
                G_fatal_error("mask read failed");
575
 
                return RLI_ERRORE;
576
 
            }
577
 
        }
578
 
        G_set_d_null_value(&precCell, 1);
579
 
        for (i = 0; i < ad->cl; i++) {  /*for each cell in the row */
580
 
 
581
 
            corrCell = buf[i + ad->x];
582
 
 
583
 
            if ((masked) && (mask_buf[i + ad->x] == 0)) {
584
 
                G_set_d_null_value(&corrCell, 1);
585
 
            }
586
 
 
587
 
            if (!(G_is_null_value(&corrCell, DCELL_TYPE))) {
588
 
                area++;
589
 
                if (i > 0)
590
 
                    precCell = buf[i - 1 + ad->x];
591
 
 
592
 
                if (j == 0)
593
 
                    G_set_d_null_value(&supCell, 1);
594
 
                else
595
 
                    supCell = buf_sup[i + ad->x];
596
 
 
597
 
                if (corrCell != precCell) {
598
 
                    if (corrCell != supCell) {
599
 
                        /*new patch */
600
 
                        if (idCorr == 0) {      /*first patch */
601
 
                            lastId = 1;
602
 
                            idCorr = 1;
603
 
                            mask_patch_corr[i] = idCorr;
604
 
                        }
605
 
                        else {  /*not first patch */
606
 
                            /* put in the tree previous value */
607
 
                            if (albero == NULL) {
608
 
                                albero = avlID_make(idCorr, uno);
609
 
                                if (albero == NULL) {
610
 
                                    G_fatal_error("avlID_make error");
611
 
                                    return RLI_ERRORE;
612
 
                                }
613
 
                                npatch++;
614
 
 
615
 
                            }
616
 
                            else {      /* tree not null */
617
 
 
618
 
                                ris = avlID_add(&albero, idCorr, uno);
619
 
                                switch (ris) {
620
 
                                case AVL_ERR:
621
 
                                    {
622
 
                                        G_fatal_error("avlID_add error");
623
 
                                        return RLI_ERRORE;
624
 
                                    }
625
 
                                case AVL_ADD:
626
 
                                    {
627
 
                                        npatch++;
628
 
                                        break;
629
 
                                    }
630
 
                                case AVL_PRES:
631
 
                                    {
632
 
                                        break;
633
 
                                    }
634
 
                                default:
635
 
                                    {
636
 
                                        G_fatal_error
637
 
                                            ("avlID_add unknown error");
638
 
                                        return RLI_ERRORE;
639
 
                                    }
640
 
                                }
641
 
                            }
642
 
 
643
 
                            lastId++;
644
 
                            idCorr = lastId;
645
 
                            mask_patch_corr[i] = idCorr;
646
 
                        }
647
 
                    }
648
 
                    else {      /*current cell and upper cell are equal */
649
 
 
650
 
                        if ((corrCell == precCell) &&
651
 
                            (mask_patch_sup[i] != mask_patch_corr[i - 1])) {
652
 
                            long r = 0;
653
 
                            long del = mask_patch_sup[i];
654
 
 
655
 
                            /*Remove one patch because it makes part of a patch already found */
656
 
                            r = avlID_sub(&albero, del);
657
 
                            if (r == 0) {
658
 
                                G_fatal_error("avlID_sub error");
659
 
                                return RLI_ERRORE;
660
 
                            }
661
 
 
662
 
                            ris = avlID_add(&albero, idCorr, uno);
663
 
                            switch (ris) {
664
 
                            case AVL_ERR:
665
 
                                {
666
 
                                    G_fatal_error("avlID_add error");
667
 
                                    return RLI_ERRORE;
668
 
                                }
669
 
                            case AVL_ADD:
670
 
                                {
671
 
                                    npatch++;
672
 
                                    break;
673
 
                                }
674
 
                            case AVL_PRES:
675
 
                                {
676
 
                                    break;
677
 
                                }
678
 
                            default:
679
 
                                {
680
 
                                    G_fatal_error("avlID_add unknown error");
681
 
                                    return RLI_ERRORE;
682
 
                                }
683
 
                            }
684
 
                            r = i;
685
 
                            while (i < ad->cl) {
686
 
                                if (mask_patch_sup[r] == del) {
687
 
                                    mask_patch_sup[r] = idCorr;
688
 
                                }
689
 
                                else {
690
 
                                    r = ad->cl + 1;
691
 
                                }
692
 
                            }
693
 
                        }
694
 
 
695
 
                        if (albero == NULL) {
696
 
                            albero = avlID_make(idCorr, uno);
697
 
                            if (albero == NULL) {
698
 
                                G_fatal_error("avlID_make error");
699
 
                                return RLI_ERRORE;
700
 
                            }
701
 
                            npatch++;
702
 
                        }
703
 
                        else {  /*tree not null */
704
 
 
705
 
                            ris = avlID_add(&albero, idCorr, uno);
706
 
                            switch (ris) {
707
 
                            case AVL_ERR:
708
 
                                {
709
 
                                    G_fatal_error("avlID_add error");
710
 
                                    return RLI_ERRORE;
711
 
                                }
712
 
                            case AVL_ADD:
713
 
                                {
714
 
                                    npatch++;
715
 
                                    break;
716
 
                                }
717
 
                            case AVL_PRES:
718
 
                                {
719
 
                                    break;
720
 
                                }
721
 
                            default:
722
 
                                {
723
 
                                    G_fatal_error("avlID_add unknown error");
724
 
                                    return RLI_ERRORE;
725
 
                                }
726
 
                            }
727
 
                        }
728
 
 
729
 
                        idCorr = mask_patch_sup[i];
730
 
                        mask_patch_corr[i] = idCorr;
731
 
                    }
732
 
                }
733
 
                else {          /*current cell and previous cell are equals */
734
 
 
735
 
 
736
 
                    if ((corrCell == supCell) &&
737
 
                        (mask_patch_sup[i] != mask_patch_corr[i - 1])) {
738
 
                        int l;
739
 
 
740
 
                        mask_patch_corr[i] = mask_patch_sup[i];
741
 
                        l = i - 1;
742
 
                        while (l >= 0) {
743
 
                            if (mask_patch_corr[l] == idCorr) {
744
 
                                mask_patch_corr[l] = mask_patch_sup[i];
745
 
                                l--;
746
 
                            }
747
 
                            else {
748
 
                                l = (-1);
749
 
                            }
750
 
                        }
751
 
                        lastId--;
752
 
                        idCorr = mask_patch_sup[i];
753
 
                    }
754
 
                    else {
755
 
                        mask_patch_corr[i] = idCorr;
756
 
                    }
757
 
 
758
 
                }
759
 
            }
760
 
            else {              /*cell null or not to consider */
761
 
 
762
 
                mask_patch_corr[i] = 0;
763
 
            }
764
 
        }
765
 
        mask_patch_sup = mask_patch_corr;
766
 
    }
767
 
 
768
 
 
769
 
 
770
 
    if (area != 0) {
771
 
        if (albero == NULL) {
772
 
            albero = avlID_make(idCorr, uno);
773
 
            if (albero == NULL) {
774
 
                G_fatal_error("avlID_make error");
775
 
                return RLI_ERRORE;
776
 
            }
777
 
            npatch++;
778
 
        }
779
 
        else {
780
 
            ris = avlID_add(&albero, idCorr, uno);
781
 
            switch (ris) {
782
 
            case AVL_ERR:
783
 
                {
784
 
                    G_fatal_error("avlID_add error");
785
 
                    return RLI_ERRORE;
786
 
                }
787
 
            case AVL_ADD:
788
 
                {
789
 
                    npatch++;
790
 
                    break;
791
 
                }
792
 
            case AVL_PRES:
793
 
                {
794
 
                    break;
795
 
                }
796
 
            default:
797
 
                {
798
 
                    G_fatal_error("avlID_add unknown error");
799
 
                    return RLI_ERRORE;
800
 
                }
801
 
            }
802
 
        }
803
 
 
804
 
 
805
 
        array = G_malloc(npatch * sizeof(avlID_tableRow));
806
 
        if (array == NULL) {
807
 
            G_fatal_error("malloc array failed");
808
 
            return RLI_ERRORE;
809
 
        }
810
 
        tot = avlID_to_array(albero, zero, array);
811
 
 
812
 
        if (tot != npatch) {
813
 
            G_warning
814
 
                ("avlID_to_array unaspected value. the result could be wrong");
815
 
            return RLI_ERRORE;
816
 
        }
817
 
 
818
 
        for (i = 0; i < npatch; i++) {
819
 
            if (array[i]->tot == 0)
820
 
                doppi++;
821
 
 
822
 
        }
823
 
        npatch = npatch - doppi;
824
 
 
825
 
        /*calculate distance */
826
 
        G_begin_distance_calculations();
827
 
        /* EW Dist at North edge */
828
 
        EW_DIST1 = G_distance(hd.east, hd.north, hd.west, hd.north);
829
 
        /* EW Dist at South Edge */
830
 
        EW_DIST2 = G_distance(hd.east, hd.south, hd.west, hd.south);
831
 
        /* NS Dist at East edge */
832
 
        NS_DIST1 = G_distance(hd.east, hd.north, hd.east, hd.south);
833
 
        /* NS Dist at West edge */
834
 
        NS_DIST2 = G_distance(hd.west, hd.north, hd.west, hd.south);
835
 
 
836
 
        areaCorrect = (((EW_DIST1 + EW_DIST2) / 2) / hd.cols) *
837
 
            (((NS_DIST1 + NS_DIST2) / 2) / hd.rows) * (area);
838
 
        indice = areaCorrect / npatch;
839
 
        G_free(array);
840
 
    }
841
 
    else
842
 
        indice = (double)(0);
843
 
 
844
 
 
845
 
    *result = indice;
846
 
 
847
 
 
848
 
    if (masked)
849
 
        G_free(mask_buf);
850
 
 
851
 
    G_free(mask_patch_corr);
852
 
 
853
 
    return RLI_OK;
854
 
}
855
 
 
856
 
 
857
 
int calculateF(int fd, area_des ad, struct Cell_head hd, double *result)
858
 
{
859
 
    FCELL *buf;
860
 
    FCELL *buf_sup;
861
 
 
862
 
    FCELL corrCell;
863
 
    FCELL precCell;
864
 
    FCELL supCell;
865
 
 
866
 
    int i, j;
867
 
    int mask_fd = -1, *mask_buf;
868
 
    int ris = 0;
869
 
    int masked = FALSE;
870
 
 
871
 
 
872
 
    long npatch = 0;
873
 
    long tot = 0;
874
 
    long zero = 0;
875
 
    long uno = 1;
876
 
    long idCorr = 0;
877
 
    long lastId = 0;
878
 
    long doppi = 0;
879
 
    long *mask_patch_sup;
880
 
    long *mask_patch_corr;
881
 
 
882
 
    double indice = 0;
883
 
    double area = 0;            /*if all cells are null area=0 */
884
 
    double areaCorrect = 0;
885
 
    double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
886
 
 
887
 
    avlID_tree albero = NULL;
888
 
 
889
 
    avlID_table *array;
890
 
 
891
 
 
892
 
 
893
 
    /* open mask if needed */
894
 
    if (ad->mask == 1) {
895
 
        if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
896
 
            return RLI_ERRORE;
897
 
        mask_buf = G_malloc(ad->cl * sizeof(int));
898
 
        if (mask_buf == NULL) {
899
 
            G_fatal_error("malloc mask_buf failed");
900
 
            return RLI_ERRORE;
901
 
        }
902
 
        masked = TRUE;
903
 
    }
904
 
 
905
 
    mask_patch_sup = G_malloc(ad->cl * sizeof(long));
906
 
    if (mask_patch_sup == NULL) {
907
 
        G_fatal_error("malloc mask_patch_sup failed");
908
 
        return RLI_ERRORE;
909
 
    }
910
 
 
911
 
    mask_patch_corr = G_malloc(ad->cl * sizeof(long));
912
 
    if (mask_patch_corr == NULL) {
913
 
        G_fatal_error("malloc mask_patch_corr failed");
914
 
        return RLI_ERRORE;
915
 
    }
916
 
 
917
 
    buf_sup = G_allocate_f_raster_buf();
918
 
    if (buf_sup == NULL) {
919
 
        G_fatal_error("malloc buf_sup failed");
920
 
        return RLI_ERRORE;
921
 
    }
922
 
 
923
 
 
924
 
    buf = G_allocate_f_raster_buf();
925
 
    if (buf == NULL) {
926
 
        G_fatal_error("malloc buf failed");
927
 
        return RLI_ERRORE;
928
 
    }
929
 
 
930
 
    G_set_f_null_value(buf_sup + ad->x, ad->cl);        /*the first time buf_sup is all null */
931
 
 
932
 
    for (i = 0; i < ad->cl; i++) {
933
 
        mask_patch_sup[i] = 0;
934
 
        mask_patch_corr[i] = 0;
935
 
    }
936
 
    /*for each raster row */
937
 
 
938
 
    for (j = 0; j < ad->rl; j++) {
939
 
        if (j > 0) {
940
 
            buf_sup = RLI_get_fcell_raster_row(fd, j - 1 + ad->y, ad);
941
 
        }
942
 
        buf = RLI_get_fcell_raster_row(fd, j + ad->y, ad);
943
 
        if (masked) {
944
 
            if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0) {
945
 
                G_fatal_error("mask read failed");
946
 
                return RLI_ERRORE;
947
 
            }
948
 
        }
949
 
        G_set_f_null_value(&precCell, 1);
950
 
        for (i = 0; i < ad->cl; i++) {  /*for each cell in the row */
951
 
 
952
 
 
953
 
            corrCell = buf[i + ad->x];
954
 
            if (((masked) && (mask_buf[i + ad->x] == 0))) {
955
 
                G_set_f_null_value(&corrCell, 1);
956
 
            }
957
 
            if (!(G_is_null_value(&corrCell, FCELL_TYPE))) {
958
 
                area++;
959
 
                if (i > 0)
960
 
                    precCell = buf[i - 1 + ad->x];
961
 
                if (j == 0)
962
 
                    G_set_f_null_value(&supCell, 1);
963
 
                else
964
 
                    supCell = buf_sup[i + ad->x];
965
 
 
966
 
 
967
 
                if (corrCell != precCell) {
968
 
                    if (corrCell != supCell) {
969
 
                        /*new patch */
970
 
                        if (idCorr == 0) {      /*first patch */
971
 
                            lastId = 1;
972
 
                            idCorr = 1;
973
 
                            mask_patch_corr[i] = idCorr;
974
 
                        }
975
 
                        else {  /*not first patch */
976
 
                            /* put in the tree previous value */
977
 
                            if (albero == NULL) {
978
 
                                albero = avlID_make(idCorr, uno);
979
 
                                if (albero == NULL) {
980
 
                                    G_fatal_error("avlID_make error");
981
 
                                    return RLI_ERRORE;
982
 
                                }
983
 
                                npatch++;
984
 
 
985
 
                            }
986
 
                            else {      /*tree not empty */
987
 
 
988
 
                                ris = avlID_add(&albero, idCorr, uno);
989
 
                                switch (ris) {
990
 
                                case AVL_ERR:
991
 
                                    {
992
 
                                        G_fatal_error("avlID_add error");
993
 
                                        return RLI_ERRORE;
994
 
                                    }
995
 
                                case AVL_ADD:
996
 
                                    {
997
 
                                        npatch++;
998
 
                                        break;
999
 
                                    }
1000
 
                                case AVL_PRES:
1001
 
                                    {
1002
 
                                        break;
1003
 
                                    }
1004
 
                                default:
1005
 
                                    {
1006
 
                                        G_fatal_error
1007
 
                                            ("avlID_add unknown error");
1008
 
                                        return RLI_ERRORE;
1009
 
                                    }
1010
 
                                }
1011
 
                            }
1012
 
 
1013
 
                            lastId++;
1014
 
                            idCorr = lastId;
1015
 
                            mask_patch_corr[i] = idCorr;
1016
 
                        }
1017
 
                    }
1018
 
                    else {      /*current cell and upper cell are equal */
1019
 
 
1020
 
                        if ((corrCell == precCell) &&
1021
 
                            (mask_patch_sup[i] != mask_patch_corr[i - 1])) {
1022
 
                            long r = 0;
1023
 
                            long del = mask_patch_sup[i];
1024
 
 
1025
 
                            r = avlID_sub(&albero, del);
1026
 
                            if (r == 0) {
1027
 
                                G_fatal_error("avlID_sub error");
1028
 
                                return RLI_ERRORE;
1029
 
                            }
1030
 
                            /*Remove one patch because it makes part of a patch already found */
1031
 
                            ris = avlID_add(&albero, idCorr, uno);
1032
 
                            switch (ris) {
1033
 
                            case AVL_ERR:
1034
 
                                {
1035
 
                                    G_fatal_error("avlID_add error");
1036
 
                                    return RLI_ERRORE;
1037
 
                                }
1038
 
                            case AVL_ADD:
1039
 
                                {
1040
 
                                    npatch++;
1041
 
                                    break;
1042
 
                                }
1043
 
                            case AVL_PRES:
1044
 
                                {
1045
 
                                    break;
1046
 
                                }
1047
 
                            default:
1048
 
                                {
1049
 
                                    G_fatal_error("avlID_add unknown error");
1050
 
                                    return RLI_ERRORE;
1051
 
                                }
1052
 
                            }
1053
 
                            r = i;
1054
 
                            while (i < ad->cl) {
1055
 
                                if (mask_patch_sup[r] == del) {
1056
 
                                    mask_patch_sup[r] = idCorr;
1057
 
                                }
1058
 
                                else {
1059
 
                                    r = ad->cl + 1;
1060
 
                                }
1061
 
                            }
1062
 
                        }
1063
 
 
1064
 
                        if (albero == NULL) {
1065
 
                            albero = avlID_make(idCorr, uno);
1066
 
                            if (albero == NULL) {
1067
 
                                G_fatal_error("avlID_make error");
1068
 
                                return RLI_ERRORE;
1069
 
                            }
1070
 
                            npatch++;
1071
 
                        }
1072
 
                        else {  /*the tree (albero) isn't null */
1073
 
 
1074
 
                            ris = avlID_add(&albero, idCorr, uno);
1075
 
                            switch (ris) {
1076
 
                            case AVL_ERR:
1077
 
                                {
1078
 
                                    G_fatal_error("avlID_add error");
1079
 
                                    return RLI_ERRORE;
1080
 
                                }
1081
 
                            case AVL_ADD:
1082
 
                                {
1083
 
                                    npatch++;
1084
 
                                    break;
1085
 
                                }
1086
 
                            case AVL_PRES:
1087
 
                                {
1088
 
                                    break;
1089
 
                                }
1090
 
                            default:
1091
 
                                {
1092
 
                                    G_fatal_error("avlID_add unknown error");
1093
 
                                    return RLI_ERRORE;
1094
 
                                }
1095
 
                            }
1096
 
                        }
1097
 
 
1098
 
                        idCorr = mask_patch_sup[i];
1099
 
                        mask_patch_corr[i] = idCorr;
1100
 
                    }
1101
 
                }
1102
 
                else {          /*current cell and previous cell are equal */
1103
 
 
1104
 
                    if ((corrCell == supCell) &&
1105
 
                        (mask_patch_sup[i] != mask_patch_corr[i - 1])) {
1106
 
                        int l;
1107
 
 
1108
 
                        mask_patch_corr[i] = mask_patch_sup[i];
1109
 
                        l = i - 1;
1110
 
                        while (l >= 0) {
1111
 
                            if (mask_patch_corr[l] == idCorr) {
1112
 
                                mask_patch_corr[l] = mask_patch_sup[i];
1113
 
                                l--;
1114
 
                            }
1115
 
                            else {
1116
 
                                l = (-1);
1117
 
                            }
1118
 
                        }
1119
 
                        lastId--;
1120
 
                        idCorr = mask_patch_sup[i];
1121
 
                    }
1122
 
                    else {
1123
 
                        mask_patch_corr[i] = idCorr;
1124
 
                    }
1125
 
 
1126
 
                }
1127
 
            }
1128
 
            else {              /*null cell or cell not to consider */
1129
 
 
1130
 
                mask_patch_corr[i] = 0;
1131
 
            }
1132
 
        }
1133
 
        mask_patch_sup = mask_patch_corr;
1134
 
    }
1135
 
 
1136
 
 
1137
 
 
1138
 
    if (area != 0) {
1139
 
        if (albero == NULL) {
1140
 
            albero = avlID_make(idCorr, uno);
1141
 
            if (albero == NULL) {
1142
 
                G_fatal_error("avlID_make error");
1143
 
                return RLI_ERRORE;
1144
 
            }
1145
 
            npatch++;
1146
 
        }
1147
 
        else {
1148
 
            ris = avlID_add(&albero, idCorr, uno);
1149
 
            switch (ris) {
1150
 
            case AVL_ERR:
1151
 
                {
1152
 
                    G_fatal_error("avlID_add error");
1153
 
                    return RLI_ERRORE;
1154
 
                }
1155
 
            case AVL_ADD:
1156
 
                {
1157
 
                    npatch++;
1158
 
                    break;
1159
 
                }
1160
 
            case AVL_PRES:
1161
 
                {
1162
 
                    break;
1163
 
                }
1164
 
            default:
1165
 
                {
1166
 
                    G_fatal_error("avlID_add unknown error");
1167
 
                    return RLI_ERRORE;
1168
 
                }
1169
 
            }
1170
 
        }
1171
 
 
1172
 
 
1173
 
        array = G_malloc(npatch * sizeof(avlID_tableRow));
1174
 
        if (array == NULL) {
1175
 
            G_fatal_error("malloc array failed");
1176
 
            return RLI_ERRORE;
1177
 
        }
1178
 
        tot = avlID_to_array(albero, zero, array);
1179
 
 
1180
 
        if (tot != npatch) {
1181
 
            G_warning
1182
 
                ("avlID_to_array unaspected value. the result could be wrong");
1183
 
            return RLI_ERRORE;
1184
 
        }
1185
 
 
1186
 
        for (i = 0; i < npatch; i++) {
1187
 
            if (array[i]->tot == 0)
1188
 
                doppi++;
1189
 
 
1190
 
        }
1191
 
        npatch = npatch - doppi;
1192
 
 
1193
 
        /*calculate distance */
1194
 
        G_begin_distance_calculations();
1195
 
        /* EW Dist at North edge */
1196
 
        EW_DIST1 = G_distance(hd.east, hd.north, hd.west, hd.north);
1197
 
        /* EW Dist at South Edge */
1198
 
        EW_DIST2 = G_distance(hd.east, hd.south, hd.west, hd.south);
1199
 
        /* NS Dist at East edge */
1200
 
        NS_DIST1 = G_distance(hd.east, hd.north, hd.east, hd.south);
1201
 
        /* NS Dist at West edge */
1202
 
        NS_DIST2 = G_distance(hd.west, hd.north, hd.west, hd.south);
1203
 
 
1204
 
        areaCorrect = (((EW_DIST1 + EW_DIST2) / 2) / hd.cols) *
1205
 
            (((NS_DIST1 + NS_DIST2) / 2) / hd.rows) * (area);
1206
 
        indice = areaCorrect / npatch;
1207
 
        G_free(array);
1208
 
    }
1209
 
    else
1210
 
        indice = (double)(0);
1211
 
 
1212
 
 
1213
 
    *result = indice;
1214
 
 
1215
 
 
1216
 
    if (masked)
1217
 
        G_free(mask_buf);
1218
 
 
1219
 
    G_free(mask_patch_corr);
 
108
 
 
109
    return RLI_OK;
 
110
}
 
111
 
 
112
 
 
113
int calculate(int fd, struct area_entry *ad, double *result)
 
114
{
 
115
    CELL *buf, *buf_sup, *buf_null;
 
116
    CELL corrCell, precCell, supCell;
 
117
    long npatch, area; 
 
118
    long pid, old_pid, new_pid, *pid_corr, *pid_sup, *ltmp;
 
119
    struct pst *pst;
 
120
    long nalloc, incr;
 
121
    int i, j, k;
 
122
    int connected;
 
123
    int mask_fd, *mask_buf, *mask_sup, *mask_tmp, masked;
 
124
    struct Cell_head hd;
 
125
 
 
126
    Rast_get_window(&hd);
 
127
 
 
128
    buf_null = Rast_allocate_c_buf();
 
129
    Rast_set_c_null_value(buf_null, Rast_window_cols());
 
130
    buf_sup = buf_null;
 
131
 
 
132
    /* initialize patch ids */
 
133
    pid_corr = G_malloc(ad->cl * sizeof(long));
 
134
    pid_sup = G_malloc(ad->cl * sizeof(long));
 
135
 
 
136
    for (j = 0; j < ad->cl; j++) {
 
137
        pid_corr[j] = 0;
 
138
        pid_sup[j] = 0;
 
139
    }
 
140
 
 
141
    /* open mask if needed */
 
142
    mask_fd = -1;
 
143
    mask_buf = mask_sup = NULL;
 
144
    masked = FALSE;
 
145
    if (ad->mask == 1) {
 
146
        if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
 
147
            return RLI_ERRORE;
 
148
        mask_buf = G_malloc(ad->cl * sizeof(int));
 
149
        if (mask_buf == NULL) {
 
150
            G_fatal_error("malloc mask_buf failed");
 
151
            return RLI_ERRORE;
 
152
        }
 
153
        mask_sup = G_malloc(ad->cl * sizeof(int));
 
154
        if (mask_sup == NULL) {
 
155
            G_fatal_error("malloc mask_buf failed");
 
156
            return RLI_ERRORE;
 
157
        }
 
158
        for (j = 0; j < ad->cl; j++)
 
159
            mask_buf[j] = 0;
 
160
 
 
161
        masked = TRUE;
 
162
    }
 
163
 
 
164
    /* calculate number of patches */
 
165
    npatch = 0;
 
166
    area = 0;
 
167
    pid = 0;
 
168
 
 
169
    /* patch size and type */
 
170
    incr = 1024;
 
171
    if (incr > ad->rl)
 
172
        incr = ad->rl;
 
173
    if (incr > ad->cl)
 
174
        incr = ad->cl;
 
175
    if (incr < 2)
 
176
        incr = 2;
 
177
    nalloc = incr;
 
178
    pst = G_malloc(nalloc * sizeof(struct pst));
 
179
    for (k = 0; k < nalloc; k++) {
 
180
        pst[k].count = 0;
 
181
    }
 
182
 
 
183
    for (i = 0; i < ad->rl; i++) {
 
184
        buf = RLI_get_cell_raster_row(fd, i + ad->y, ad);
 
185
        if (i > 0) {
 
186
            buf_sup = RLI_get_cell_raster_row(fd, i - 1 + ad->y, ad);
 
187
        }
 
188
 
 
189
        if (masked) {
 
190
            mask_tmp = mask_sup;
 
191
            mask_sup = mask_buf;
 
192
            mask_buf = mask_tmp;
 
193
            if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0)
 
194
                return 0;
 
195
        }
 
196
 
 
197
        ltmp = pid_sup;
 
198
        pid_sup = pid_corr;
 
199
        pid_corr = ltmp;
 
200
 
 
201
        Rast_set_c_null_value(&precCell, 1);
 
202
 
 
203
        connected = 0;
 
204
        for (j = 0; j < ad->cl; j++) {
 
205
            pid_corr[j] = 0;
 
206
            
 
207
            corrCell = buf[j + ad->x];
 
208
            if (masked && (mask_buf[j] == 0)) {
 
209
                Rast_set_c_null_value(&corrCell, 1);
 
210
            }
 
211
 
 
212
            if (Rast_is_c_null_value(&corrCell)) {
 
213
                connected = 0;
 
214
                precCell = corrCell;
 
215
                continue;
 
216
            }
 
217
            
 
218
            area++;
 
219
            
 
220
            supCell = buf_sup[j + ad->x];
 
221
            if (masked && (mask_sup[j] == 0)) {
 
222
                Rast_set_c_null_value(&supCell, 1);
 
223
            }
 
224
 
 
225
            if (!Rast_is_c_null_value(&precCell) && corrCell == precCell) {
 
226
                pid_corr[j] = pid_corr[j - 1];
 
227
                connected = 1;
 
228
                pst[pid_corr[j]].count++;
 
229
            }
 
230
            else {
 
231
                connected = 0;
 
232
            }
 
233
 
 
234
            if (!Rast_is_c_null_value(&supCell) && corrCell == supCell) {
 
235
 
 
236
                if (pid_corr[j] != pid_sup[j]) {
 
237
                    /* connect or merge */
 
238
                    /* after r.clump */
 
239
                    if (connected) {
 
240
                        npatch--;
 
241
 
 
242
                        if (npatch == 0) {
 
243
                            G_fatal_error("npatch == 0 at row %d, col %d", i, j);
 
244
                        }
 
245
                    }
 
246
 
 
247
                    old_pid = pid_corr[j];
 
248
                    new_pid = pid_sup[j];
 
249
                    pid_corr[j] = new_pid;
 
250
                    if (old_pid > 0) {
 
251
                        /* merge */
 
252
                        /* update left side of the current row */
 
253
                        for (k = 0; k < j; k++) {
 
254
                            if (pid_corr[k] == old_pid)
 
255
                                pid_corr[k] = new_pid;
 
256
                        }
 
257
                        /* update right side of the previous row */
 
258
                        for (k = j + 1; k < ad->cl; k++) {
 
259
                            if (pid_sup[k] == old_pid)
 
260
                                pid_sup[k] = new_pid;
 
261
                        }
 
262
                        pst[new_pid].count += pst[old_pid].count;
 
263
                        pst[old_pid].count = 0;
 
264
                        
 
265
                        if (old_pid == pid)
 
266
                            pid--;
 
267
                    }
 
268
                    else {
 
269
                        pst[new_pid].count++;
 
270
                    }
 
271
                }
 
272
                connected = 1;
 
273
            }
 
274
 
 
275
            if (!connected) {
 
276
                /* start new patch */
 
277
                npatch++;
 
278
                pid++;
 
279
                pid_corr[j] = pid;
 
280
 
 
281
                if (pid >= nalloc) {
 
282
                    pst = (struct pst *)G_realloc(pst, (pid + incr) * sizeof(struct pst));
 
283
 
 
284
                    for (k = nalloc; k < pid + incr; k++)
 
285
                        pst[k].count = 0;
 
286
                        
 
287
                    nalloc = pid + incr;
 
288
                }
 
289
 
 
290
                pst[pid].count = 1;
 
291
                pst[pid].type.t = CELL_TYPE;
 
292
                pst[pid].type.val.c = corrCell;
 
293
            }
 
294
            precCell = corrCell;
 
295
        }
 
296
    }
 
297
 
 
298
    if (npatch > 0) {
 
299
        double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
 
300
        double area_units;
 
301
 
 
302
        /* calculate distance */
 
303
        G_begin_distance_calculations();
 
304
        /* EW Dist at North edge */
 
305
        EW_DIST1 = G_distance(hd.east, hd.north, hd.west, hd.north);
 
306
        /* EW Dist at South Edge */
 
307
        EW_DIST2 = G_distance(hd.east, hd.south, hd.west, hd.south);
 
308
        /* NS Dist at East edge */
 
309
        NS_DIST1 = G_distance(hd.east, hd.north, hd.east, hd.south);
 
310
        /* NS Dist at West edge */
 
311
        NS_DIST2 = G_distance(hd.west, hd.north, hd.west, hd.south);
 
312
 
 
313
        area_units = (((EW_DIST1 + EW_DIST2) / 2) / hd.cols) *
 
314
            (((NS_DIST1 + NS_DIST2) / 2) / hd.rows) * area;
 
315
 
 
316
        *result = area_units / (npatch * 10000);
 
317
    }
 
318
    else {
 
319
        *result = 0;
 
320
    }
 
321
 
 
322
    if (masked) {
 
323
        close(mask_fd);
 
324
        G_free(mask_buf);
 
325
        G_free(mask_sup);
 
326
    }
 
327
    G_free(buf_null);
 
328
    G_free(pid_corr);
 
329
    G_free(pid_sup);
 
330
    G_free(pst);
 
331
 
 
332
    return RLI_OK;
 
333
}
 
334
 
 
335
 
 
336
int calculateD(int fd, struct area_entry *ad, double *result)
 
337
{
 
338
    DCELL *buf, *buf_sup, *buf_null;
 
339
    DCELL corrCell, precCell, supCell;
 
340
    long npatch, area; 
 
341
    long pid, old_pid, new_pid, *pid_corr, *pid_sup, *ltmp;
 
342
    struct pst *pst;
 
343
    long nalloc, incr;
 
344
    int i, j, k;
 
345
    int connected;
 
346
    int mask_fd, *mask_buf, *mask_sup, *mask_tmp, masked;
 
347
    struct Cell_head hd;
 
348
 
 
349
    Rast_get_window(&hd);
 
350
 
 
351
    buf_null = Rast_allocate_d_buf();
 
352
    Rast_set_d_null_value(buf_null, Rast_window_cols());
 
353
    buf_sup = buf_null;
 
354
 
 
355
    /* initialize patch ids */
 
356
    pid_corr = G_malloc(ad->cl * sizeof(long));
 
357
    pid_sup = G_malloc(ad->cl * sizeof(long));
 
358
 
 
359
    for (j = 0; j < ad->cl; j++) {
 
360
        pid_corr[j] = 0;
 
361
        pid_sup[j] = 0;
 
362
    }
 
363
 
 
364
    /* open mask if needed */
 
365
    mask_fd = -1;
 
366
    mask_buf = mask_sup = NULL;
 
367
    masked = FALSE;
 
368
    if (ad->mask == 1) {
 
369
        if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
 
370
            return RLI_ERRORE;
 
371
        mask_buf = G_malloc(ad->cl * sizeof(int));
 
372
        if (mask_buf == NULL) {
 
373
            G_fatal_error("malloc mask_buf failed");
 
374
            return RLI_ERRORE;
 
375
        }
 
376
        mask_sup = G_malloc(ad->cl * sizeof(int));
 
377
        if (mask_sup == NULL) {
 
378
            G_fatal_error("malloc mask_buf failed");
 
379
            return RLI_ERRORE;
 
380
        }
 
381
        for (j = 0; j < ad->cl; j++)
 
382
            mask_buf[j] = 0;
 
383
 
 
384
        masked = TRUE;
 
385
    }
 
386
 
 
387
    /* calculate number of patches */
 
388
    npatch = 0;
 
389
    area = 0;
 
390
    pid = 0;
 
391
 
 
392
    /* patch size and type */
 
393
    incr = 1024;
 
394
    if (incr > ad->rl)
 
395
        incr = ad->rl;
 
396
    if (incr > ad->cl)
 
397
        incr = ad->cl;
 
398
    if (incr < 2)
 
399
        incr = 2;
 
400
    nalloc = incr;
 
401
    pst = G_malloc(nalloc * sizeof(struct pst));
 
402
    for (k = 0; k < nalloc; k++) {
 
403
        pst[k].count = 0;
 
404
    }
 
405
 
 
406
    for (i = 0; i < ad->rl; i++) {
 
407
        buf = RLI_get_dcell_raster_row(fd, i + ad->y, ad);
 
408
        if (i > 0) {
 
409
            buf_sup = RLI_get_dcell_raster_row(fd, i - 1 + ad->y, ad);
 
410
        }
 
411
 
 
412
        if (masked) {
 
413
            mask_tmp = mask_sup;
 
414
            mask_sup = mask_buf;
 
415
            mask_buf = mask_tmp;
 
416
            if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0)
 
417
                return 0;
 
418
        }
 
419
        
 
420
        ltmp = pid_sup;
 
421
        pid_sup = pid_corr;
 
422
        pid_corr = ltmp;
 
423
 
 
424
        Rast_set_d_null_value(&precCell, 1);
 
425
 
 
426
        connected = 0;
 
427
        for (j = 0; j < ad->cl; j++) {
 
428
            pid_corr[j] = 0;
 
429
            
 
430
            corrCell = buf[j + ad->x];
 
431
            if (masked && (mask_buf[j] == 0)) {
 
432
                Rast_set_d_null_value(&corrCell, 1);
 
433
            }
 
434
 
 
435
            if (Rast_is_d_null_value(&corrCell)) {
 
436
                connected = 0;
 
437
                precCell = corrCell;
 
438
                continue;
 
439
            }
 
440
            
 
441
            area++;
 
442
            
 
443
            supCell = buf_sup[j + ad->x];
 
444
            if (masked && (mask_sup[j] == 0)) {
 
445
                Rast_set_d_null_value(&supCell, 1);
 
446
            }
 
447
 
 
448
            if (!Rast_is_d_null_value(&precCell) && corrCell == precCell) {
 
449
                pid_corr[j] = pid_corr[j - 1];
 
450
                connected = 1;
 
451
                pst[pid_corr[j]].count++;
 
452
            }
 
453
            else {
 
454
                connected = 0;
 
455
            }
 
456
 
 
457
            if (!Rast_is_d_null_value(&supCell) && corrCell == supCell) {
 
458
 
 
459
                if (pid_corr[j] != pid_sup[j]) {
 
460
                    /* connect or merge */
 
461
                    /* after r.clump */
 
462
                    if (connected) {
 
463
                        npatch--;
 
464
 
 
465
                        if (npatch == 0) {
 
466
                            G_fatal_error("npatch == 0 at row %d, col %d", i, j);
 
467
                        }
 
468
                    }
 
469
 
 
470
                    old_pid = pid_corr[j];
 
471
                    new_pid = pid_sup[j];
 
472
                    pid_corr[j] = new_pid;
 
473
                    if (old_pid > 0) {
 
474
                        /* merge */
 
475
                        /* update left side of the current row */
 
476
                        for (k = 0; k < j; k++) {
 
477
                            if (pid_corr[k] == old_pid)
 
478
                                pid_corr[k] = new_pid;
 
479
                        }
 
480
                        /* update right side of the previous row */
 
481
                        for (k = j + 1; k < ad->cl; k++) {
 
482
                            if (pid_sup[k] == old_pid)
 
483
                                pid_sup[k] = new_pid;
 
484
                        }
 
485
                        pst[new_pid].count += pst[old_pid].count;
 
486
                        pst[old_pid].count = 0;
 
487
                        
 
488
                        if (old_pid == pid)
 
489
                            pid--;
 
490
                    }
 
491
                    else {
 
492
                        pst[new_pid].count++;
 
493
                    }
 
494
                }
 
495
                connected = 1;
 
496
            }
 
497
 
 
498
            if (!connected) {
 
499
                /* start new patch */
 
500
                npatch++;
 
501
                pid++;
 
502
                pid_corr[j] = pid;
 
503
 
 
504
                if (pid >= nalloc) {
 
505
                    pst = (struct pst *)G_realloc(pst, (pid + incr) * sizeof(struct pst));
 
506
 
 
507
                    for (k = nalloc; k < pid + incr; k++)
 
508
                        pst[k].count = 0;
 
509
                        
 
510
                    nalloc = pid + incr;
 
511
                }
 
512
 
 
513
                pst[pid].count = 1;
 
514
                pst[pid].type.t = CELL_TYPE;
 
515
                pst[pid].type.val.c = corrCell;
 
516
            }
 
517
            precCell = corrCell;
 
518
        }
 
519
    }
 
520
 
 
521
    if (npatch > 0) {
 
522
        double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
 
523
        double area_units;
 
524
 
 
525
        /* calculate distance */
 
526
        G_begin_distance_calculations();
 
527
        /* EW Dist at North edge */
 
528
        EW_DIST1 = G_distance(hd.east, hd.north, hd.west, hd.north);
 
529
        /* EW Dist at South Edge */
 
530
        EW_DIST2 = G_distance(hd.east, hd.south, hd.west, hd.south);
 
531
        /* NS Dist at East edge */
 
532
        NS_DIST1 = G_distance(hd.east, hd.north, hd.east, hd.south);
 
533
        /* NS Dist at West edge */
 
534
        NS_DIST2 = G_distance(hd.west, hd.north, hd.west, hd.south);
 
535
 
 
536
        area_units = (((EW_DIST1 + EW_DIST2) / 2) / hd.cols) *
 
537
            (((NS_DIST1 + NS_DIST2) / 2) / hd.rows) * area;
 
538
 
 
539
        *result = area_units / (npatch * 10000);
 
540
    }
 
541
    else {
 
542
        *result = 0;
 
543
    }
 
544
 
 
545
    if (masked) {
 
546
        close(mask_fd);
 
547
        G_free(mask_buf);
 
548
        G_free(mask_sup);
 
549
    }
 
550
    G_free(buf_null);
 
551
    G_free(pid_corr);
 
552
    G_free(pid_sup);
 
553
    G_free(pst);
 
554
 
 
555
    return RLI_OK;
 
556
}
 
557
 
 
558
 
 
559
int calculateF(int fd, struct area_entry *ad, double *result)
 
560
{
 
561
    FCELL *buf, *buf_sup, *buf_null;
 
562
    FCELL corrCell, precCell, supCell;
 
563
    long npatch, area; 
 
564
    long pid, old_pid, new_pid, *pid_corr, *pid_sup, *ltmp;
 
565
    struct pst *pst;
 
566
    long nalloc, incr;
 
567
    int i, j, k;
 
568
    int connected;
 
569
    int mask_fd, *mask_buf, *mask_sup, *mask_tmp, masked;
 
570
    struct Cell_head hd;
 
571
 
 
572
    Rast_get_window(&hd);
 
573
 
 
574
    buf_null = Rast_allocate_f_buf();
 
575
    Rast_set_f_null_value(buf_null, Rast_window_cols());
 
576
    buf_sup = buf_null;
 
577
 
 
578
    /* initialize patch ids */
 
579
    pid_corr = G_malloc(ad->cl * sizeof(long));
 
580
    pid_sup = G_malloc(ad->cl * sizeof(long));
 
581
 
 
582
    for (j = 0; j < ad->cl; j++) {
 
583
        pid_corr[j] = 0;
 
584
        pid_sup[j] = 0;
 
585
    }
 
586
 
 
587
    /* open mask if needed */
 
588
    mask_fd = -1;
 
589
    mask_buf = mask_sup = NULL;
 
590
    masked = FALSE;
 
591
    if (ad->mask == 1) {
 
592
        if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
 
593
            return RLI_ERRORE;
 
594
        mask_buf = G_malloc(ad->cl * sizeof(int));
 
595
        if (mask_buf == NULL) {
 
596
            G_fatal_error("malloc mask_buf failed");
 
597
            return RLI_ERRORE;
 
598
        }
 
599
        mask_sup = G_malloc(ad->cl * sizeof(int));
 
600
        if (mask_sup == NULL) {
 
601
            G_fatal_error("malloc mask_buf failed");
 
602
            return RLI_ERRORE;
 
603
        }
 
604
        for (j = 0; j < ad->cl; j++)
 
605
            mask_buf[j] = 0;
 
606
 
 
607
        masked = TRUE;
 
608
    }
 
609
 
 
610
    /* calculate number of patches */
 
611
    npatch = 0;
 
612
    area = 0;
 
613
    pid = 0;
 
614
 
 
615
    /* patch size and type */
 
616
    incr = 1024;
 
617
    if (incr > ad->rl)
 
618
        incr = ad->rl;
 
619
    if (incr > ad->cl)
 
620
        incr = ad->cl;
 
621
    if (incr < 2)
 
622
        incr = 2;
 
623
    nalloc = incr;
 
624
    pst = G_malloc(nalloc * sizeof(struct pst));
 
625
    for (k = 0; k < nalloc; k++) {
 
626
        pst[k].count = 0;
 
627
    }
 
628
 
 
629
    for (i = 0; i < ad->rl; i++) {
 
630
        buf = RLI_get_fcell_raster_row(fd, i + ad->y, ad);
 
631
        if (i > 0) {
 
632
            buf_sup = RLI_get_fcell_raster_row(fd, i - 1 + ad->y, ad);
 
633
        }
 
634
 
 
635
        if (masked) {
 
636
            mask_tmp = mask_sup;
 
637
            mask_sup = mask_buf;
 
638
            mask_buf = mask_tmp;
 
639
            if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0)
 
640
                return 0;
 
641
        }
 
642
        
 
643
        ltmp = pid_sup;
 
644
        pid_sup = pid_corr;
 
645
        pid_corr = ltmp;
 
646
 
 
647
        Rast_set_f_null_value(&precCell, 1);
 
648
 
 
649
        connected = 0;
 
650
        for (j = 0; j < ad->cl; j++) {
 
651
            pid_corr[j] = 0;
 
652
            
 
653
            corrCell = buf[j + ad->x];
 
654
            if (masked && (mask_buf[j] == 0)) {
 
655
                Rast_set_f_null_value(&corrCell, 1);
 
656
            }
 
657
 
 
658
            if (Rast_is_f_null_value(&corrCell)) {
 
659
                connected = 0;
 
660
                precCell = corrCell;
 
661
                continue;
 
662
            }
 
663
            
 
664
            area++;
 
665
            
 
666
            supCell = buf_sup[j + ad->x];
 
667
            if (masked && (mask_sup[j] == 0)) {
 
668
                Rast_set_f_null_value(&supCell, 1);
 
669
            }
 
670
 
 
671
            if (!Rast_is_f_null_value(&precCell) && corrCell == precCell) {
 
672
                pid_corr[j] = pid_corr[j - 1];
 
673
                connected = 1;
 
674
                pst[pid_corr[j]].count++;
 
675
            }
 
676
            else {
 
677
                connected = 0;
 
678
            }
 
679
 
 
680
            if (!Rast_is_f_null_value(&supCell) && corrCell == supCell) {
 
681
 
 
682
                if (pid_corr[j] != pid_sup[j]) {
 
683
                    /* connect or merge */
 
684
                    /* after r.clump */
 
685
                    if (connected) {
 
686
                        npatch--;
 
687
 
 
688
                        if (npatch == 0) {
 
689
                            G_fatal_error("npatch == 0 at row %d, col %d", i, j);
 
690
                        }
 
691
                    }
 
692
 
 
693
                    old_pid = pid_corr[j];
 
694
                    new_pid = pid_sup[j];
 
695
                    pid_corr[j] = new_pid;
 
696
                    if (old_pid > 0) {
 
697
                        /* merge */
 
698
                        /* update left side of the current row */
 
699
                        for (k = 0; k < j; k++) {
 
700
                            if (pid_corr[k] == old_pid)
 
701
                                pid_corr[k] = new_pid;
 
702
                        }
 
703
                        /* update right side of the previous row */
 
704
                        for (k = j + 1; k < ad->cl; k++) {
 
705
                            if (pid_sup[k] == old_pid)
 
706
                                pid_sup[k] = new_pid;
 
707
                        }
 
708
                        pst[new_pid].count += pst[old_pid].count;
 
709
                        pst[old_pid].count = 0;
 
710
                        
 
711
                        if (old_pid == pid)
 
712
                            pid--;
 
713
                    }
 
714
                    else {
 
715
                        pst[new_pid].count++;
 
716
                    }
 
717
                }
 
718
                connected = 1;
 
719
            }
 
720
 
 
721
            if (!connected) {
 
722
                /* start new patch */
 
723
                npatch++;
 
724
                pid++;
 
725
                pid_corr[j] = pid;
 
726
 
 
727
                if (pid >= nalloc) {
 
728
                    pst = (struct pst *)G_realloc(pst, (pid + incr) * sizeof(struct pst));
 
729
 
 
730
                    for (k = nalloc; k < pid + incr; k++)
 
731
                        pst[k].count = 0;
 
732
                        
 
733
                    nalloc = pid + incr;
 
734
                }
 
735
 
 
736
                pst[pid].count = 1;
 
737
                pst[pid].type.t = CELL_TYPE;
 
738
                pst[pid].type.val.c = corrCell;
 
739
            }
 
740
            precCell = corrCell;
 
741
        }
 
742
    }
 
743
 
 
744
    if (npatch > 0) {
 
745
        double EW_DIST1, EW_DIST2, NS_DIST1, NS_DIST2;
 
746
        double area_units;
 
747
 
 
748
        /* calculate distance */
 
749
        G_begin_distance_calculations();
 
750
        /* EW Dist at North edge */
 
751
        EW_DIST1 = G_distance(hd.east, hd.north, hd.west, hd.north);
 
752
        /* EW Dist at South Edge */
 
753
        EW_DIST2 = G_distance(hd.east, hd.south, hd.west, hd.south);
 
754
        /* NS Dist at East edge */
 
755
        NS_DIST1 = G_distance(hd.east, hd.north, hd.east, hd.south);
 
756
        /* NS Dist at West edge */
 
757
        NS_DIST2 = G_distance(hd.west, hd.north, hd.west, hd.south);
 
758
 
 
759
        area_units = (((EW_DIST1 + EW_DIST2) / 2) / hd.cols) *
 
760
            (((NS_DIST1 + NS_DIST2) / 2) / hd.rows) * area;
 
761
 
 
762
        *result = area_units / (npatch * 10000);
 
763
    }
 
764
    else {
 
765
        *result = 0;
 
766
    }
 
767
 
 
768
    if (masked) {
 
769
        close(mask_fd);
 
770
        G_free(mask_buf);
 
771
        G_free(mask_sup);
 
772
    }
 
773
    G_free(buf_null);
 
774
    G_free(pid_corr);
 
775
    G_free(pid_sup);
 
776
    G_free(pst);
1220
777
 
1221
778
    return RLI_OK;
1222
779
}