GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
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 {
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  denom -= count * xbar * xbar;
54 
55  if (which >= REGRESSION_COEFF_DET || which == REGRESSION_T) {
56  denom2 = 0.0;
57  for (i = 0; i < n; i++)
58  if (!Rast_is_d_null_value(&values[i]))
59  denom2 += values[i] * values[i];
60  denom2 -= count * ybar * ybar;
61  Rsq = (numer * numer) / (denom * denom2);
62  }
63 
64  switch (which) {
65  case REGRESSION_SLOPE:
66  *result = numer / denom;
67  break;
68  case REGRESSION_OFFSET:
69  *result = ybar - xbar * numer / denom;
70  break;
72  *result = Rsq;
73  break;
74  case REGRESSION_T:
75  *result = sqrt(Rsq * (count - 2) / (1 - Rsq));
76  break;
77  default:
78  Rast_set_d_null_value(result, 1);
79  break;
80  }
81 
82  /* Check for NaN */
83  if (*result != *result)
84  Rast_set_d_null_value(result, 1);
85 }
86 
87 void c_reg_m(DCELL * result, DCELL * values, int n, const void *closure)
88 {
89  regression(result, values, n, REGRESSION_SLOPE);
90 }
91 
92 void c_reg_c(DCELL * result, DCELL * values, int n, const void *closure)
93 {
94  regression(result, values, n, REGRESSION_OFFSET);
95 }
96 
97 void c_reg_r2(DCELL * result, DCELL * values, int n, const void *closure)
98 {
99  regression(result, values, n, REGRESSION_COEFF_DET);
100 }
101 
102 void c_reg_t(DCELL * result, DCELL * values, int n, const void *closure)
103 {
104  regression(result, values, n, REGRESSION_T);
105 }
106 
107 static void regression_w(DCELL * result, DCELL(*values)[2], int n, int which)
108 {
109  DCELL xsum, ysum;
110  DCELL xbar, ybar;
111  DCELL numer, denom, denom2;
112  DCELL Rsq;
113  int count;
114  int i;
115 
116  xsum = ysum = 0.0;
117  count = 0;
118 
119  for (i = 0; i < n; i++) {
120  if (Rast_is_d_null_value(&values[i][0]))
121  continue;
122 
123  xsum += i * values[i][1];
124  ysum += values[i][0] * values[i][1];
125  count += values[i][1];
126  }
127 
128  if (count < 2) {
129  Rast_set_d_null_value(result, 1);
130  return;
131  }
132 
133  xbar = xsum / count;
134  ybar = ysum / count;
135 
136  numer = 0.0;
137  for (i = 0; i < n; i++)
138  if (!Rast_is_d_null_value(&values[i][0]))
139  numer += i * values[i][0] * values[i][1];
140  numer -= count * xbar * ybar;
141 
142  denom = 0.0;
143  for (i = 0; i < n; i++)
144  if (!Rast_is_d_null_value(&values[i][0]))
145  denom += (DCELL) i * i * values[i][1];
146 
147  denom -= count * xbar * xbar;
148 
149  if (which == REGRESSION_COEFF_DET || which == REGRESSION_T) {
150  denom2 = 0.0;
151  for (i = 0; i < n; i++)
152  if (!Rast_is_d_null_value(&values[i][0]))
153  denom2 += values[i][0] * values[i][0] * values[i][1];
154  denom2 -= count * ybar * ybar;
155  Rsq = (numer * numer) / (denom * denom2);
156  }
157 
158  switch (which) {
159  case REGRESSION_SLOPE:
160  *result = numer / denom;
161  break;
162  case REGRESSION_OFFSET:
163  *result = ybar - xbar * numer / denom;
164  break;
166  *result = Rsq;
167  break;
168  case REGRESSION_T:
169  *result = sqrt(Rsq * (count - 2) / (1 - Rsq));
170  break;
171  default:
172  Rast_set_d_null_value(result, 1);
173  break;
174  }
175 
176  /* Check for NaN */
177  if (*result != *result)
178  Rast_set_d_null_value(result, 1);
179 }
180 
181 void w_reg_m(DCELL * result, DCELL(*values)[2], int n, const void *closure)
182 {
183  regression_w(result, values, n, REGRESSION_SLOPE);
184 }
185 
186 void w_reg_c(DCELL * result, DCELL(*values)[2], int n, const void *closure)
187 {
188  regression_w(result, values, n, REGRESSION_OFFSET);
189 }
190 
191 void w_reg_r2(DCELL * result, DCELL(*values)[2], int n, const void *closure)
192 {
193  regression_w(result, values, n, REGRESSION_COEFF_DET);
194 }
195 
196 void w_reg_t(DCELL * result, DCELL(*values)[2], int n, const void *closure)
197 {
198  regression_w(result, values, n, REGRESSION_T);
199 }
void w_reg_c(DCELL *result, DCELL(*values)[2], int n, const void *closure)
Definition: c_reg.c:186
void c_reg_r2(DCELL *result, DCELL *values, int n, const void *closure)
Definition: c_reg.c:97
#define Rast_is_d_null_value(dcellVal)
Definition: defs/raster.h:420
double DCELL
Definition: gis.h:603
int count
void w_reg_m(DCELL *result, DCELL(*values)[2], int n, const void *closure)
Definition: c_reg.c:181
void c_reg_c(DCELL *result, DCELL *values, int n, const void *closure)
Definition: c_reg.c:92
void w_reg_r2(DCELL *result, DCELL(*values)[2], int n, const void *closure)
Definition: c_reg.c:191
void w_reg_t(DCELL *result, DCELL(*values)[2], int n, const void *closure)
Definition: c_reg.c:196
void c_reg_t(DCELL *result, DCELL *values, int n, const void *closure)
Definition: c_reg.c:102
void c_reg_m(DCELL *result, DCELL *values, int n, const void *closure)
Definition: c_reg.c:87
void Rast_set_d_null_value(DCELL *, int)
To set a number of DCELL raster values to NULL.
Definition: null_val.c:155