padrange.c 24 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279
  1. /*
  2. * \brief calculates range of patch area size
  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. */
  11. #include <grass/gis.h>
  12. #include <grass/raster.h>
  13. #include <grass/glocale.h>
  14. #include <stdlib.h>
  15. #include <fcntl.h>
  16. #include <math.h>
  17. #include "../r.li.daemon/defs.h"
  18. #include "../r.li.daemon/avlDefs.h"
  19. #include "../r.li.daemon/avlID.h"
  20. #include "../r.li.daemon/GenericCell.h"
  21. #include "../r.li.daemon/daemon.h"
  22. int calculate(int fd, struct area_entry *ad, double *result);
  23. int calculateD(int fd, struct area_entry *ad, double *result);
  24. int calculateF(int fd, struct area_entry *ad, double *result);
  25. int main(int argc, char *argv[])
  26. {
  27. struct Option *raster, *conf, *output;
  28. struct GModule *module;
  29. G_gisinit(argv[0]);
  30. module = G_define_module();
  31. module->description =
  32. _("Calculates range of patch area size on a raster map");
  33. G_add_keyword(_("raster"));
  34. G_add_keyword(_("landscape structure analysis"));
  35. G_add_keyword(_("patch index"));
  36. /* define options */
  37. raster = G_define_standard_option(G_OPT_R_INPUT);
  38. conf = G_define_standard_option(G_OPT_F_INPUT);
  39. conf->key = "config";
  40. conf->description = _("Configuration file");
  41. conf->required = YES;
  42. output = G_define_standard_option(G_OPT_R_OUTPUT);
  43. if (G_parser(argc, argv))
  44. exit(EXIT_FAILURE);
  45. return calculateIndex(conf->answer, patchAreaDistributionRANGE, NULL,
  46. raster->answer, output->answer);
  47. }
  48. int patchAreaDistributionRANGE(int fd, char **par, struct area_entry *ad,
  49. double *result)
  50. {
  51. double indice = 0;
  52. struct Cell_head hd;
  53. int ris = RLI_OK;
  54. Rast_get_cellhd(ad->raster, "", &hd);
  55. switch (ad->data_type) {
  56. case CELL_TYPE:
  57. {
  58. ris = calculate(fd, ad, &indice);
  59. break;
  60. }
  61. case DCELL_TYPE:
  62. {
  63. ris = calculateD(fd, ad, &indice);
  64. break;
  65. }
  66. case FCELL_TYPE:
  67. {
  68. ris = calculateF(fd, ad, &indice);
  69. break;
  70. }
  71. default:
  72. {
  73. G_fatal_error("data type unknown");
  74. return RLI_ERRORE;
  75. }
  76. }
  77. if (ris != RLI_OK) {
  78. *result = -1;
  79. return RLI_ERRORE;
  80. }
  81. *result = indice;
  82. return RLI_OK;
  83. }
  84. int calculate(int fd, struct area_entry *ad, double *result)
  85. {
  86. CELL *buf;
  87. CELL *buf_sup;
  88. CELL corrCell;
  89. CELL precCell;
  90. CELL supCell;
  91. int i, j;
  92. int mask_fd = -1, *mask_buf;
  93. int ris = 0;
  94. int masked = FALSE;
  95. int areaPatch = 0; /*if all cells are null areaPatch=0 */
  96. int area_max = 0, area_min = 0;
  97. long npatch = 0;
  98. long tot = 0;
  99. long zero = 0;
  100. long totCorr = 0;
  101. long idCorr = 0;
  102. long lastId = 0;
  103. long *mask_patch_sup;
  104. long *mask_patch_corr;
  105. double indice = 0;
  106. double rk = 0;
  107. avlID_tree albero = NULL;
  108. avlID_table *array = NULL;
  109. generic_cell gc;
  110. gc.t = CELL_TYPE;
  111. /* open mask if needed */
  112. if (ad->mask == 1) {
  113. if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
  114. return RLI_ERRORE;
  115. mask_buf = G_malloc(ad->cl * sizeof(int));
  116. if (mask_buf == NULL) {
  117. G_fatal_error("malloc mask_buf failed");
  118. return RLI_ERRORE;
  119. }
  120. masked = TRUE;
  121. }
  122. mask_patch_sup = G_malloc(ad->cl * sizeof(long));
  123. if (mask_patch_sup == NULL) {
  124. G_fatal_error("malloc mask_patch_sup failed");
  125. return RLI_ERRORE;
  126. }
  127. mask_patch_corr = G_malloc(ad->cl * sizeof(long));
  128. if (mask_patch_corr == NULL) {
  129. G_fatal_error("malloc mask_patch_corr failed");
  130. return RLI_ERRORE;
  131. }
  132. buf_sup = Rast_allocate_c_buf();
  133. if (buf_sup == NULL) {
  134. G_fatal_error("malloc buf_sup failed");
  135. return RLI_ERRORE;
  136. }
  137. buf = Rast_allocate_c_buf();
  138. if (buf == NULL) {
  139. G_fatal_error("malloc buf failed");
  140. return RLI_ERRORE;
  141. }
  142. Rast_set_c_null_value(buf_sup + ad->x, ad->cl); /*the first time buf_sup is all null */
  143. for (i = 0; i < ad->cl; i++) {
  144. mask_patch_sup[i] = 0;
  145. mask_patch_corr[i] = 0;
  146. }
  147. for (j = 0; j < ad->rl; j++)
  148. /*for each raster row */
  149. {
  150. if (j > 0) {
  151. buf_sup = RLI_get_cell_raster_row(fd, j - 1 + ad->y, ad);
  152. }
  153. buf = RLI_get_cell_raster_row(fd, j + ad->y, ad);
  154. if (masked) {
  155. if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0) {
  156. G_fatal_error("mask read failed");
  157. return RLI_ERRORE;
  158. }
  159. }
  160. Rast_set_c_null_value(&precCell, 1);
  161. for (i = 0; i < ad->cl; i++)
  162. /* for each cell in the row */
  163. {
  164. corrCell = buf[i + ad->x];
  165. if (masked && mask_buf[i + ad->x] == 0) {
  166. Rast_set_c_null_value(&corrCell, 1);
  167. }
  168. if (!(Rast_is_null_value(&corrCell, gc.t))) {
  169. areaPatch++;
  170. if (i > 0)
  171. precCell = buf[i - 1 + ad->x];
  172. if (j == 0)
  173. Rast_set_c_null_value(&supCell, 1);
  174. else
  175. supCell = buf_sup[i + ad->x];
  176. if (corrCell != precCell)
  177. /* ?
  178. * 1 2
  179. * */
  180. {
  181. if (corrCell != supCell) {
  182. /* 3
  183. * 1 2
  184. * */
  185. /*new patch */
  186. if (idCorr == 0)
  187. /*first found patch */
  188. {
  189. lastId = 1;
  190. idCorr = 1;
  191. totCorr = 1;
  192. mask_patch_corr[i] = idCorr;
  193. }
  194. else
  195. /*not first patch */
  196. /* put in the tree the previous value */
  197. {
  198. if (albero == NULL) {
  199. albero = avlID_make(idCorr, totCorr);
  200. if (albero == NULL) {
  201. G_fatal_error("avlID_make error");
  202. return RLI_ERRORE;
  203. }
  204. npatch++;
  205. }
  206. else
  207. /*tree not empty */
  208. {
  209. ris = avlID_add(&albero, idCorr, totCorr);
  210. switch (ris) {
  211. case AVL_ERR:
  212. {
  213. G_fatal_error("avlID_add error");
  214. return RLI_ERRORE;
  215. }
  216. case AVL_ADD:
  217. {
  218. npatch++;
  219. break;
  220. }
  221. case AVL_PRES:
  222. {
  223. break;
  224. }
  225. default:
  226. {
  227. G_fatal_error
  228. ("avlID_add unknown error");
  229. return RLI_ERRORE;
  230. }
  231. }
  232. }
  233. totCorr = 1;
  234. lastId++;
  235. idCorr = lastId;
  236. mask_patch_corr[i] = idCorr;
  237. }
  238. }
  239. else
  240. /* current cell and upper cell are equal */
  241. /* 2
  242. * 1 2
  243. * */
  244. {
  245. if (albero == NULL) {
  246. albero = avlID_make(idCorr, totCorr);
  247. if (albero == NULL) {
  248. G_fatal_error("avlID_make error");
  249. return RLI_ERRORE;
  250. }
  251. npatch++;
  252. }
  253. else { /*tree not null */
  254. ris = avlID_add(&albero, idCorr, totCorr);
  255. switch (ris) {
  256. case AVL_ERR:
  257. {
  258. G_fatal_error("avlID_add error");
  259. return RLI_ERRORE;
  260. }
  261. case AVL_ADD:
  262. {
  263. npatch++;
  264. break;
  265. }
  266. case AVL_PRES:
  267. {
  268. break;
  269. }
  270. default:
  271. {
  272. G_fatal_error("avlID_add unknown error");
  273. return RLI_ERRORE;
  274. }
  275. }
  276. }
  277. idCorr = mask_patch_sup[i];
  278. mask_patch_corr[i] = idCorr;
  279. totCorr = 1;
  280. }
  281. }
  282. else { /*current cell and previuos cell are equal */
  283. /* ?
  284. * 1 1
  285. */
  286. if (corrCell == supCell) { /*current cell and upper cell are equal */
  287. /* 1
  288. * 1 1
  289. */
  290. if (mask_patch_sup[i] != mask_patch_corr[i - 1]) {
  291. long r = 0;
  292. long del = mask_patch_sup[i];
  293. r = avlID_sub(&albero, del); /*r=number of cell of patch removed */
  294. if (r == 0) {
  295. G_fatal_error("avlID_sub error");
  296. return RLI_ERRORE;
  297. }
  298. /*Remove one patch because it makes part of a patch already found */
  299. ris = avlID_add(&albero, idCorr, r);
  300. switch (ris) {
  301. case AVL_ERR:
  302. {
  303. G_fatal_error("avlID_add error");
  304. return RLI_ERRORE;
  305. }
  306. case AVL_ADD:
  307. {
  308. npatch++;
  309. break;
  310. }
  311. case AVL_PRES:
  312. {
  313. break;
  314. }
  315. default:
  316. {
  317. G_fatal_error("avlID_add unknown error");
  318. return RLI_ERRORE;
  319. }
  320. }
  321. r = i;
  322. while (r < ad->cl) {
  323. if (mask_patch_sup[r] == del) {
  324. mask_patch_sup[r] = idCorr;
  325. }
  326. r++;
  327. }
  328. mask_patch_corr[i] = idCorr;
  329. }
  330. else {
  331. mask_patch_corr[i] = idCorr;
  332. }
  333. }
  334. else { /*current cell and upper cell are not equal */
  335. /* 2
  336. * 1 1
  337. */
  338. mask_patch_corr[i] = idCorr;
  339. }
  340. totCorr++;
  341. }
  342. }
  343. else { /*cell is null or is not to consider */
  344. mask_patch_corr[i] = 0;
  345. }
  346. }
  347. {
  348. int ii;
  349. long c;
  350. for (ii = 0; ii < ad->cl; ii++) {
  351. c = mask_patch_corr[ii];
  352. mask_patch_sup[ii] = c;
  353. mask_patch_corr[ii] = 0;
  354. }
  355. }
  356. }
  357. if (areaPatch != 0) {
  358. if (albero == NULL) {
  359. albero = avlID_make(idCorr, totCorr);
  360. if (albero == NULL) {
  361. G_fatal_error("avlID_make error");
  362. return RLI_ERRORE;
  363. }
  364. npatch++;
  365. }
  366. else {
  367. ris = avlID_add(&albero, idCorr, totCorr);
  368. switch (ris) {
  369. case AVL_ERR:
  370. {
  371. G_fatal_error("avlID_add error");
  372. return RLI_ERRORE;
  373. }
  374. case AVL_ADD:
  375. {
  376. npatch++;
  377. break;
  378. }
  379. case AVL_PRES:
  380. {
  381. break;
  382. }
  383. default:
  384. {
  385. G_fatal_error("avlID_add unknown error");
  386. return RLI_ERRORE;
  387. }
  388. }
  389. }
  390. array = G_malloc(npatch * sizeof(avlID_tableRow));
  391. if (array == NULL) {
  392. G_fatal_error("malloc array failed");
  393. return RLI_ERRORE;
  394. }
  395. tot = avlID_to_array(albero, zero, array);
  396. if (tot != npatch) {
  397. G_warning
  398. ("avlID_to_array unaspected value. the result could be wrong");
  399. return RLI_ERRORE;
  400. }
  401. for (i = 0; i < npatch; i++) {
  402. if (array[i]->tot != 0) {
  403. if (array[i]->tot > area_max) {
  404. area_max = array[i]->tot;
  405. }
  406. if (array[i]->tot < area_min || area_min == 0) {
  407. area_min = array[i]->tot;
  408. }
  409. }
  410. }
  411. rk = area_max - area_min;
  412. indice = rk;
  413. G_free(array);
  414. }
  415. else
  416. indice = (double)(-1);
  417. if (masked)
  418. G_free(mask_buf);
  419. G_free(mask_patch_sup);
  420. *result = indice;
  421. G_free(buf_sup);
  422. return RLI_OK;
  423. }
  424. int calculateD(int fd, struct area_entry *ad, double *result)
  425. {
  426. DCELL *buf;
  427. DCELL *buf_sup;
  428. DCELL corrCell;
  429. DCELL precCell;
  430. DCELL supCell;
  431. int i, j;
  432. int mask_fd = -1, *mask_buf;
  433. int ris = 0;
  434. int masked = FALSE;
  435. int areaPatch = 0; /*if all cells are null areaPatch=0 */
  436. int area_max = 0, area_min = 0;
  437. long npatch = 0;
  438. long tot = 0;
  439. long zero = 0;
  440. long totCorr = 0;
  441. long idCorr = 0;
  442. long lastId = 0;
  443. long *mask_patch_sup;
  444. long *mask_patch_corr;
  445. double indice = 0;
  446. double rk = 0;
  447. avlID_tree albero = NULL;
  448. avlID_table *array = NULL;
  449. generic_cell gc;
  450. gc.t = DCELL_TYPE;
  451. /* open mask if needed */
  452. if (ad->mask == 1) {
  453. if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
  454. return RLI_ERRORE;
  455. mask_buf = G_malloc(ad->cl * sizeof(int));
  456. if (mask_buf == NULL) {
  457. G_fatal_error("malloc mask_buf failed");
  458. return RLI_ERRORE;
  459. }
  460. masked = TRUE;
  461. }
  462. mask_patch_sup = G_malloc(ad->cl * sizeof(long));
  463. if (mask_patch_sup == NULL) {
  464. G_fatal_error("malloc mask_patch_sup failed");
  465. return RLI_ERRORE;
  466. }
  467. mask_patch_corr = G_malloc(ad->cl * sizeof(long));
  468. if (mask_patch_corr == NULL) {
  469. G_fatal_error("malloc mask_patch_corr failed");
  470. return RLI_ERRORE;
  471. }
  472. buf_sup = Rast_allocate_d_buf();
  473. if (buf_sup == NULL) {
  474. G_fatal_error("malloc buf_sup failed");
  475. return RLI_ERRORE;
  476. }
  477. buf = Rast_allocate_d_buf();
  478. if (buf == NULL) {
  479. G_fatal_error("malloc buf failed");
  480. return RLI_ERRORE;
  481. }
  482. Rast_set_d_null_value(buf_sup + ad->x, ad->cl); /*the first time buf_sup is all null */
  483. for (i = 0; i < ad->cl; i++) {
  484. mask_patch_sup[i] = 0;
  485. mask_patch_corr[i] = 0;
  486. }
  487. for (j = 0; j < ad->rl; j++)
  488. /*for each raster row */
  489. {
  490. if (j > 0) {
  491. buf_sup = RLI_get_dcell_raster_row(fd, j - 1 + ad->y, ad);
  492. }
  493. buf = RLI_get_dcell_raster_row(fd, j + ad->y, ad);
  494. if (masked) {
  495. if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0) {
  496. G_fatal_error("mask read failed");
  497. return RLI_ERRORE;
  498. }
  499. }
  500. Rast_set_d_null_value(&precCell, 1);
  501. for (i = 0; i < ad->cl; i++) {
  502. /* for each dcell in the row */
  503. corrCell = buf[i + ad->x];
  504. if (masked && mask_buf[i + ad->x] == 0) {
  505. Rast_set_d_null_value(&corrCell, 1);
  506. }
  507. if (!(Rast_is_null_value(&corrCell, gc.t))) {
  508. areaPatch++;
  509. if (i > 0)
  510. precCell = buf[i - 1 + ad->x];
  511. if (j == 0)
  512. Rast_set_d_null_value(&supCell, 1);
  513. else
  514. supCell = buf_sup[i + ad->x];
  515. if (corrCell != precCell)
  516. /* ?
  517. * 1 2
  518. * */
  519. {
  520. if (corrCell != supCell) {
  521. /* 3
  522. * 1 2
  523. * */
  524. /*new patch */
  525. if (idCorr == 0) { /*first found patch */
  526. lastId = 1;
  527. idCorr = 1;
  528. totCorr = 1;
  529. mask_patch_corr[i] = idCorr;
  530. }
  531. else
  532. /*not first patch */
  533. /* put in the tree the previous value */
  534. {
  535. if (albero == NULL) {
  536. albero = avlID_make(idCorr, totCorr);
  537. if (albero == NULL) {
  538. G_fatal_error("avlID_make error");
  539. return RLI_ERRORE;
  540. }
  541. npatch++;
  542. }
  543. else
  544. /*tree not empty */
  545. {
  546. ris = avlID_add(&albero, idCorr, totCorr);
  547. switch (ris) {
  548. case AVL_ERR:
  549. {
  550. G_fatal_error("avlID_add error");
  551. return RLI_ERRORE;
  552. }
  553. case AVL_ADD:
  554. {
  555. npatch++;
  556. break;
  557. }
  558. case AVL_PRES:
  559. {
  560. break;
  561. }
  562. default:
  563. {
  564. G_fatal_error
  565. ("avlID_add unknown error");
  566. return RLI_ERRORE;
  567. }
  568. }
  569. }
  570. totCorr = 1;
  571. lastId++;
  572. idCorr = lastId;
  573. mask_patch_corr[i] = idCorr;
  574. }
  575. }
  576. else
  577. /* current cell and upper cell are equal */
  578. /* 2
  579. * 1 2
  580. * */
  581. {
  582. if (albero == NULL) {
  583. albero = avlID_make(idCorr, totCorr);
  584. if (albero == NULL) {
  585. G_fatal_error("avlID_make error");
  586. return RLI_ERRORE;
  587. }
  588. npatch++;
  589. }
  590. else { /*tree not null */
  591. ris = avlID_add(&albero, idCorr, totCorr);
  592. switch (ris) {
  593. case AVL_ERR:
  594. {
  595. G_fatal_error("avlID_add error");
  596. return RLI_ERRORE;
  597. }
  598. case AVL_ADD:
  599. {
  600. npatch++;
  601. break;
  602. }
  603. case AVL_PRES:
  604. {
  605. break;
  606. }
  607. default:
  608. {
  609. G_fatal_error("avlID_add unknown error");
  610. return RLI_ERRORE;
  611. }
  612. }
  613. }
  614. idCorr = mask_patch_sup[i];
  615. mask_patch_corr[i] = idCorr;
  616. totCorr = 1;
  617. }
  618. }
  619. else { /*current cell and previuos cell are equal */
  620. /* ?
  621. * 1 1
  622. */
  623. if (corrCell == supCell) { /*current cell and upper cell are equal */
  624. /* 1
  625. * 1 1
  626. */
  627. if (mask_patch_sup[i] != mask_patch_corr[i - 1]) {
  628. long r = 0;
  629. long del = mask_patch_sup[i];
  630. r = avlID_sub(&albero, del); /*r=number of cell of patch removed */
  631. if (r == 0) {
  632. G_fatal_error("avlID_sub error");
  633. return RLI_ERRORE;
  634. }
  635. /*Remove one patch because it makes part of a patch already found */
  636. ris = avlID_add(&albero, idCorr, r);
  637. switch (ris) {
  638. case AVL_ERR:
  639. {
  640. G_fatal_error("avlID_add error");
  641. return RLI_ERRORE;
  642. }
  643. case AVL_ADD:
  644. {
  645. npatch++;
  646. break;
  647. }
  648. case AVL_PRES:
  649. {
  650. break;
  651. }
  652. default:
  653. {
  654. G_fatal_error("avlID_add unknown error");
  655. return RLI_ERRORE;
  656. }
  657. }
  658. r = i;
  659. while (r < ad->cl) {
  660. if (mask_patch_sup[r] == del) {
  661. mask_patch_sup[r] = idCorr;
  662. }
  663. r++;
  664. }
  665. mask_patch_corr[i] = idCorr;
  666. }
  667. else {
  668. mask_patch_corr[i] = idCorr;
  669. }
  670. }
  671. else { /*current cell and upper cell are not equal */
  672. /* 2
  673. * 1 1
  674. */
  675. mask_patch_corr[i] = idCorr;
  676. }
  677. totCorr++;
  678. }
  679. }
  680. else { /*cell is null or is not to consider */
  681. mask_patch_corr[i] = 0;
  682. }
  683. }
  684. {
  685. int ii;
  686. long c;
  687. for (ii = 0; ii < ad->cl; ii++) {
  688. c = mask_patch_corr[ii];
  689. mask_patch_sup[ii] = c;
  690. mask_patch_corr[ii] = 0;
  691. }
  692. }
  693. }
  694. if (areaPatch != 0) {
  695. if (albero == NULL) {
  696. albero = avlID_make(idCorr, totCorr);
  697. if (albero == NULL) {
  698. G_fatal_error("avlID_make error");
  699. return RLI_ERRORE;
  700. }
  701. npatch++;
  702. }
  703. else {
  704. ris = avlID_add(&albero, idCorr, totCorr);
  705. switch (ris) {
  706. case AVL_ERR:
  707. {
  708. G_fatal_error("avlID_add error");
  709. return RLI_ERRORE;
  710. }
  711. case AVL_ADD:
  712. {
  713. npatch++;
  714. break;
  715. }
  716. case AVL_PRES:
  717. {
  718. break;
  719. }
  720. default:
  721. {
  722. G_fatal_error("avlID_add unknown error");
  723. return RLI_ERRORE;
  724. }
  725. }
  726. }
  727. array = G_malloc(npatch * sizeof(avlID_tableRow));
  728. if (array == NULL) {
  729. G_fatal_error("malloc array failed");
  730. return RLI_ERRORE;
  731. }
  732. tot = avlID_to_array(albero, zero, array);
  733. if (tot != npatch) {
  734. G_warning
  735. ("avlID_to_array unaspected value. the result could be wrong");
  736. return RLI_ERRORE;
  737. }
  738. for (i = 0; i < npatch; i++) {
  739. if (array[i]->tot != 0) {
  740. if (array[i]->tot > area_max) {
  741. area_max = array[i]->tot;
  742. }
  743. if (array[i]->tot < area_min || area_min == 0) {
  744. area_min = array[i]->tot;
  745. }
  746. }
  747. }
  748. rk = area_max - area_min;
  749. indice = rk;
  750. G_free(array);
  751. }
  752. else
  753. indice = (double)(-1);
  754. if (masked)
  755. G_free(mask_buf);
  756. G_free(mask_patch_sup);
  757. *result = indice;
  758. return RLI_OK;
  759. }
  760. int calculateF(int fd, struct area_entry *ad, double *result)
  761. {
  762. FCELL *buf;
  763. FCELL *buf_sup;
  764. FCELL corrCell;
  765. FCELL precCell;
  766. FCELL supCell;
  767. int i, j;
  768. int mask_fd = -1, *mask_buf;
  769. int ris = 0;
  770. int masked = FALSE;
  771. int areaPatch = 0; /*if all cells are null areaPatch=0 */
  772. int area_max = 0, area_min = 0;
  773. long npatch = 0;
  774. long tot = 0;
  775. long zero = 0;
  776. long totCorr = 0;
  777. long idCorr = 0;
  778. long lastId = 0;
  779. long *mask_patch_sup;
  780. long *mask_patch_corr;
  781. double indice = 0;
  782. double rk = 0;
  783. avlID_tree albero = NULL;
  784. avlID_table *array = NULL;
  785. generic_cell gc;
  786. gc.t = FCELL_TYPE;
  787. /* open mask if needed */
  788. if (ad->mask == 1) {
  789. if ((mask_fd = open(ad->mask_name, O_RDONLY, 0755)) < 0)
  790. return RLI_ERRORE;
  791. mask_buf = G_malloc(ad->cl * sizeof(int));
  792. if (mask_buf == NULL) {
  793. G_fatal_error("malloc mask_buf failed");
  794. return RLI_ERRORE;
  795. }
  796. masked = TRUE;
  797. }
  798. mask_patch_sup = G_malloc(ad->cl * sizeof(long));
  799. if (mask_patch_sup == NULL) {
  800. G_fatal_error("malloc mask_patch_sup failed");
  801. return RLI_ERRORE;
  802. }
  803. mask_patch_corr = G_malloc(ad->cl * sizeof(long));
  804. if (mask_patch_corr == NULL) {
  805. G_fatal_error("malloc mask_patch_corr failed");
  806. return RLI_ERRORE;
  807. }
  808. buf_sup = Rast_allocate_f_buf();
  809. if (buf_sup == NULL) {
  810. G_fatal_error("malloc buf_sup failed");
  811. return RLI_ERRORE;
  812. }
  813. buf = Rast_allocate_f_buf();
  814. if (buf == NULL) {
  815. G_fatal_error("malloc buf failed");
  816. return RLI_ERRORE;
  817. }
  818. Rast_set_f_null_value(buf_sup + ad->x, ad->cl); /*the first time buf_sup is all null */
  819. for (i = 0; i < ad->cl; i++) {
  820. mask_patch_sup[i] = 0;
  821. mask_patch_corr[i] = 0;
  822. }
  823. for (j = 0; j < ad->rl; j++)
  824. /*for each raster row */
  825. {
  826. if (j > 0) {
  827. buf_sup = RLI_get_fcell_raster_row(fd, j - 1 + ad->y, ad);
  828. }
  829. buf = RLI_get_fcell_raster_row(fd, j + ad->y, ad);
  830. if (masked) {
  831. if (read(mask_fd, mask_buf, (ad->cl * sizeof(int))) < 0) {
  832. G_fatal_error("mask read failed");
  833. return RLI_ERRORE;
  834. }
  835. }
  836. Rast_set_f_null_value(&precCell, 1);
  837. for (i = 0; i < ad->cl; i++) {
  838. /* for each fcell in the row */
  839. corrCell = buf[i + ad->x];
  840. if (masked && mask_buf[i + ad->x] == 0) {
  841. Rast_set_f_null_value(&corrCell, 1);
  842. }
  843. if (!(Rast_is_null_value(&corrCell, gc.t))) {
  844. areaPatch++;
  845. if (i > 0)
  846. precCell = buf[i - 1 + ad->x];
  847. if (j == 0)
  848. Rast_set_f_null_value(&supCell, 1);
  849. else
  850. supCell = buf_sup[i + ad->x];
  851. if (corrCell != precCell)
  852. /* ?
  853. * 1 2
  854. * */
  855. {
  856. if (corrCell != supCell) {
  857. /* 3
  858. * 1 2
  859. * */
  860. /*new patch */
  861. if (idCorr == 0) { /*first found patch */
  862. lastId = 1;
  863. idCorr = 1;
  864. totCorr = 1;
  865. mask_patch_corr[i] = idCorr;
  866. }
  867. else
  868. /*not first patch */
  869. /* put in the tree the previous value */
  870. {
  871. if (albero == NULL) {
  872. albero = avlID_make(idCorr, totCorr);
  873. if (albero == NULL) {
  874. G_fatal_error("avlID_make error");
  875. return RLI_ERRORE;
  876. }
  877. npatch++;
  878. }
  879. else
  880. /*tree not empty */
  881. {
  882. ris = avlID_add(&albero, idCorr, totCorr);
  883. switch (ris) {
  884. case AVL_ERR:
  885. {
  886. G_fatal_error("avlID_add error");
  887. return RLI_ERRORE;
  888. }
  889. case AVL_ADD:
  890. {
  891. npatch++;
  892. break;
  893. }
  894. case AVL_PRES:
  895. {
  896. break;
  897. }
  898. default:
  899. {
  900. G_fatal_error
  901. ("avlID_add unknown error");
  902. return RLI_ERRORE;
  903. }
  904. }
  905. }
  906. totCorr = 1;
  907. lastId++;
  908. idCorr = lastId;
  909. mask_patch_corr[i] = idCorr;
  910. }
  911. }
  912. else
  913. /* current cell and upper cell are equal */
  914. /* 2
  915. * 1 2
  916. * */
  917. {
  918. if (albero == NULL) {
  919. albero = avlID_make(idCorr, totCorr);
  920. if (albero == NULL) {
  921. G_fatal_error("avlID_make error");
  922. return RLI_ERRORE;
  923. }
  924. npatch++;
  925. }
  926. else { /*tree not null */
  927. ris = avlID_add(&albero, idCorr, totCorr);
  928. switch (ris) {
  929. case AVL_ERR:
  930. {
  931. G_fatal_error("avlID_add error");
  932. return RLI_ERRORE;
  933. }
  934. case AVL_ADD:
  935. {
  936. npatch++;
  937. break;
  938. }
  939. case AVL_PRES:
  940. {
  941. break;
  942. }
  943. default:
  944. {
  945. G_fatal_error("avlID_add unknown error");
  946. return RLI_ERRORE;
  947. }
  948. }
  949. }
  950. idCorr = mask_patch_sup[i];
  951. mask_patch_corr[i] = idCorr;
  952. totCorr = 1;
  953. }
  954. }
  955. else { /*current cell and previuos cell are equal */
  956. /* ?
  957. * 1 1
  958. */
  959. if (corrCell == supCell) { /*current cell and upper cell are equal */
  960. /* 1
  961. * 1 1
  962. */
  963. if (mask_patch_sup[i] != mask_patch_corr[i - 1]) {
  964. long r = 0;
  965. long del = mask_patch_sup[i];
  966. r = avlID_sub(&albero, del); /*r=number of cell of patch removed */
  967. if (r == 0) {
  968. G_fatal_error("avlID_sub error");
  969. return RLI_ERRORE;
  970. }
  971. /*Remove one patch because it makes part of a patch already found */
  972. ris = avlID_add(&albero, idCorr, r);
  973. switch (ris) {
  974. case AVL_ERR:
  975. {
  976. G_fatal_error("avlID_add error");
  977. return RLI_ERRORE;
  978. }
  979. case AVL_ADD:
  980. {
  981. npatch++;
  982. break;
  983. }
  984. case AVL_PRES:
  985. {
  986. break;
  987. }
  988. default:
  989. {
  990. G_fatal_error("avlID_add unknown error");
  991. return RLI_ERRORE;
  992. }
  993. }
  994. r = i;
  995. while (r < ad->cl) {
  996. if (mask_patch_sup[r] == del) {
  997. mask_patch_sup[r] = idCorr;
  998. }
  999. r++;
  1000. }
  1001. mask_patch_corr[i] = idCorr;
  1002. }
  1003. else {
  1004. mask_patch_corr[i] = idCorr;
  1005. }
  1006. }
  1007. else { /*current cell and upper cell are not equal */
  1008. /* 2
  1009. * 1 1
  1010. */
  1011. mask_patch_corr[i] = idCorr;
  1012. }
  1013. totCorr++;
  1014. }
  1015. }
  1016. else { /*cell is null or is not to consider */
  1017. mask_patch_corr[i] = 0;
  1018. }
  1019. }
  1020. {
  1021. int ii;
  1022. long c;
  1023. for (ii = 0; ii < ad->cl; ii++) {
  1024. c = mask_patch_corr[ii];
  1025. mask_patch_sup[ii] = c;
  1026. mask_patch_corr[ii] = 0;
  1027. }
  1028. }
  1029. }
  1030. if (areaPatch != 0) {
  1031. if (albero == NULL) {
  1032. albero = avlID_make(idCorr, totCorr);
  1033. if (albero == NULL) {
  1034. G_fatal_error("avlID_make error");
  1035. return RLI_ERRORE;
  1036. }
  1037. npatch++;
  1038. }
  1039. else {
  1040. ris = avlID_add(&albero, idCorr, totCorr);
  1041. switch (ris) {
  1042. case AVL_ERR:
  1043. {
  1044. G_fatal_error("avlID_add error");
  1045. return RLI_ERRORE;
  1046. }
  1047. case AVL_ADD:
  1048. {
  1049. npatch++;
  1050. break;
  1051. }
  1052. case AVL_PRES:
  1053. {
  1054. break;
  1055. }
  1056. default:
  1057. {
  1058. G_fatal_error("avlID_add unknown error");
  1059. return RLI_ERRORE;
  1060. }
  1061. }
  1062. }
  1063. array = G_malloc(npatch * sizeof(avlID_tableRow));
  1064. if (array == NULL) {
  1065. G_fatal_error("malloc array failed");
  1066. return RLI_ERRORE;
  1067. }
  1068. tot = avlID_to_array(albero, zero, array);
  1069. if (tot != npatch) {
  1070. G_warning
  1071. ("avlID_to_array unaspected value. the result could be wrong");
  1072. return RLI_ERRORE;
  1073. }
  1074. for (i = 0; i < npatch; i++) {
  1075. if (array[i]->tot != 0) {
  1076. if (array[i]->tot > area_max) {
  1077. area_max = array[i]->tot;
  1078. }
  1079. if (array[i]->tot < area_min || area_min == 0) {
  1080. area_min = array[i]->tot;
  1081. }
  1082. }
  1083. }
  1084. rk = area_max - area_min;
  1085. indice = rk;
  1086. G_free(array);
  1087. }
  1088. else
  1089. indice = (double)(-1);
  1090. if (masked)
  1091. G_free(mask_buf);
  1092. G_free(mask_patch_sup);
  1093. *result = indice;
  1094. return RLI_OK;
  1095. }