1 /*
2 Gri - A language for scientific graphics programming
3 Copyright (C) 2008 Daniel Kelley
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; version 3 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License along
16 with this program; if not, write to the Free Software Foundation, Inc.,
17 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 */
19
20 #include <math.h>
21 #include <string.h>
22 #include "extern.hh"
23 #include "image_ex.hh"
24
25 bool mathCmd();
26 static bool column_math(double *Ptr, unsigned int n, int operator_position = 1);
27 static bool image_math(void);
28 static bool grid_data_math(void);
29 extern void highpass_image(void);
30 extern void lowpass_image(void);
31
32 bool
mathCmd()33 mathCmd()
34 {
35 if (_nword < 3) {
36 err("Proper format: `.x. = 10.0', `.x. *= 2', or `.x. = { rpn ... }'");
37 return false;
38 }
39 double number, old_value = 0.0;
40 // Variable
41 if (is_var(_word[0])) {
42 if (!getdnum(_word[2], &number))
43 return false;
44 if (word_is(1, "=")) {
45 // `.variable. = .value.'
46 PUT_VAR(_word[0], number);
47 } else if (word_is(1, "+=")) {
48 get_var(_word[0], &old_value);
49 old_value += number;
50 PUT_VAR(_word[0], old_value);
51 } else if (word_is(1, "-=")) {
52 get_var(_word[0], &old_value);
53 old_value -= number;
54 PUT_VAR(_word[0], old_value);
55 } else if (word_is(1, "*=")) {
56 get_var(_word[0], &old_value);
57 old_value *= number;
58 PUT_VAR(_word[0], old_value);
59 } else if (word_is(1, "/=")) {
60 get_var(_word[0], &old_value);
61 old_value /= number;
62 PUT_VAR(_word[0], old_value);
63 } else if (word_is(1, "^=")) {
64 get_var(_word[0], &old_value);
65 old_value = pow(old_value, number);
66 PUT_VAR(_word[0], old_value);
67 } else if (word_is(1, "_=")) {
68 get_var(_word[0], &old_value);
69 if (number < 0.0) {
70 err("Cannot do log to negative base");
71 return false;
72 }
73 if (!gr_missing(old_value) && old_value > 0.0) {
74 old_value = log(old_value) / log(number);
75 PUT_VAR(_word[0], old_value);
76 } else {
77 PUT_VAR(_word[0], gr_currentmissingvalue());
78 }
79 } else {
80 err("`\\", _word[1], "' is not a known operator for variables", "'\\");
81 return false;
82 }
83 } else if (_nword > 3 && word_is(0, "grid") && word_is(1, "data")) {
84 grid_data_math();
85 } else if (_nword > 3 && word_is(0, "grid") && word_is(1, "x")) {
86 if (_num_xmatrix_data)
87 column_math(_xmatrix, _num_xmatrix_data, 2);
88 else {
89 warning("Can't do math on xmatrix with no data in it already.");
90 return true;
91 }
92 } else if (_nword > 3 && word_is(0, "grid") && word_is(1, "y")) {
93 if (_num_ymatrix_data)
94 column_math(_ymatrix, _num_ymatrix_data, 2);
95 else {
96 warning("Can't do math on ymatrix with no data in it already.");
97 return true;
98 }
99 } else if (word_is(0, "u")) {
100 if (_colU.size())
101 column_math(_colU.begin(), _colU.size());
102 else {
103 warning("Can't do math on u-column with no data in it already.");
104 return true;
105 }
106 } else if (word_is(0, "v")) {
107 if (_colV.size())
108 column_math(_colV.begin(), _colV.size());
109 else {
110 warning("Can't do math on v-column with no data in it already.");
111 return true;
112 }
113 } else if (word_is(0, "x")) {
114 if (_colX.size())
115 column_math(_colX.begin(), _colX.size());
116 else {
117 warning("Can't do math on x-column with no data in it already.");
118 return true;
119 }
120 } else if (word_is(0, "y")) {
121 if (_colY.size())
122 column_math(_colY.begin(), _colY.size());
123 else {
124 warning("Can't do math on y-column with no data in it already.");
125 return true;
126 }
127 } else if (word_is(0, "z")) {
128 if (_colZ.size())
129 column_math(_colZ.begin(), _colZ.size());
130 else {
131 warning("Can't do math on z-column with no data in it already.");
132 return true;
133 }
134 } else if (word_is(0, "image")) {
135 image_math();
136 } else {
137 err("Cannot do mathematical operation on unrecognized item `\\",
138 _word[0],
139 "'.",
140 "\\");
141 return false;
142 }
143 return true;
144 }
145
146 static bool
column_math(double * Ptr,unsigned int n,int operator_position)147 column_math(double *Ptr, unsigned int n, int operator_position /*= 1*/)
148 {
149 Require(n > 0, gr_Error("Trying to do column-math on non-existent column."));
150 // Ensure enough space
151 if (word_is(0, "x")) {_colX.setDepth(n); Ptr = _colX.begin();}
152 else if (word_is(0, "y")) {_colY.setDepth(n); Ptr = _colY.begin();}
153 else if (word_is(0, "u")) {_colU.setDepth(n); Ptr = _colU.begin();}
154 else if (word_is(0, "v")) {_colV.setDepth(n); Ptr = _colV.begin();}
155 else if (word_is(0, "z")) {_colZ.setDepth(n); Ptr = _colZ.begin();}
156 else if (word_is(0, "WEIGHT")) {_colWEIGHT.setDepth(n); Ptr = _colWEIGHT.begin();}
157
158 // special case of `y = x...'
159 if (_nword == 3 && word_is(0, "y") && word_is(1, "=") && word_is(2, "x")) {
160 for (unsigned int i = 0; i < n; i++)
161 _colY[i] = _colX[i];
162 if (word_is(0, "y") && _need_y_axis && !_user_set_y_axis)
163 create_y_scale();
164 return true;
165 }
166
167 double number;
168 if (!getdnum(_word[operator_position + 1], &number))
169 return false;
170
171 // Special case of assignment
172 if (word_is(operator_position, "=")) {
173 // Make sure column can hold the data
174 for (unsigned int i = 0; i < n; i++, Ptr++) {
175 if (!gr_missing(*Ptr)) {
176 *Ptr = number;
177 }
178 }
179 if (word_is(0, "y") && _need_y_axis && !_user_set_y_axis)
180 create_y_scale();
181 return true;
182 }
183
184 // Do modification (OP=) cases
185 if (word_is(operator_position, "+=")) {
186 for (unsigned int i = 0; i < n; i++, Ptr++)
187 if (!gr_missing(*Ptr))
188 *Ptr += number;
189 } else if (word_is(operator_position, "-=")) {
190 for (unsigned int i = 0; i < n; i++, Ptr++)
191 if (!gr_missing(*Ptr))
192 *Ptr -= number;
193 } else if (word_is(operator_position, "*=")) {
194 for (unsigned int i = 0; i < n; i++, Ptr++)
195 if (!gr_missing(*Ptr))
196 *Ptr *= number;
197 } else if (word_is(operator_position, "/=")) {
198 for (unsigned int i = 0; i < n; i++, Ptr++)
199 if (!gr_missing(*Ptr))
200 *Ptr /= number;
201 } else if (word_is(operator_position, "^=")) {
202 for (unsigned int i = 0; i < n; i++, Ptr++)
203 if (!gr_missing(*Ptr))
204 *Ptr = pow(*Ptr, number);
205 } else if (word_is(operator_position, "_=")) {
206 if (number < 0.0) {
207 err("Cannot do log to negative base");
208 return false;
209 }
210 double lbase = log(number);
211 for (unsigned int i = 0; i < n; i++, Ptr++)
212 if (!gr_missing(*Ptr) && *Ptr > 0.0)
213 *Ptr = log(*Ptr) / lbase;
214 else
215 *Ptr = gr_currentmissingvalue();
216 } else {
217 err("`\\", _word[operator_position], "' not a known operator", "'\\");
218 return false;
219 }
220 if (word_is(0, "x") && _need_x_axis && !_user_set_x_axis)
221 create_x_scale();
222 if (word_is(0, "y") && _need_y_axis && !_user_set_y_axis)
223 create_y_scale();
224 return true;
225 }
226
227 static bool
grid_data_math()228 grid_data_math()
229 {
230 if (_nword < 4) {
231 err("No number for operation.");
232 return false;
233 }
234 if (!_grid_exists) {
235 err("First `read grid data'");
236 return false;
237 }
238 double number;
239 if (!getdnum(_word[3], &number)) {
240 err("Can't read number");
241 return false;
242 }
243 unsigned int row, col;
244 if (word_is(2, "=")) {
245 for (col = 0; col < _num_xmatrix_data; col++) {
246 for (row = 0; row < _num_ymatrix_data; row++) {
247 _f_xy(col, row) = number;
248 _legit_xy(col, row) = true;
249 }
250 }
251 _f_min = _f_max = number;
252 return true;
253 } else if (word_is(2, "+=")) {
254 for (col = 0; col < _num_xmatrix_data; col++)
255 for (row = 0; row < _num_ymatrix_data; row++)
256 _f_xy(col, row) += number;
257 _f_min += number;
258 _f_max += number;
259 return true;
260 } else if (word_is(2, "-=")) {
261 for (col = 0; col < _num_xmatrix_data; col++)
262 for (row = 0; row < _num_ymatrix_data; row++)
263 _f_xy(col, row) -= number;
264 _f_min -= number;
265 _f_max -= number;
266 return true;
267 } else if (word_is(2, "*=")) {
268 for (col = 0; col < _num_xmatrix_data; col++)
269 for (row = 0; row < _num_ymatrix_data; row++)
270 _f_xy(col, row) *= number;
271 _f_min *= number;
272 _f_max *= number;
273 return true;
274 } else if (word_is(2, "/=")) {
275 if (number == 0.0) {
276 warning("Can't divide by zero");
277 return true;
278 }
279 for (col = 0; col < _num_xmatrix_data; col++)
280 for (row = 0; row < _num_ymatrix_data; row++)
281 _f_xy(col, row) /= number;
282 _f_min /= number;
283 _f_max /= number;
284 return true;
285 } else if (word_is(2, "^=")) {
286 for (col = 0; col < _num_xmatrix_data; col++)
287 for (row = 0; row < _num_ymatrix_data; row++)
288 _f_xy(col, row) = pow(_f_xy(col, row), number);
289 matrix_limits(&_f_min, &_f_max);
290 return true;
291 } else if (word_is(2, "_=")) {
292 if (number < 0.0) {
293 err("Cannot do log to negative base");
294 return false;
295 }
296 double lbase = log(number);
297 for (col = 0; col < _num_xmatrix_data; col++)
298 for (row = 0; row < _num_ymatrix_data; row++)
299 if (!gr_missing(_f_xy(col, row)) && _f_xy(col, row) > 0.0)
300 _f_xy(col, row) = log(_f_xy(col, row)) / lbase;
301 else
302 _f_xy(col, row) = gr_currentmissingvalue();
303 matrix_limits(&_f_min, &_f_max);
304 return true;
305 } else {
306 err("`grid data' given unknown operator `\\", _word[2], ";", "\\");
307 return false;
308 }
309 }
310
311 // `image [grayscale|colorscale] OPERATOR OPERAND'
312 static bool
image_math()313 image_math()
314 {
315 unsigned int n;
316 unsigned char *Ptr;
317 double number, number_scaled=0.0;
318 switch (_nword) {
319 case 3:
320 // `image OPERATOR OPERAND'
321 if (word_is(1, "=")) {
322 if (word_is(2, "highpass"))
323 highpass_image();
324 else if (word_is(2, "lowpass"))
325 lowpass_image();
326 else {
327 err("`image =' what?");
328 return false;
329 }
330 return true;
331 }
332 if (!getdnum(_word[2], &number))
333 return false;
334 if (!_image.storage_exists) {
335 err("First `read image' or `convert grid to image'");
336 return false;
337 }
338 // Convert to image units
339 number_scaled *= 255.0 / (_image255 - _image0);
340 Ptr = _image.image;
341 n = _image.ras_length;
342 if (word_is(1, "+=")) {
343 for (unsigned int i = 0; i < n; i++, Ptr++)
344 *Ptr = int(*Ptr + number_scaled);
345 } else if (word_is(1, "-=")) {
346 for (unsigned int i = 0; i < n; i++, Ptr++)
347 *Ptr = int(*Ptr - number_scaled);
348 } else if (word_is(1, "*=")) {
349 for (unsigned int i = 0; i < n; i++, Ptr++)
350 *Ptr = int(*Ptr * number);
351 } else if (word_is(1, "/=")) {
352 for (unsigned int i = 0; i < n; i++, Ptr++)
353 *Ptr = int(*Ptr / number);
354 } else if (word_is(1, "^=")) {
355 for (unsigned int i = 0; i < n; i++, Ptr++)
356 *Ptr = int(pow(double(*Ptr), number));
357 } else if (word_is(1, "_=")) {
358 if (number < 0.0) {
359 err("Cannot do log to negative base");
360 return false;
361 }
362 double lbase = log(number);
363 for (unsigned int i = 0; i < n; i++, Ptr++)
364 if (!gr_missing(*Ptr) && *Ptr > 0)
365 *Ptr = int(log(double(*Ptr)) / lbase);
366 else
367 *Ptr = 0; // BUG -- should be something else
368 } else {
369 err("Cannot use operator `\\", _word[1], "' on images.", "'\\");
370 return false;
371 }
372 break;
373 case 4:
374 if (word_is(1, "colorscale") || word_is(1, "colourscale")) {
375 n = 3 * 256;
376 } else if (word_is(1, "greyscale") || word_is(1, "grayscale")) {
377 n = 256;
378 } else {
379 err("Second word must be `colorscale', `colourscale', `grayscale' or `greyscale'");
380 return false;
381 }
382 if (!getdnum(_word[3], &number))
383 return false;
384 Ptr = _imageTransform;
385 if (word_is(2, "+=")) {
386 for (unsigned int i = 0; i < n; i++, Ptr++) {
387 *Ptr = (int) (255.0 * (double(*Ptr) / 255.0 + number));
388 }
389 } else if (word_is(2, "-=")) {
390 for (unsigned int i = 0; i < n; i++, Ptr++) {
391 *Ptr = (int) (255.0 * (double(*Ptr) / 255.0 - number));
392 }
393 } else if (word_is(2, "*=")) {
394 for (unsigned int i = 0; i < n; i++, Ptr++) {
395 *Ptr = (int) (255.0 * (double(*Ptr) / 255.0 * number));
396 }
397 } else if (word_is(2, "/=")) {
398 for (unsigned int i = 0; i < n; i++, Ptr++) {
399 *Ptr = (int) (255.0 * (double(*Ptr) / 255.0 / number));
400 }
401 } else if (word_is(2, "^=")) {
402 for (unsigned int i = 0; i < n; i++, Ptr++) {
403 *Ptr = (int) (255.0 * pow(double(*Ptr) / 255.0, number));
404 }
405 } else if (word_is(2, "_=")) {
406 if (number < 0.0) {
407 err("Cannot do log to negative base");
408 return false;
409 }
410 double lbase = log(number);
411 for (unsigned int i = 0; i < n; i++, Ptr++)
412 *Ptr = int(255.0 * log(double(*Ptr) / 255.0) / lbase);
413 } else {
414 err("Cannot use operator `\\", _word[2], "' on image colorscale|grayscale", "'\\");
415 return false;
416 }
417 break;
418 default:
419 NUMBER_WORDS_ERROR;
420 return false;
421 }
422 return true;
423 }
424