GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-dcc4425b9d
c_reg.c
Go to the documentation of this file.
1 #include <math.h>
2 
3 #include <grass/gis.h>
4 #include <grass/raster.h>
5 
6 enum {
11  REGRESSION_P = 4
12 };
13 
14 static void regression(DCELL *result, DCELL *values, int n, int which)
15 {
16  DCELL xsum, ysum;
17  DCELL xbar, ybar;
18  DCELL numer, denom, denom2;
19  DCELL Rsq;
20  int count;
21  int i;
22 
23  xsum = ysum = 0.0;
24  count = 0;
25 
26  for (i = 0; i < n; i++) {
27  if (Rast_is_d_null_value(&values[i]))
28  continue;
29 
30  xsum += i;
31  ysum += values[i];
32  count++;
33  }
34 
35  if (count < 2) {
36  Rast_set_d_null_value(result, 1);
37  return;
38  }
39 
40  xbar = xsum / count;
41  ybar = ysum / count;
42 
43  numer = 0.0;
44  for (i = 0; i < n; i++)
45  if (!Rast_is_d_null_value(&values[i]))
46  numer += i * values[i];
47  numer -= count * xbar * ybar;
48 
49  denom = 0.0;
50  for (i = 0; i < n; i++)
51  if (!Rast_is_d_null_value(&values[i]))
52  denom += (DCELL)i * i;
53 
54  denom -= count * xbar * xbar;
55 
56  if (which >= REGRESSION_COEFF_DET || which == REGRESSION_T) {
57  denom2 = 0.0;
58  for (i = 0; i < n; i++)
59  if (!Rast_is_d_null_value(&values[i]))
60  denom2 += values[i] * values[i];
61  denom2 -= count * ybar * ybar;
62  Rsq = (numer * numer) / (denom * denom2);
63  }
64 
65  switch (which) {
66  case REGRESSION_SLOPE:
67  *result = numer / denom;
68  break;
69  case REGRESSION_OFFSET:
70  *result = ybar - xbar * numer / denom;
71  break;
73  *result = Rsq;
74  break;
75  case REGRESSION_T:
76  *result = sqrt(Rsq * (count - 2) / (1 - Rsq));
77  break;
78  default:
79  Rast_set_d_null_value(result, 1);
80  break;
81  }
82 
83  /* Check for NaN */
84  if (*result != *result)
85  Rast_set_d_null_value(result, 1);
86 }
87 
88 void c_reg_m(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
89 {
90  regression(result, values, n, REGRESSION_SLOPE);
91 }
92 
93 void c_reg_c(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
94 {
95  regression(result, values, n, REGRESSION_OFFSET);
96 }
97 
98 void c_reg_r2(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
99 {
100  regression(result, values, n, REGRESSION_COEFF_DET);
101 }
102 
103 void c_reg_t(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
104 {
105  regression(result, values, n, REGRESSION_T);
106 }
107 
108 static void regression_w(DCELL *result, DCELL (*values)[2], int n, int which)
109 {
110  DCELL xsum, ysum;
111  DCELL xbar, ybar;
112  DCELL numer, denom, denom2;
113  DCELL Rsq;
114  int count;
115  int i;
116 
117  xsum = ysum = 0.0;
118  count = 0;
119 
120  for (i = 0; i < n; i++) {
121  if (Rast_is_d_null_value(&values[i][0]))
122  continue;
123 
124  xsum += i * values[i][1];
125  ysum += values[i][0] * values[i][1];
126  count += values[i][1];
127  }
128 
129  if (count < 2) {
130  Rast_set_d_null_value(result, 1);
131  return;
132  }
133 
134  xbar = xsum / count;
135  ybar = ysum / count;
136 
137  numer = 0.0;
138  for (i = 0; i < n; i++)
139  if (!Rast_is_d_null_value(&values[i][0]))
140  numer += i * values[i][0] * values[i][1];
141  numer -= count * xbar * ybar;
142 
143  denom = 0.0;
144  for (i = 0; i < n; i++)
145  if (!Rast_is_d_null_value(&values[i][0]))
146  denom += (DCELL)i * i * values[i][1];
147 
148  denom -= count * xbar * xbar;
149 
150  if (which == REGRESSION_COEFF_DET || which == REGRESSION_T) {
151  denom2 = 0.0;
152  for (i = 0; i < n; i++)
153  if (!Rast_is_d_null_value(&values[i][0]))
154  denom2 += values[i][0] * values[i][0] * values[i][1];
155  denom2 -= count * ybar * ybar;
156  Rsq = (numer * numer) / (denom * denom2);
157  }
158 
159  switch (which) {
160  case REGRESSION_SLOPE:
161  *result = numer / denom;
162  break;
163  case REGRESSION_OFFSET:
164  *result = ybar - xbar * numer / denom;
165  break;
167  *result = Rsq;
168  break;
169  case REGRESSION_T:
170  *result = sqrt(Rsq * (count - 2) / (1 - Rsq));
171  break;
172  default:
173  Rast_set_d_null_value(result, 1);
174  break;
175  }
176 
177  /* Check for NaN */
178  if (*result != *result)
179  Rast_set_d_null_value(result, 1);
180 }
181 
182 void w_reg_m(DCELL *result, DCELL (*values)[2], int n,
183  const void *closure UNUSED)
184 {
185  regression_w(result, values, n, REGRESSION_SLOPE);
186 }
187 
188 void w_reg_c(DCELL *result, DCELL (*values)[2], int n,
189  const void *closure UNUSED)
190 {
191  regression_w(result, values, n, REGRESSION_OFFSET);
192 }
193 
194 void w_reg_r2(DCELL *result, DCELL (*values)[2], int n,
195  const void *closure UNUSED)
196 {
197  regression_w(result, values, n, REGRESSION_COEFF_DET);
198 }
199 
200 void w_reg_t(DCELL *result, DCELL (*values)[2], int n,
201  const void *closure UNUSED)
202 {
203  regression_w(result, values, n, REGRESSION_T);
204 }
void w_reg_m(DCELL *result, DCELL(*values)[2], int n, const void *closure UNUSED)
Definition: c_reg.c:182
void c_reg_c(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
Definition: c_reg.c:93
void w_reg_c(DCELL *result, DCELL(*values)[2], int n, const void *closure UNUSED)
Definition: c_reg.c:188
void c_reg_r2(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
Definition: c_reg.c:98
@ REGRESSION_P
Definition: c_reg.c:11
@ REGRESSION_COEFF_DET
Definition: c_reg.c:9
@ REGRESSION_OFFSET
Definition: c_reg.c:8
@ REGRESSION_T
Definition: c_reg.c:10
@ REGRESSION_SLOPE
Definition: c_reg.c:7
void w_reg_t(DCELL *result, DCELL(*values)[2], int n, const void *closure UNUSED)
Definition: c_reg.c:200
void c_reg_t(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
Definition: c_reg.c:103
void c_reg_m(DCELL *result, DCELL *values, int n, const void *closure UNUSED)
Definition: c_reg.c:88
void w_reg_r2(DCELL *result, DCELL(*values)[2], int n, const void *closure UNUSED)
Definition: c_reg.c:194
void Rast_set_d_null_value(DCELL *, int)
To set a number of DCELL raster values to NULL.
Definition: null_val.c:153
#define Rast_is_d_null_value(dcellVal)
Definition: defs/raster.h:405
double DCELL
Definition: gis.h:626
#define UNUSED
A macro for an attribute, if attached to a variable, indicating that the variable is not used.
Definition: gis.h:47
int count