1 #include "csf.h"
2 #include "csfimpl.h"
3 #include <math.h>
4
5 /* round up and down come up with a number such
6 * we have rounded to integer multiplication
7 * anyway this should hold:
8 POSTCOND(RoundUp( 5.3, 4) == 8 );
9 POSTCOND(RoundUp( 4 , 4) == 8 );
10 POSTCOND(RoundUp(-5.3, 4) == -4 );
11 POSTCOND(RoundUp(-4, 4) == 0 );
12 POSTCOND(RoundDown( 5.3, 4) == 4 );
13 POSTCOND(RoundDown( 4 , 4) == 0 );
14 POSTCOND(RoundDown(-5.3, 4) == -8 );
15 POSTCOND(RoundDown(-4 , 4) == -8 );
16 */
17
RoundDown(double v,double dfRound)18 static double RoundDown(
19 double v,
20 double dfRound)
21 {
22 double rVal = fmod(v, dfRound);
23 double x;
24 if(rVal == 0)
25 return v-dfRound;
26 if (v < 0)
27 x = v-dfRound-rVal;
28 else
29 x = v-rVal;
30 return x;
31 }
32
RoundUp(double v,double dfRound)33 static double RoundUp(
34 double v,
35 double dfRound)
36 {
37 double rVal = fmod(v, dfRound);
38 if(rVal == 0)
39 return v+dfRound;
40 if (v < 0)
41 return v-rVal;
42 else
43 return v+dfRound-rVal;
44 }
45
46 /* compute (xUL,yUL) and nrRows, nrCols from some coordinates
47 * RcomputeExtend computes parameters to create a raster maps
48 * from minimum and maximum x and y coordinates, projection information,
49 * cellsize and units. The resulting parameters are computed that the
50 * smallest raster map can be created that will include the two
51 * coordinates given, assuming a default angle of 0.
52 * Which coordinates are the maximum or minimum are
53 * determined by the function itself.
54 */
RcomputeExtend(REAL8 * xUL,REAL8 * yUL,size_t * nrRows,size_t * nrCols,double x_1,double y_1,double x_2,double y_2,CSF_PT projection,REAL8 cellSize,double rounding)55 void RcomputeExtend(
56 REAL8 *xUL, /* write-only, resulting xUL */
57 REAL8 *yUL, /* write-only, resulting yUL */
58 size_t *nrRows, /* write-only, resulting nrRows */
59 size_t *nrCols, /* write-only, resulting nrCols */
60 double x_1, /* first x-coordinate */
61 double y_1, /* first y-coordinate */
62 double x_2, /* second x-coordinate */
63 double y_2, /* second y-coordinate */
64 CSF_PT projection, /* required projection */
65 REAL8 cellSize, /* required cellsize, > 0 */
66 double rounding) /* assure that (xUL/rounding), (yUL/rounding)
67 * (xLL/rounding) and (yLL/rounding) will
68 * will all be an integers values > 0
69 */
70 {
71 /*
72 * xUL ______
73 | |
74 | |
75 | |
76 ------
77
78 */
79 double yLL,xUR = x_1 > x_2 ? x_1 : x_2;
80 *xUL = x_1 < x_2 ? x_1 : x_2;
81 *xUL = RoundDown(*xUL, rounding); /* Round down */
82 xUR = RoundUp( xUR, rounding); /* Round up */
83 POSTCOND(*xUL <= xUR);
84 *nrCols = (size_t)ceil((xUR - *xUL)/cellSize);
85 if (projection == PT_YINCT2B)
86 {
87 yLL = y_1 > y_2 ? y_1 : y_2; /* highest value at bottom */
88 *yUL = y_1 < y_2 ? y_1 : y_2; /* lowest value at top */
89 *yUL = RoundDown(*yUL, rounding);
90 yLL = RoundUp( yLL, rounding);
91 }
92 else
93 {
94 yLL = y_1 < y_2 ? y_1 : y_2; /* lowest value at bottom */
95 *yUL = y_1 > y_2 ? y_1 : y_2; /* highest value at top */
96 *yUL = RoundUp( *yUL, rounding);
97 yLL = RoundDown( yLL, rounding);
98 }
99 *nrRows = (size_t)ceil(fabs(yLL - *yUL)/cellSize);
100 }
101
102