GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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] =
10  { 255, 128, 192, 224, 240, 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) ||
145  (*(c1 + 3) != 0))
146  e += 1023 - 127;
147  c = e / 16;
148 
149  *c2++ = (sign | c);
150 
151  c1++;
152 
153  c = e % 16;
154  *c2 = (c << 4);
155  *c2++ |= ((*c1 & (unsigned char)127) >> 3);
156 
157  *c2 = ((*c1++ & (unsigned char)7) << 5);
158  *c2++ |= (*c1 >> 3);
159 
160  *c2 = ((*c1++ & (unsigned char)7) << 5);
161  *c2++ |= (*c1 >> 3);
162 
163  *c2++ = ((*c1 & (unsigned char)7) << 5);
164 
165  *c2++ = (unsigned char)0;
166  *c2++ = (unsigned char)0;
167  *c2 = (unsigned char)0;
168 }
169 
170 /*---------------------------------------------------------------------------*/
171 
172 static int Rast3d_compareFloats(float *f1, int p1, float *f2, int p2)
173 {
174  unsigned char *c1, *c2;
175  float xdrf1, xdrf2;
176 
179 
180  Rast3d_float2xdrFloat(f1, &xdrf1);
181  Rast3d_float2xdrFloat(f2, &xdrf2);
182 
183  c1 = (unsigned char *)&xdrf1;
184  c2 = (unsigned char *)&xdrf2;
185 
186  /* printf ("%d %d (%d %d %d %d) (%d %d %d %d) %d\n", p1, p2, *c1, *(c1 + 1), *(c1 + 2), *(c1 + 3), *c2, *(c2 + 1), *(c2 + 2), *(c2 + 3), *f1 == *f2); */
187 
188  if ((p1 != -1) && (p1 < 23) && ((p1 < p2) || (p2 == -1)))
189  Rast3d_truncFloat(&xdrf2, p1);
190  if ((p2 != -1) && (p2 < 23) && ((p2 < p1) || (p1 == -1)))
191  Rast3d_truncFloat(&xdrf1, p2);
192 
193  /* printf ("%d %d (%d %d %d %d) (%d %d %d %d) %d\n", p1, p2, *c1, *(c1 + 1), *(c1 + 2), *(c1 + 3), *c2, *(c2 + 1), *(c2 + 2), *(c2 + 3), *f1 == *f2); */
194 
195  return (*c1 == *c2) && (*(c1 + 1) == *(c2 + 1)) &&
196  (*(c1 + 2) == *(c2 + 2))
197  && (*(c1 + 3) == *(c2 + 3));
198 }
199 
200 
201 /*---------------------------------------------------------------------------*/
202 
203 static int Rast3d_compareDoubles(double *d1, int p1, double *d2, int p2)
204 {
205  unsigned char *c1, *c2;
206  double xdrd1, xdrd2;
207 
210 
211  Rast3d_double2xdrDouble(d1, &xdrd1);
212  Rast3d_double2xdrDouble(d2, &xdrd2);
213 
214  c1 = (unsigned char *)&xdrd1;
215  c2 = (unsigned char *)&xdrd2;
216 
217  /* printf ("%d %d (%d %d %d %d %d %d %d %d) (%d %d %d %d %d %d %d %d)\n", p1, p2, *c1, *(c1 + 1), *(c1 + 2), *(c1 + 3), *(c1 + 4), *(c1 + 5), *(c1 + 6), *(c1 + 7), *c2, *(c2 + 1), *(c2 + 2), *(c2 + 3), *(c2 + 4), *(c2 + 5), *(c2 + 6), *(c2 + 7)); */
218 
219  if ((p1 != -1) && (p1 < 52) && ((p1 < p2) || (p2 == -1)))
220  Rast3d_truncDouble(&xdrd2, p1);
221  if ((p2 != -1) && (p2 < 52) && ((p2 < p1) || (p1 == -1)))
222  Rast3d_truncDouble(&xdrd1, p2);
223 
224  /* printf ("%d %d (%d %d %d %d %d %d %d %d) (%d %d %d %d %d %d %d %d)\n", p1, p2, *c1, *(c1 + 1), *(c1 + 2), *(c1 + 3), *(c1 + 4), *(c1 + 5), *(c1 + 6), *(c1 + 7), *c2, *(c2 + 1), *(c2 + 2), *(c2 + 3), *(c2 + 4), *(c2 + 5), *(c2 + 6), *(c2 + 7)); */
225 
226  return (*c1 == *c2) && (*(c1 + 1) == *(c2 + 1)) &&
227  (*(c1 + 2) == *(c2 + 2))
228  && (*(c1 + 3) == *(c2 + 3)) && (*(c1 + 4) == *(c2 + 4))
229  && (*(c1 + 5) == *(c2 + 5)) && (*(c1 + 6) == *(c2 + 6))
230  && (*(c1 + 7) == *(c2 + 7));
231 }
232 
233 
234 /*---------------------------------------------------------------------------*/
235 
236 static int Rast3d_compareFloatDouble(float *f, int p1, double *d, int p2)
237 {
238  unsigned char *c1, *c2;
239  float xdrf, fTmp;
240  double xdrd, xdrd2, dTmp;
241 
244 
245  /* need this since assigning a double to a float actually may change the */
246  /* bit pattern. an example (in xdr format) is the double */
247  /* (63 237 133 81 81 108 3 32) which truncated to 23 bits precision should */
248  /* become (63 237 133 81 64 0 0 0). however assigned to a float it becomes */
249  /* (63 237 133 81 96 0 0 0). */
250  fTmp = *d;
251  dTmp = fTmp;
252 
253  Rast3d_float2xdrFloat(f, &xdrf);
254  Rast3d_float2Double(&xdrf, &xdrd2);
255  Rast3d_double2xdrDouble(&dTmp, &xdrd);
256 
257  c1 = (unsigned char *)&xdrd2;
258  c2 = (unsigned char *)&xdrd;
259 
260  /* printf ("%d %d (%d %d %d %d) (%d %d %d %d %d %d %d %d) (%d %d %d %d %d %d %d %d)\n", p1, p2, *((unsigned char *) &xdrf), *(((unsigned char *) &xdrf) + 1), *(((unsigned char *) &xdrf) + 2), *(((unsigned char *) &xdrf) + 3), *c1, *(c1 + 1), *(c1 + 2), *(c1 + 3), *(c1 + 4), *(c1 + 5), *(c1 + 6), *(c1 + 7), *c2, *(c2 + 1), *(c2 + 2), *(c2 + 3), *(c2 + 4), *(c2 + 5), *(c2 + 6), *(c2 + 7)); */
261 
262 
263  if (((p1 != -1) && ((p1 < p2) || (p2 == -1))) ||
264  ((p1 == -1) && ((p2 > 23) || (p2 == -1))))
265  Rast3d_truncDouble(&xdrd, (p1 != -1 ? p1 : 23));
266  if ((p2 != -1) && (p2 < 23) && ((p2 < p1) || (p1 == -1)))
267  Rast3d_truncDouble(&xdrd2, p2);
268 
269  /* printf ("%d %d (%d %d %d %d) (%d %d %d %d %d %d %d %d) (%d %d %d %d %d %d %d %d)\n", p1, p2, *((unsigned char *) &xdrf), *(((unsigned char *) &xdrf) + 1), *(((unsigned char *) &xdrf) + 2), *(((unsigned char *) &xdrf) + 3), *c1, *(c1 + 1), *(c1 + 2), *(c1 + 3), *(c1 + 4), *(c1 + 5), *(c1 + 6), *(c1 + 7), *c2, *(c2 + 1), *(c2 + 2), *(c2 + 3), *(c2 + 4), *(c2 + 5), *(c2 + 6), *(c2 + 7)); */
270 
271  return (*c1 == *c2) && (*(c1 + 1) == *(c2 + 1)) &&
272  (*(c1 + 2) == *(c2 + 2))
273  && (*(c1 + 3) == *(c2 + 3)) && (*(c1 + 4) == *(c2 + 4))
274  && (*(c1 + 5) == *(c2 + 5)) && (*(c1 + 6) == *(c2 + 6))
275  && (*(c1 + 7) == *(c2 + 7));
276 }
277 
278 /*---------------------------------------------------------------------------*/
279 
280 static void compareFilesNocache(void *map, void *map2)
281 {
282  double n1 = 0, n2 = 0;
283  double *n1p, *n2p;
284  float *f1p, *f2p;
285  int x, y, z, correct;
286  int p1, p2;
287  int tileX, tileY, tileZ, typeIntern, typeIntern2;
288  int nx, ny, nz;
289 
290  p1 = Rast3d_tile_precision_map(map);
291  p2 = Rast3d_tile_precision_map(map2);
292 
293  Rast3d_get_tile_dimensions_map(map, &tileX, &tileY, &tileZ);
294  Rast3d_get_nof_tiles_map(map2, &nx, &ny, &nz);
295  typeIntern = Rast3d_tile_type_map(map);
296  typeIntern2 = Rast3d_tile_type_map(map2);
297 
298  n1p = &n1;
299  f1p = (float *)&n1;
300  n2p = &n2;
301  f2p = (float *)&n2;
302 
303  for (z = 0; z < nz * tileZ; z++) {
304  printf("comparing: z = %d\n", z);
305 
306  for (y = 0; y < ny * tileY; y++) {
307  for (x = 0; x < nx * tileX; x++) {
308 
309  Rast3d_get_block(map, x, y, z, 1, 1, 1, n1p, typeIntern);
310  Rast3d_get_block(map2, x, y, z, 1, 1, 1, n2p, typeIntern2);
311 
312  if (typeIntern == FCELL_TYPE) {
313  if (typeIntern2 == FCELL_TYPE)
314  correct = Rast3d_compareFloats(f1p, p1, f2p, p2);
315  else
316  correct = Rast3d_compareFloatDouble(f1p, p1, n2p, p2);
317  }
318  else {
319  if (typeIntern2 == FCELL_TYPE)
320  correct = Rast3d_compareFloatDouble(f2p, p2, n1p, p1);
321  else
322  correct = Rast3d_compareDoubles(n1p, p1, n2p, p2);
323  }
324 
325  if (!correct) {
326  int xTile, yTile, zTile, xOffs, yOffs, zOffs;
327 
328  Rast3d_coord2tile_coord(map2, x, y, z, &xTile, &yTile, &zTile,
329  &xOffs, &yOffs, &zOffs);
330  printf("(%d %d %d) (%d %d %d) (%d %d %d) %.20f %.20f\n",
331  x, y, z, xTile, yTile, zTile, xOffs, yOffs, zOffs,
332  *n1p, *n2p);
334  ("compareFilesNocache: files don't match\n");
335  }
336  }
337  }
338  }
339 
340  printf("Files are identical up to precision.\n");
341 }
342 
343 /*---------------------------------------------------------------------------*/
344 
345 
346 /*!
347  * \brief
348  *
349  * Compares the cell-values of file <em>f1</em> in mapset
350  * <em>mapset1</em> and file <em>f2</em> in mapset <em>mapset2</em>.
351  * The values are compared up to precision.
352  * Terminates in error if the files don't match.
353  * This function uses the more advanced features of the cache.
354  * The source code can be found in <em>filecompare.c</em>.
355  *
356  * \param f1
357  * \param mapset1
358  * \param f2
359  * \param mapset2
360  * \return void
361  */
362 
363 void
364 Rast3d_compare_files(const char *f1, const char *mapset1, const char *f2,
365  const char *mapset2)
366 {
367  void *map, *map2;
368  double n1 = 0, n2 = 0;
369  double *n1p, *n2p;
370  float *f1p, *f2p;
371  int x, y, z, correct;
372  int p1, p2;
373  int rows, cols, depths;
374  int tileX, tileY, tileZ, typeIntern, typeIntern2, tileX2, tileY2, tileZ2;
375  int nx, ny, nz;
376 
377  printf("\nComparing %s and %s\n", f1, f2);
378 
381  if (map == NULL)
382  Rast3d_fatal_error("Rast3d_compare_files: error in Rast3d_open_cell_old");
383 
384  Rast3d_print_header(map);
385 
386  map2 = Rast3d_open_cell_old(f2, mapset2, RASTER3D_DEFAULT_WINDOW,
388  if (map2 == NULL)
389  Rast3d_fatal_error("Rast3d_compare_files: error in Rast3d_open_cell_old");
390 
391  Rast3d_print_header(map2);
392 
393  typeIntern = Rast3d_tile_type_map(map);
394  typeIntern2 = Rast3d_tile_type_map(map2);
395 
396  p1 = Rast3d_tile_precision_map(map);
397  p2 = Rast3d_tile_precision_map(map2);
398 
399  Rast3d_get_tile_dimensions_map(map, &tileX, &tileY, &tileZ);
400  Rast3d_get_tile_dimensions_map(map2, &tileX2, &tileY2, &tileZ2);
401  Rast3d_get_nof_tiles_map(map2, &nx, &ny, &nz);
402  Rast3d_get_coords_map(map, &rows, &cols, &depths);
403 
404  if ((!Rast3d_tile_use_cache_map(map)) || (!Rast3d_tile_use_cache_map(map2))) {
405  compareFilesNocache(map, map2);
406  Rast3d_close(map);
407  Rast3d_close(map2);
408  return;
409  }
410 
411  n1p = &n1;
412  f1p = (float *)&n1;
413  n2p = &n2;
414  f2p = (float *)&n2;
415 
416  Rast3d_autolock_on(map);
417  Rast3d_autolock_on(map2);
418  Rast3d_min_unlocked(map, cols / tileX + 1);
419 
420  Rast3d_get_coords_map(map2, &rows, &cols, &depths);
421  Rast3d_min_unlocked(map2, cols / tileX + 1);
422 
423  Rast3d_get_coords_map(map, &rows, &cols, &depths);
424  for (z = 0; z < depths; z++) {
425  printf("comparing: z = %d\n", z);
426 
427  if ((z % tileZ) == 0) {
428  if (!Rast3d_unlock_all(map))
429  Rast3d_fatal_error("Rast3d_compare_files: error in Rast3d_unlock_all");
430  }
431  if ((z % tileZ2) == 0) {
432  if (!Rast3d_unlock_all(map2))
433  Rast3d_fatal_error("Rast3d_compare_files: error in Rast3d_unlock_all");
434  }
435 
436  for (y = 0; y < rows; y++) {
437  for (x = 0; x < cols; x++) {
438  Rast3d_get_value_region(map, x, y, z, n1p, typeIntern);
439  Rast3d_get_value_region(map2, x, y, z, n2p, typeIntern2);
440 
441  Rast3d_is_null_value_num(n1p, typeIntern);
442  Rast3d_is_null_value_num(n2p, typeIntern2);
443 
444  if (typeIntern == FCELL_TYPE) {
445  if (typeIntern2 == FCELL_TYPE)
446  correct = Rast3d_compareFloats(f1p, p1, f2p, p2);
447  else
448  correct = Rast3d_compareFloatDouble(f1p, p1, n2p, p2);
449  }
450  else {
451  if (typeIntern2 == FCELL_TYPE)
452  correct = Rast3d_compareFloatDouble(f2p, p2, n1p, p1);
453  else
454  correct = Rast3d_compareDoubles(n1p, p1, n2p, p2);
455  }
456 
457  if (!correct) {
458  int xTile, yTile, zTile, xOffs, yOffs, zOffs;
459 
460  Rast3d_coord2tile_coord(map2, x, y, z, &xTile, &yTile, &zTile,
461  &xOffs, &yOffs, &zOffs);
462  Rast3d_fatal_error("Rast3d_compare_files: files don't match\n");
463  }
464  }
465  }
466  }
467 
468  printf("Files are identical up to precision.\n");
469  Rast3d_close(map);
470  Rast3d_close(map2);
471 }
#define RASTER3D_USE_CACHE_DEFAULT
Definition: raster3d.h:20
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:77
void Rast3d_get_coords_map(RASTER3D_Map *, int *, int *, int *)
Returns the size of the region of map in cells.
Definition: headerinfo.c:19
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:386
void Rast3d_get_tile_dimensions_map(RASTER3D_Map *, int *, int *, int *)
Returns the tile dimensions used for map.
Definition: headerinfo.c:142
int Rast3d_tile_use_cache_map(RASTER3D_Map *)
Returns 1 if map uses cache, returns 0 otherwise.
Definition: headerinfo.c:315
#define RASTER3D_DEFAULT_WINDOW
Definition: raster3d.h:29
int Rast3d_close(RASTER3D_Map *)
Close 3D raster map files.
#define NULL
Definition: ccmath.h:32
void Rast3d_autolock_on(RASTER3D_Map *)
Turns autolock mode on.
Definition: tileread.c:337
#define x
void Rast3d_fatal_error(const char *,...) __attribute__((format(printf
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:105
int Rast3d_tile_precision_map(RASTER3D_Map *)
Returns the precision used to store map.
Definition: headerinfo.c:298
#define DCELL_TYPE
Definition: raster.h:13
void G_xdr_put_double(void *, const double *)
Definition: gis/xdr.c:88
void G_xdr_put_float(void *, const float *)
Definition: gis/xdr.c:78
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:256
#define RASTER3D_TILE_SAME_AS_FILE
Definition: raster3d.h:12
void Rast3d_print_header(RASTER3D_Map *)
Prints the header information of map.
Definition: headerinfo.c:331
#define FCELL_TYPE
Definition: raster.h:12
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:136
int Rast3d_is_null_value_num(const void *, int)
Definition: null.c:12
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. The values are compared up to precision. Terminates in error if the files don&#39;t match. This function uses the more advanced features of the cache. The source code can be found in filecompare.c.
Definition: filecompare.c:364
dglInt32_t sign(dglInt32_t x)
Definition: flow.c:25
int Rast3d_unlock_all(RASTER3D_Map *)
Unlocks every tile in cache of map.
Definition: tileread.c:312
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:51
int Rast3d_tile_type_map(RASTER3D_Map *)
Returns the type in which tiles of map are stored in memory.
Definition: headerinfo.c:161