init.c 8.5 KB

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