GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-f241532c70
filecompare.c
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <sys/types.h>
4 #include <unistd.h>
5 #include <grass/raster3d.h>
6 
7 /*--------------------------------------------------------------------------*/
8 
9 static unsigned char clearMask[9] = {255, 128, 192, 224, 240,
10  248, 252, 254, 255};
11 
12 /*---------------------------------------------------------------------------*/
13 
14 static void Rast3d_float2xdrFloat(const float *f, float *xdrf)
15 {
16  G_xdr_put_float(xdrf, f);
17 }
18 
19 /*---------------------------------------------------------------------------*/
20 
21 static void Rast3d_double2xdrDouble(const double *d, double *xdrd)
22 {
23  G_xdr_put_double(xdrd, d);
24 }
25 
26 /*---------------------------------------------------------------------------*/
27 
28 static void Rast3d_truncFloat(float *f, int p)
29 {
30  unsigned char *c;
31 
32  if ((p == -1) || (p >= 23))
33  return;
34 
35  c = (unsigned char *)f;
36 
37  c++;
38  if (p <= 7) {
39  *c++ &= clearMask[(p + 1) % 8];
40  *c++ = 0;
41  *c = 0;
42  return;
43  }
44 
45  c++;
46  if (p <= 15) {
47  *c++ &= clearMask[(p + 1) % 8];
48  *c = 0;
49  return;
50  }
51 
52  c++;
53  *c &= clearMask[(p + 1) % 8];
54  return;
55 }
56 
57 /*---------------------------------------------------------------------------*/
58 
59 static void Rast3d_truncDouble(double *d, int p)
60 {
61  unsigned char *c;
62 
63  if ((p == -1) || (p >= 52))
64  return;
65 
66  c = (unsigned char *)d;
67 
68  c++;
69  if (p <= 4) {
70  *c++ &= clearMask[(p + 4) % 8];
71  *c++ = 0;
72  *c++ = 0;
73  *c++ = 0;
74  *c++ = 0;
75  *c++ = 0;
76  *c = 0;
77  return;
78  }
79 
80  c++;
81  if (p <= 12) {
82  *c++ &= clearMask[(p + 4) % 8];
83  *c++ = 0;
84  *c++ = 0;
85  *c++ = 0;
86  *c++ = 0;
87  *c = 0;
88  return;
89  }
90 
91  c++;
92  if (p <= 20) {
93  *c++ &= clearMask[(p + 4) % 8];
94  *c++ = 0;
95  *c++ = 0;
96  *c++ = 0;
97  *c = 0;
98  return;
99  }
100 
101  c++;
102  if (p <= 28) {
103  *c++ &= clearMask[(p + 4) % 8];
104  *c++ = 0;
105  *c++ = 0;
106  *c = 0;
107  return;
108  }
109 
110  c++;
111  if (p <= 36) {
112  *c++ &= clearMask[(p + 4) % 8];
113  *c++ = 0;
114  *c = 0;
115  return;
116  }
117 
118  c++;
119  if (p <= 44) {
120  *c++ &= clearMask[(p + 4) % 8];
121  *c = 0;
122  return;
123  }
124 
125  c++;
126  *c &= clearMask[(p + 4) % 8];
127  return;
128 }
129 
130 /*---------------------------------------------------------------------------*/
131 
132 static void Rast3d_float2Double(float *f, double *d)
133 {
134  unsigned char *c1, *c2, sign, c;
135  int e;
136 
137  c1 = (unsigned char *)f;
138  c2 = (unsigned char *)d;
139 
140  sign = (*c1 & (unsigned char)128);
141  e = (((*c1 & (unsigned char)127) << 1) |
142  ((*(c1 + 1) & (unsigned char)128) >> 7));
143 
144  if ((*c1 != 0) || (*(c1 + 1) != 0) || (*(c1 + 2) != 0) || (*(c1 + 3) != 0))
145  e += 1023 - 127;
146  c = e / 16;
147 
148  *c2++ = (sign | c);
149 
150  c1++;
151 
152  c = e % 16;
153  *c2 = (c << 4);
154  *c2++ |= ((*c1 & (unsigned char)127) >> 3);
155 
156  *c2 = ((*c1++ & (unsigned char)7) << 5);
157  *c2++ |= (*c1 >> 3);
158 
159  *c2 = ((*c1++ & (unsigned char)7) << 5);
160  *c2++ |= (*c1 >> 3);
161 
162  *c2++ = ((*c1 & (unsigned char)7) << 5);
163 
164  *c2++ = (unsigned char)0;
165  *c2++ = (unsigned char)0;
166  *c2 = (unsigned char)0;
167 }
168 
169 /*---------------------------------------------------------------------------*/
170 
171 static int Rast3d_compareFloats(float *f1, int p1, float *f2, int p2)
172 {
173  unsigned char *c1, *c2;
174  float xdrf1, xdrf2;
175 
178 
179  Rast3d_float2xdrFloat(f1, &xdrf1);
180  Rast3d_float2xdrFloat(f2, &xdrf2);
181 
182  c1 = (unsigned char *)&xdrf1;
183  c2 = (unsigned char *)&xdrf2;
184 
185  /* printf ("%d %d (%d %d %d %d) (%d %d %d %d) %d\n", p1, p2, *c1, *(c1 +
186  * 1), *(c1 + 2), *(c1 + 3), *c2, *(c2 + 1), *(c2 + 2), *(c2 + 3), *f1 ==
187  * *f2); */
188 
189  if ((p1 != -1) && (p1 < 23) && ((p1 < p2) || (p2 == -1)))
190  Rast3d_truncFloat(&xdrf2, p1);
191  if ((p2 != -1) && (p2 < 23) && ((p2 < p1) || (p1 == -1)))
192  Rast3d_truncFloat(&xdrf1, p2);
193 
194  /* printf ("%d %d (%d %d %d %d) (%d %d %d %d) %d\n", p1, p2, *c1, *(c1 +
195  * 1), *(c1 + 2), *(c1 + 3), *c2, *(c2 + 1), *(c2 + 2), *(c2 + 3), *f1 ==
196  * *f2); */
197 
198  return (*c1 == *c2) && (*(c1 + 1) == *(c2 + 1)) &&
199  (*(c1 + 2) == *(c2 + 2)) && (*(c1 + 3) == *(c2 + 3));
200 }
201 
202 /*---------------------------------------------------------------------------*/
203 
204 static int Rast3d_compareDoubles(double *d1, int p1, double *d2, int p2)
205 {
206  unsigned char *c1, *c2;
207  double xdrd1, xdrd2;
208 
211 
212  Rast3d_double2xdrDouble(d1, &xdrd1);
213  Rast3d_double2xdrDouble(d2, &xdrd2);
214 
215  c1 = (unsigned char *)&xdrd1;
216  c2 = (unsigned char *)&xdrd2;
217 
218  /* printf ("%d %d (%d %d %d %d %d %d %d %d) (%d %d %d %d %d %d %d %d)\n",
219  * p1, p2, *c1, *(c1 + 1), *(c1 + 2), *(c1 + 3), *(c1 + 4), *(c1 + 5), *(c1
220  * + 6), *(c1 + 7), *c2, *(c2 + 1), *(c2 + 2), *(c2 + 3), *(c2 + 4), *(c2 +
221  * 5), *(c2 + 6), *(c2 + 7)); */
222 
223  if ((p1 != -1) && (p1 < 52) && ((p1 < p2) || (p2 == -1)))
224  Rast3d_truncDouble(&xdrd2, p1);
225  if ((p2 != -1) && (p2 < 52) && ((p2 < p1) || (p1 == -1)))
226  Rast3d_truncDouble(&xdrd1, p2);
227 
228  /* printf ("%d %d (%d %d %d %d %d %d %d %d) (%d %d %d %d %d %d %d %d)\n",
229  * p1, p2, *c1, *(c1 + 1), *(c1 + 2), *(c1 + 3), *(c1 + 4), *(c1 + 5), *(c1
230  * + 6), *(c1 + 7), *c2, *(c2 + 1), *(c2 + 2), *(c2 + 3), *(c2 + 4), *(c2 +
231  * 5), *(c2 + 6), *(c2 + 7)); */
232 
233  return (*c1 == *c2) && (*(c1 + 1) == *(c2 + 1)) &&
234  (*(c1 + 2) == *(c2 + 2)) && (*(c1 + 3) == *(c2 + 3)) &&
235  (*(c1 + 4) == *(c2 + 4)) && (*(c1 + 5) == *(c2 + 5)) &&
236  (*(c1 + 6) == *(c2 + 6)) && (*(c1 + 7) == *(c2 + 7));
237 }
238 
239 /*---------------------------------------------------------------------------*/
240 
241 static int Rast3d_compareFloatDouble(float *f, int p1, double *d, int p2)
242 {
243  unsigned char *c1, *c2;
244  float xdrf, fTmp;
245  double xdrd, xdrd2, dTmp;
246 
249 
250  /* need this since assigning a double to a float actually may change the */
251  /* bit pattern. an example (in xdr format) is the double */
252  /* (63 237 133 81 81 108 3 32) which truncated to 23 bits precision should
253  */
254  /* become (63 237 133 81 64 0 0 0). however assigned to a float it becomes
255  */
256  /* (63 237 133 81 96 0 0 0). */
257  fTmp = *d;
258  dTmp = fTmp;
259 
260  Rast3d_float2xdrFloat(f, &xdrf);
261  Rast3d_float2Double(&xdrf, &xdrd2);
262  Rast3d_double2xdrDouble(&dTmp, &xdrd);
263 
264  c1 = (unsigned char *)&xdrd2;
265  c2 = (unsigned char *)&xdrd;
266 
267  /* printf ("%d %d (%d %d %d %d) (%d %d %d %d %d %d %d %d) (%d %d %d %d
268  * %d %d %d %d)\n", p1, p2, *((unsigned char *) &xdrf), *(((unsigned char *)
269  * &xdrf) + 1), *(((unsigned char *) &xdrf) + 2), *(((unsigned char *)
270  * &xdrf) + 3), *c1, *(c1 + 1), *(c1 + 2), *(c1 + 3), *(c1 + 4), *(c1 + 5),
271  * *(c1 + 6), *(c1 + 7), *c2, *(c2 + 1), *(c2 + 2), *(c2 + 3), *(c2 + 4),
272  * *(c2 + 5), *(c2 + 6), *(c2 + 7)); */
273 
274  if (((p1 != -1) && ((p1 < p2) || (p2 == -1))) ||
275  ((p1 == -1) && ((p2 > 23) || (p2 == -1))))
276  Rast3d_truncDouble(&xdrd, (p1 != -1 ? p1 : 23));
277  if ((p2 != -1) && (p2 < 23) && ((p2 < p1) || (p1 == -1)))
278  Rast3d_truncDouble(&xdrd2, p2);
279 
280  /* printf ("%d %d (%d %d %d %d) (%d %d %d %d %d %d %d %d) (%d %d %d %d %d
281  * %d %d %d)\n", p1, p2, *((unsigned char *) &xdrf), *(((unsigned char *)
282  * &xdrf) + 1), *(((unsigned char *) &xdrf) + 2), *(((unsigned char *)
283  * &xdrf) + 3), *c1, *(c1 + 1), *(c1 + 2), *(c1 + 3), *(c1 + 4), *(c1 + 5),
284  * *(c1 + 6), *(c1 + 7), *c2, *(c2 + 1), *(c2 + 2), *(c2 + 3), *(c2 + 4),
285  * *(c2 + 5), *(c2 + 6), *(c2 + 7)); */
286 
287  return (*c1 == *c2) && (*(c1 + 1) == *(c2 + 1)) &&
288  (*(c1 + 2) == *(c2 + 2)) && (*(c1 + 3) == *(c2 + 3)) &&
289  (*(c1 + 4) == *(c2 + 4)) && (*(c1 + 5) == *(c2 + 5)) &&
290  (*(c1 + 6) == *(c2 + 6)) && (*(c1 + 7) == *(c2 + 7));
291 }
292 
293 /*---------------------------------------------------------------------------*/
294 
295 static void compareFilesNocache(void *map, void *map2)
296 {
297  double n1 = 0, n2 = 0;
298  double *n1p, *n2p;
299  float *f1p, *f2p;
300  int x, y, z, correct;
301  int p1, p2;
302  int tileX, tileY, tileZ, typeIntern, typeIntern2;
303  int nx, ny, nz;
304 
305  p1 = Rast3d_tile_precision_map(map);
306  p2 = Rast3d_tile_precision_map(map2);
307 
308  Rast3d_get_tile_dimensions_map(map, &tileX, &tileY, &tileZ);
309  Rast3d_get_nof_tiles_map(map2, &nx, &ny, &nz);
310  typeIntern = Rast3d_tile_type_map(map);
311  typeIntern2 = Rast3d_tile_type_map(map2);
312 
313  n1p = &n1;
314  f1p = (float *)&n1;
315  n2p = &n2;
316  f2p = (float *)&n2;
317 
318  for (z = 0; z < nz * tileZ; z++) {
319  printf("comparing: z = %d\n", z);
320 
321  for (y = 0; y < ny * tileY; y++) {
322  for (x = 0; x < nx * tileX; x++) {
323 
324  Rast3d_get_block(map, x, y, z, 1, 1, 1, n1p, typeIntern);
325  Rast3d_get_block(map2, x, y, z, 1, 1, 1, n2p, typeIntern2);
326 
327  if (typeIntern == FCELL_TYPE) {
328  if (typeIntern2 == FCELL_TYPE)
329  correct = Rast3d_compareFloats(f1p, p1, f2p, p2);
330  else
331  correct = Rast3d_compareFloatDouble(f1p, p1, n2p, p2);
332  }
333  else {
334  if (typeIntern2 == FCELL_TYPE)
335  correct = Rast3d_compareFloatDouble(f2p, p2, n1p, p1);
336  else
337  correct = Rast3d_compareDoubles(n1p, p1, n2p, p2);
338  }
339 
340  if (!correct) {
341  int xTile, yTile, zTile, xOffs, yOffs, zOffs;
342 
343  Rast3d_coord2tile_coord(map2, x, y, z, &xTile, &yTile,
344  &zTile, &xOffs, &yOffs, &zOffs);
345  printf("(%d %d %d) (%d %d %d) (%d %d %d) %.20f %.20f\n", x,
346  y, z, xTile, yTile, zTile, xOffs, yOffs, zOffs, *n1p,
347  *n2p);
349  "compareFilesNocache: files don't match\n");
350  }
351  }
352  }
353  }
354 
355  printf("Files are identical up to precision.\n");
356 }
357 
358 /*---------------------------------------------------------------------------*/
359 
360 /*!
361  * \brief
362  *
363  * Compares the cell-values of file <em>f1</em> in mapset
364  * <em>mapset1</em> and file <em>f2</em> in mapset <em>mapset2</em>.
365  * The values are compared up to precision.
366  * Terminates in error if the files don't match.
367  * This function uses the more advanced features of the cache.
368  * The source code can be found in <em>filecompare.c</em>.
369  *
370  * \param f1
371  * \param mapset1
372  * \param f2
373  * \param mapset2
374  * \return void
375  */
376 
377 void Rast3d_compare_files(const char *f1, const char *mapset1, const char *f2,
378  const char *mapset2)
379 {
380  void *map, *map2;
381  double n1 = 0, n2 = 0;
382  double *n1p, *n2p;
383  float *f1p, *f2p;
384  int x, y, z, correct;
385  int p1, p2;
386  int rows, cols, depths;
387  int tileX, tileY, tileZ, typeIntern, typeIntern2, tileX2, tileY2, tileZ2;
388  int nx, ny, nz;
389 
390  printf("\nComparing %s and %s\n", f1, f2);
391 
395  if (map == NULL)
397  "Rast3d_compare_files: error in Rast3d_open_cell_old");
398 
399  Rast3d_print_header(map);
400 
401  map2 = Rast3d_open_cell_old(f2, mapset2, RASTER3D_DEFAULT_WINDOW,
404  if (map2 == NULL)
406  "Rast3d_compare_files: error in Rast3d_open_cell_old");
407 
408  Rast3d_print_header(map2);
409 
410  typeIntern = Rast3d_tile_type_map(map);
411  typeIntern2 = Rast3d_tile_type_map(map2);
412 
413  p1 = Rast3d_tile_precision_map(map);
414  p2 = Rast3d_tile_precision_map(map2);
415 
416  Rast3d_get_tile_dimensions_map(map, &tileX, &tileY, &tileZ);
417  Rast3d_get_tile_dimensions_map(map2, &tileX2, &tileY2, &tileZ2);
418  Rast3d_get_nof_tiles_map(map2, &nx, &ny, &nz);
419  Rast3d_get_coords_map(map, &rows, &cols, &depths);
420 
421  if ((!Rast3d_tile_use_cache_map(map)) ||
422  (!Rast3d_tile_use_cache_map(map2))) {
423  compareFilesNocache(map, map2);
424  Rast3d_close(map);
425  Rast3d_close(map2);
426  return;
427  }
428 
429  n1p = &n1;
430  f1p = (float *)&n1;
431  n2p = &n2;
432  f2p = (float *)&n2;
433 
434  Rast3d_autolock_on(map);
435  Rast3d_autolock_on(map2);
436  Rast3d_min_unlocked(map, cols / tileX + 1);
437 
438  Rast3d_get_coords_map(map2, &rows, &cols, &depths);
439  Rast3d_min_unlocked(map2, cols / tileX + 1);
440 
441  Rast3d_get_coords_map(map, &rows, &cols, &depths);
442  for (z = 0; z < depths; z++) {
443  printf("comparing: z = %d\n", z);
444 
445  if ((z % tileZ) == 0) {
446  if (!Rast3d_unlock_all(map))
448  "Rast3d_compare_files: error in Rast3d_unlock_all");
449  }
450  if ((z % tileZ2) == 0) {
451  if (!Rast3d_unlock_all(map2))
453  "Rast3d_compare_files: error in Rast3d_unlock_all");
454  }
455 
456  for (y = 0; y < rows; y++) {
457  for (x = 0; x < cols; x++) {
458  Rast3d_get_value_region(map, x, y, z, n1p, typeIntern);
459  Rast3d_get_value_region(map2, x, y, z, n2p, typeIntern2);
460 
461  Rast3d_is_null_value_num(n1p, typeIntern);
462  Rast3d_is_null_value_num(n2p, typeIntern2);
463 
464  if (typeIntern == FCELL_TYPE) {
465  if (typeIntern2 == FCELL_TYPE)
466  correct = Rast3d_compareFloats(f1p, p1, f2p, p2);
467  else
468  correct = Rast3d_compareFloatDouble(f1p, p1, n2p, p2);
469  }
470  else {
471  if (typeIntern2 == FCELL_TYPE)
472  correct = Rast3d_compareFloatDouble(f2p, p2, n1p, p1);
473  else
474  correct = Rast3d_compareDoubles(n1p, p1, n2p, p2);
475  }
476 
477  if (!correct) {
478  int xTile, yTile, zTile, xOffs, yOffs, zOffs;
479 
480  Rast3d_coord2tile_coord(map2, x, y, z, &xTile, &yTile,
481  &zTile, &xOffs, &yOffs, &zOffs);
483  "Rast3d_compare_files: files don't match\n");
484  }
485  }
486  }
487  }
488 
489  printf("Files are identical up to precision.\n");
490  Rast3d_close(map);
491  Rast3d_close(map2);
492 }
#define NULL
Definition: ccmath.h:32
void G_xdr_put_float(void *, const float *)
Definition: gis/xdr.c:84
void G_xdr_put_double(void *, const double *)
Definition: gis/xdr.c:94
int Rast3d_tile_type_map(RASTER3D_Map *)
Returns the type in which tiles of map are stored in memory.
Definition: headerinfo.c:156
void Rast3d_get_nof_tiles_map(RASTER3D_Map *, int *, int *, int *)
Returns the dimensions of the tile-cube used to tile the region of map. These numbers include partial...
Definition: headerinfo.c:50
void Rast3d_autolock_on(RASTER3D_Map *)
Turns autolock mode on.
Definition: tileread.c:340
void * Rast3d_open_cell_old(const char *, const char *, RASTER3D_Region *, int, int)
Opens existing g3d-file name in mapset. Tiles are stored in memory with type which must be any of FCE...
Definition: raster3d/open.c:80
void Rast3d_get_tile_dimensions_map(RASTER3D_Map *, int *, int *, int *)
Returns the tile dimensions used for map.
Definition: headerinfo.c:138
int Rast3d_unlock_all(RASTER3D_Map *)
Unlocks every tile in cache of map.
Definition: tileread.c:315
void Rast3d_min_unlocked(RASTER3D_Map *, int)
Sets the minimum number of unlocked tiles to minUnlocked. This function should be used in combination...
Definition: tileread.c:389
int Rast3d_is_null_value_num(const void *, int)
Definition: null.c:12
int Rast3d_tile_precision_map(RASTER3D_Map *)
Returns the precision used to store map.
Definition: headerinfo.c:292
void Rast3d_get_coords_map(RASTER3D_Map *, int *, int *, int *)
Returns the size of the region of map in cells.
Definition: headerinfo.c:18
void Rast3d_get_block(RASTER3D_Map *, int, int, int, int, int, int, void *, int)
Copies the cells contained in the block (cube) with vertices (x0, y0, z0) and (x0 + nx - 1,...
Definition: getblock.c:102
int Rast3d_tile_use_cache_map(RASTER3D_Map *)
Returns 1 if map uses cache, returns 0 otherwise.
Definition: headerinfo.c:308
void Rast3d_print_header(RASTER3D_Map *)
Prints the header information of map.
Definition: headerinfo.c:322
void Rast3d_get_value_region(RASTER3D_Map *, int, int, int, void *, int)
Returns in *value the cell-value of the cell with region-coordinate (x, y, z). The value returned is ...
Definition: getvalue.c:257
void Rast3d_coord2tile_coord(RASTER3D_Map *, int, int, int, int *, int *, int *, int *, int *, int *)
Converts cell-coordinates (x, y, z) into tile-coordinates (xTile, yTile, zTile) and the coordinate of...
Definition: tilemath.c:129
void Rast3d_fatal_error(const char *,...) __attribute__((format(printf
int Rast3d_close(RASTER3D_Map *)
Close 3D raster map files.
void Rast3d_compare_files(const char *f1, const char *mapset1, const char *f2, const char *mapset2)
Compares the cell-values of file f1 in mapset mapset1 and file f2 in mapset mapset2....
Definition: filecompare.c:377
dglInt32_t sign(dglInt32_t x)
Definition: flow.c:25
#define RASTER3D_DEFAULT_WINDOW
Definition: raster3d.h:29
#define RASTER3D_USE_CACHE_DEFAULT
Definition: raster3d.h:19
#define RASTER3D_TILE_SAME_AS_FILE
Definition: raster3d.h:11
#define FCELL_TYPE
Definition: raster.h:12
#define DCELL_TYPE
Definition: raster.h:13
#define x