init.c 8.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303
  1. /* init.c */
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include <string.h>
  5. #include <unistd.h>
  6. #include <sys/types.h>
  7. #include <grass/gis.h>
  8. #include <grass/raster.h>
  9. #include <grass/glocale.h>
  10. #include "ransurf.h"
  11. #include "local_proto.h"
  12. void Init(void)
  13. {
  14. struct Cell_head Region;
  15. int row, col, i, j, NumWeight, NumDist, NumExp;
  16. char *Name, *Number, String[80];
  17. double MinRes;
  18. G_debug(2, "Init");
  19. Rs = Rast_window_rows();
  20. Cs = Rast_window_cols();
  21. RSurface = (double **)G_malloc(Rs * sizeof(double *));
  22. for (row = 0; row < Rs; row++)
  23. RSurface[row] = (double *)G_malloc(Cs * sizeof(double));
  24. G_get_set_window(&Region);
  25. EW = Region.ew_res;
  26. NS = Region.ns_res;
  27. if (EW < NS)
  28. MinRes = EW;
  29. else
  30. MinRes = NS;
  31. if (NULL == G_find_file("cell", "MASK", G_mapset())) {
  32. MapCount = Rs * Cs;
  33. FDM = -1;
  34. }
  35. else {
  36. FDM = Rast_open_old("MASK", G_mapset());
  37. {
  38. MapCount = 0;
  39. CellBuffer = Rast_allocate_c_buf();
  40. for (row = 0; row < Rs; row++) {
  41. Rast_get_c_row_nomask(FDM, CellBuffer, row);
  42. for (col = 0; col < Cs; col++) {
  43. if (CellBuffer[col])
  44. MapCount++;
  45. }
  46. }
  47. }
  48. }
  49. if (Uniform->answer)
  50. sprintf(Buf, "Uni. R. S.");
  51. else
  52. sprintf(Buf, "Dist. R. S.");
  53. if (!range_high_stuff->answer)
  54. High = 255;
  55. else {
  56. High = atoi(range_high_stuff->answer);
  57. sprintf(String, " high=%d", High);
  58. strcat(Buf, String);
  59. }
  60. if (1 >= High)
  61. G_fatal_error(_("High (%d) must be greater than 1"), High);
  62. CatInfo.NumCat = High;
  63. NumMaps = 0;
  64. OutNames = (char **)G_malloc(sizeof(char *));
  65. for (i = 0; (Name = Output->answers[i]); i++) {
  66. for (j = i - 1; j >= 0; j--) {
  67. if (strcmp(OutNames[j], Name) == 0)
  68. G_fatal_error
  69. (_("Rastar map <%s> repeated, maps must be unique"),
  70. Name);
  71. }
  72. OutNames = (char **)G_realloc(OutNames, (i + 1) * sizeof(char *));
  73. OutNames[i] = (char *)G_malloc(sizeof(char) * (strlen(Name) + 1));
  74. strcpy(OutNames[i], Name);
  75. NumMaps++;
  76. }
  77. if (NumMaps == 0)
  78. G_fatal_error(_("Output raster map required"));
  79. Theory = 0;
  80. NumSeeds = 0;
  81. Seeds = (int *)G_malloc(NumMaps * sizeof(int));
  82. Seed = -1;
  83. if (!SeedStuff->answers) {
  84. for (i = 0; i < NumMaps; i++) {
  85. Seeds[i] = -1;
  86. }
  87. }
  88. else {
  89. for (i = 0; (Number = SeedStuff->answers[i]) && i < NumMaps; i++) {
  90. sscanf(Number, "%d", &(Seeds[i]));
  91. } /* /for */
  92. } /* /else */
  93. CellBuffer = Rast_allocate_c_buf();
  94. CatInfo.NumValue = (int *)G_malloc(CatInfo.NumCat * sizeof(int));
  95. CatInfo.Average = (double *)G_malloc(CatInfo.NumCat * sizeof(double));
  96. CatInfo.Min = (double *)G_malloc(CatInfo.NumCat * sizeof(double));
  97. CatInfo.Max = (double *)G_malloc(CatInfo.NumCat * sizeof(double));
  98. NumFilters = 1;
  99. AllFilters = (FILTER *) G_malloc(sizeof(FILTER));
  100. /*
  101. if( ! Distance->answers) {
  102. NumDist = 1;
  103. AllFilters[ 0].MaxDist = MinRes / 2.0;
  104. } else {
  105. NumDist = 0;
  106. for (i = 0; Number = Distance->answers[i]; i++) {
  107. if(NumDist == NumFilters) {
  108. AllFilters = (FILTER *) G_realloc( AllFilters,
  109. ++NumFilters * sizeof( FILTER));
  110. }
  111. sscanf( Number, "%lf", &(AllFilters[NumDist].MaxDist));
  112. if (AllFilters[NumDist].MaxDist < 0.0)
  113. G_fatal_error("%s: distance value[%d]: [%lf] must be >= 0.0",
  114. G_program_name(), NumDist,
  115. AllFilters[NumDist].MaxDist);
  116. NumDist++;
  117. }
  118. }
  119. */
  120. NumDist = 0;
  121. if (Distance->answer) {
  122. sscanf(Distance->answer, "%lf", &(AllFilters[NumDist].MaxDist));
  123. if (AllFilters[NumDist].MaxDist < 0.0)
  124. G_fatal_error(_("Distance value (%d): %lf must be >= 0.0"),
  125. NumDist, AllFilters[NumDist].MaxDist);
  126. NumDist++;
  127. }
  128. /*
  129. NumExp = 0;
  130. if( Exponent->answers) {
  131. for (i = 0; Number = Exponent->answers[i]; i++) {
  132. if(NumExp == NumFilters)
  133. AllFilters = (FILTER *) G_realloc(
  134. AllFilters, ++NumFilters * sizeof( FILTER));
  135. sscanf( Number, "%lf", &(AllFilters[NumExp].Exp));
  136. DOUBLE(AllFilters[NumExp].Exp);
  137. if (AllFilters[NumExp].Exp <= 0.0)
  138. G_fatal_error("%s: exponent value [%lf] must be > 0.0",
  139. G_program_name(), AllFilters[NumExp].Exp);
  140. AllFilters[ NumExp].Exp = 1.0 / AllFilters[ NumExp].Exp;
  141. NumExp++;
  142. }
  143. }
  144. INT(NumExp);
  145. */
  146. NumExp = 0;
  147. if (Exponent->answer) {
  148. sscanf(Exponent->answer, "%lf", &(AllFilters[NumExp].Exp));
  149. if (AllFilters[NumExp].Exp <= 0.0)
  150. G_fatal_error(_("Exponent value (%lf) must be > 0.0"),
  151. AllFilters[NumExp].Exp);
  152. NumExp++;
  153. }
  154. NumWeight = 0;
  155. /*
  156. if( Weight->answers) {
  157. for (i = 0; Number = Weight->answers[i]; i++) {
  158. if(NumWeight == NumFilters)
  159. AllFilters = (FILTER *) G_realloc(
  160. AllFilters, ++NumFilters * sizeof( FILTER));
  161. sscanf( Numbers, "%lf", &(AllFilters[NumWeight].Mult));
  162. if (AllFilters[NumWeight].Mult == 0.0)
  163. G_fatal_error("%s: weight value [%lf] must not be 0.0",
  164. G_program_name(), AllFilters[NumWeight].Mult);
  165. NumWeight++;
  166. }
  167. }
  168. */
  169. if (Weight->answer) {
  170. sscanf(Weight->answer, "%lf", &(AllFilters[NumWeight].Mult));
  171. if (AllFilters[NumWeight].Mult > AllFilters[NumWeight].MaxDist)
  172. G_fatal_error(_("Flat value (%lf) must be less than distance value (%lf)"),
  173. AllFilters[NumWeight].Mult,
  174. AllFilters[NumWeight].MaxDist);
  175. NumWeight++;
  176. }
  177. if (NumDist > 0) {
  178. sprintf(String, " dist=");
  179. strcat(Buf, String);
  180. for (i = 0; i < NumDist - 1; i++) {
  181. sprintf(String, "%.*lf,",
  182. Digits(AllFilters[i].MaxDist, 6), AllFilters[i].MaxDist);
  183. strcat(Buf, String);
  184. }
  185. sprintf(String, "%.*lf",
  186. Digits(AllFilters[i].MaxDist, 6), AllFilters[i].MaxDist);
  187. strcat(Buf, String);
  188. }
  189. if (NumDist > 1 && NumDist < NumFilters)
  190. G_fatal_error(_("Must have a distance value for each filter"));
  191. if (NumDist == 0) {
  192. AllFilters[0].MaxDist = MinRes / 4.0;
  193. }
  194. if (NumDist < NumFilters) {
  195. for (NumDist = 1; NumDist < NumFilters; NumDist++)
  196. AllFilters[NumDist].MaxDist = AllFilters[0].MaxDist;
  197. }
  198. for (NumDist = 0; NumDist < NumFilters; NumDist++) {
  199. if (AllFilters[NumDist].MaxDist < MinRes)
  200. AllFilters[NumDist].MaxDist = MinRes * 0.5;
  201. else
  202. AllFilters[NumDist].MaxDist *= 0.5;
  203. }
  204. if (NumExp > 0) {
  205. sprintf(String, " exp=");
  206. strcat(Buf, String);
  207. for (i = 0; i < NumExp - 1; i++) {
  208. sprintf(String, "%.*lf,",
  209. Digits(AllFilters[i].Exp, 6), AllFilters[i].Exp);
  210. strcat(Buf, String);
  211. }
  212. sprintf(String, "%.*lf",
  213. Digits(AllFilters[i].Exp, 6), AllFilters[i].Exp);
  214. strcat(Buf, String);
  215. }
  216. if (NumExp > 1 && NumExp < NumFilters)
  217. G_fatal_error(_("Must have a exponent value for each filter"));
  218. if (NumWeight > 0) {
  219. sprintf(String, " flat=");
  220. strcat(Buf, String);
  221. for (i = 0; i < NumWeight - 1; i++) {
  222. sprintf(String, "%.*lf,",
  223. Digits(AllFilters[i].Mult, 6), AllFilters[i].Mult);
  224. strcat(Buf, String);
  225. G_debug(3, "(AllFilters[i].Mult):%.12lf", AllFilters[i].Mult);
  226. }
  227. sprintf(String, "%.*lf",
  228. Digits(AllFilters[i].Mult, 6), AllFilters[i].Mult);
  229. strcat(Buf, String);
  230. }
  231. if (NumWeight > 1 && NumWeight < NumFilters)
  232. G_fatal_error(_("Must have a weight value for each filter"));
  233. if (NumExp == 1) {
  234. for (NumExp = 1; NumExp < NumFilters; NumExp++)
  235. AllFilters[NumExp].Exp = AllFilters[0].Exp;
  236. }
  237. if (NumExp == 0) {
  238. for (NumExp = 0; NumExp < NumFilters; NumExp++)
  239. AllFilters[NumExp].Exp = 1.0;
  240. }
  241. if (NumWeight == 0) {
  242. for (NumWeight = 0; NumWeight < NumFilters; NumWeight++)
  243. AllFilters[NumWeight].Mult = 0.0;
  244. }
  245. AllMaxDist = 0.0;
  246. for (i = 0; i < NumFilters; i++) {
  247. if (AllMaxDist < AllFilters[i].MaxDist)
  248. AllMaxDist = AllFilters[i].MaxDist;
  249. AllFilters[i].MaxSq = AllFilters[i].MaxDist * AllFilters[i].MaxDist;
  250. G_debug(3, "(i):%d", i);
  251. G_debug(3, "(AllFilters[i].Mult):%.12lf", AllFilters[i].Mult);
  252. G_debug(3, "(AllFilters[i].MaxDist):%.12lf", AllFilters[i].MaxDist);
  253. G_debug(3, "(AllFilters[i].MaxSq):%.12lf", AllFilters[i].MaxSq);
  254. G_debug(3, "(AllFilters[i].Exp):%.12lf", AllFilters[i].Exp);
  255. }
  256. BigF.RowPlus = AllMaxDist / NS;
  257. BigF.ColPlus = AllMaxDist / EW;
  258. BigF.NumR = BigF.RowPlus * 2 + 1;
  259. BigF.NumC = BigF.ColPlus * 2 + 1;
  260. BigF.LowBF = (int *)G_malloc(BigF.NumR * sizeof(int));
  261. BigF.HihBF = (int *)G_malloc(BigF.NumR * sizeof(int));
  262. BigF.F = (double **)G_malloc(BigF.NumR * sizeof(double *));
  263. for (i = 0; i < BigF.NumR; i++)
  264. BigF.F[i] = (double *)G_malloc(BigF.NumC * sizeof(double));
  265. AllMaxDist *= 2.0;
  266. }