1 /* Graphics_grey.cpp
2 *
3 * Copyright (C) 1992-2008,2011,2012,2015-2020 Paul Boersma
4 *
5 * This code 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; either version 2 of the License, or (at
8 * your option) any later version.
9 *
10 * This code is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 * See the GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this work. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 #include <stdlib.h>
20 #include "melder.h"
21 #include "Graphics.h"
22
23 //#define MAXGREYSIDE 12
24 #define MAXGREYSIDE 1000
25 #define MAXGREYPATH (2 * MAXGREYSIDE * (MAXGREYSIDE - 1) + 2)
26 #define MAXGREYEDGECONTOURS (2 * (MAXGREYSIDE - 1))
27 #define MAXGREYCLOSEDCONTOURS (MAXGREYPATH / 4 + 1)
28 #define MAXGREYEDGEPOINTS (4 * MAXGREYSIDE)
29
30 typedef struct {
31 integer numberOfPoints;
32 integer beginRow, beginCol, beginOri;
33 integer endRow, endCol, endOri;
34 integer lowerGrey, upperGrey;
35 autoVEC x, y;
36 } structEdgeContour, *EdgeContour;
37
EdgeContour_create(integer numberOfPoints)38 static EdgeContour EdgeContour_create (integer numberOfPoints) {
39 EdgeContour result = Melder_calloc (structEdgeContour, 1);
40 result -> numberOfPoints = numberOfPoints;
41 result -> x = zero_VEC (numberOfPoints);
42 result -> y = zero_VEC (numberOfPoints);
43 return result; // LEAK
44 }
EdgeContour_delete(EdgeContour e)45 static void EdgeContour_delete (EdgeContour e) {
46 e -> x. reset();
47 Melder_free (e);
48 }
49
50 typedef struct {
51 integer numberOfPoints;
52 int grey, drawn;
53 double xmin, xmax, ymin, ymax;
54 autoVEC x, y;
55 } structClosedContour, *ClosedContour;
56
ClosedContour_create(integer numberOfPoints)57 static ClosedContour ClosedContour_create (integer numberOfPoints) {
58 ClosedContour result = Melder_calloc (structClosedContour, 1);
59 result -> numberOfPoints = numberOfPoints;
60 result -> x = zero_VEC (numberOfPoints);
61 result -> y = zero_VEC (numberOfPoints);
62 return result; // LEAK
63 }
ClosedContour_delete(ClosedContour c)64 static void ClosedContour_delete (ClosedContour c) {
65 c -> x. reset();
66 Melder_free (c);
67 }
68
69 typedef struct {
70 integer ori, iContour, start, usedAsEntry, grey;
71 double val;
72 } structEdgePoint, *EdgePoint;
73
74 static Graphics theGraphics;
75 static integer numberOfEdgeContours;
76 static EdgeContour *edgeContours;
77 static integer numberOfEdgePoints;
78 static structEdgePoint *edgePoints;
79 static integer numberOfClosedContours;
80 static ClosedContour *closedContours;
81
82 static integer numberOfPoints;
83 static integer row1, row2, col1, col2;
84 static integer iBorder, numberOfBorders;
85 static autoBOOLMAT right, below;
86 static constMATVU data;
87 static double *border;
88 static autoVEC x, y;
89 static double dx, dy, xoff, yoff;
90
empty(integer row,integer col,integer ori)91 static int empty (integer row, integer col, integer ori)
92 {
93 if (ori == 3) { row ++; ori = 1; }
94 if (ori == 2) { col ++; ori = 4; }
95 if (ori == 1)
96 return (data [row] [col] < border [iBorder]) !=
97 (data [row] [col + 1] < border [iBorder]) &&
98 ! right [row - row1 + 1] [col - col1 + 1];
99 else /* ori == 4 */
100 return (data [row] [col] < border [iBorder]) !=
101 (data [row + 1] [col] < border [iBorder]) &&
102 ! below [row - row1 + 1] [col - col1 + 1];
103 }
note(integer row,integer col,integer ori)104 static integer note (integer row, integer col, integer ori)
105 {
106 ++ numberOfPoints;
107 Melder_assert (numberOfPoints <= MAXGREYPATH);
108 if (ori == 3) { row ++; ori = 1; }
109 if (ori == 2) { col ++; ori = 4; }
110 if (ori == 1)
111 {
112 right [row - row1 + 1] [col - col1 + 1] = true;
113 x [numberOfPoints] = xoff + (col + (border [iBorder] - data [row] [col]) /
114 (data [row] [col + 1] - data [row] [col])) * dx;
115 y [numberOfPoints] = yoff + row * dy;
116 }
117 else /* ori == 4 */
118 {
119 below [row - row1 + 1] [col - col1 + 1] = true;
120 x [numberOfPoints] = xoff + col * dx;
121 y [numberOfPoints] = yoff + (row + (border [iBorder] - data [row] [col]) /
122 (data [row + 1] [col] - data [row] [col])) * dy;
123 }
124 return 1;
125 }
126
fillGrey(integer numberOfPoints,constVECVU const & x,constVECVU const & y,int igrey)127 static void fillGrey (integer numberOfPoints, constVECVU const& x, constVECVU const& y, int igrey)
128 /* "igrey" is in between 1 and numberOfBorders + 1. */
129 {
130 Graphics_setGrey (theGraphics, 1.0 - (igrey - 1.0) / numberOfBorders);
131 Graphics_fillArea (theGraphics, numberOfPoints, & x [1], & y [1]);
132 }
133
makeEdgeContour(integer row0,integer col0,integer ori0)134 static void makeEdgeContour (integer row0, integer col0, integer ori0) {
135 numberOfPoints = 0;
136 integer row = row0, col = col0, ori = ori0;
137 note (row0, col0, ori0);
138
139 bool edge = false;
140 do {
141 bool clockwise = ! (ori & 1);
142 do { /* Preference for contours perpendicular to x == y. */
143 ori = (clockwise ? ori : ori + 2) % 4 + 1;
144 } while (! empty (row, col, ori));
145 switch (ori) {
146 case 1: edge = row == row1; break;
147 case 2: edge = col == col2 - 1; break;
148 case 3: edge = row == row2 - 1; break;
149 case 4: edge = col == col1; break;
150 }
151 if (! edge) {
152 switch (ori) {
153 case 1: row --; break;
154 case 2: col ++; break;
155 case 3: row ++; break;
156 case 4: col --; break;
157 }
158 ori = (ori + 1) % 4 + 1;
159 }
160 note (row, col, ori);
161 }
162 while (! edge);
163 Melder_assert (numberOfEdgeContours < MAXGREYEDGECONTOURS * numberOfBorders);
164
165 EdgeContour e = edgeContours [++ numberOfEdgeContours] = EdgeContour_create (numberOfPoints);
166 e -> beginRow = row0;
167 e -> beginCol = col0;
168 e -> beginOri = ori0;
169 e -> endRow = row;
170 e -> endCol = col;
171 e -> endOri = ori;
172
173 bool up = false;
174 switch (ori0) {
175 case 1: up = data [row0] [col0 + 1] > data [row0] [col0]; break;
176 case 2: up = data [row0 + 1] [col0 + 1] > data [row0] [col0 + 1]; break;
177 case 3: up = data [row0 + 1] [col0] > data [row0 + 1] [col0 + 1]; break;
178 case 4: up = data [row0] [col0] > data [row0 + 1] [col0]; break;
179 }
180 if (up) { e -> lowerGrey = iBorder; e -> upperGrey = iBorder + 1; }
181 else { e -> lowerGrey = iBorder + 1; e -> upperGrey = iBorder; }
182 for (integer iPoint = 1; iPoint <= numberOfPoints; iPoint ++) {
183 e -> x [iPoint] = x [iPoint];
184 e -> y [iPoint] = y [iPoint];
185 }
186 }
187
makeClosedContour(integer row0,integer col0,integer ori0)188 static void makeClosedContour (integer row0, integer col0, integer ori0) {
189 double x1, y1;
190 ClosedContour c;
191
192 numberOfPoints = 0;
193 integer row = row0, col = col0, ori = ori0;
194 do {
195 bool clockwise = ! (ori & 1);
196 do { /* Preference for contours perpendicular to x == y. */
197 ori = (clockwise ? ori : ori + 2) % 4 + 1;
198 } while (! empty (row, col, ori));
199 switch (ori) {
200 case 1: row --; break;
201 case 2: col ++; break;
202 case 3: row ++; break;
203 case 4: col --; break;
204 }
205 ori = (ori + 1) % 4 + 1;
206 note (row, col, ori);
207 }
208 while (row != row0 || col != col0 || ori != ori0);
209 Melder_assert (numberOfClosedContours < MAXGREYCLOSEDCONTOURS * numberOfBorders);
210 c = closedContours [++ numberOfClosedContours] = ClosedContour_create (numberOfPoints);
211
212 /* Find a point just inside or outside the contour. */
213 /* Find out whether it is above or below the contour. */
214
215 x1 = x [numberOfPoints];
216 y1 = y [numberOfPoints];
217 bool up = false;
218 if (ori == 3) { row ++; ori = 1; }
219 else if (ori == 2) { col ++; ori = 4; }
220 if (ori == 1) {
221 if (x1 > xoff + (col + 0.5) * dx)
222 x1 -= 0.01 * dx;
223 else
224 x1 += 0.01 * dx;
225 up = data [row] [col] + ((x1 - xoff) / dx - col) *
226 (data [row] [col + 1] - data [row] [col]) > border [iBorder];
227 } else { /* ori == 4 */
228 if (y1 > yoff + (row + 0.5) * dy)
229 y1 -= 0.01 * dy;
230 else
231 y1 += 0.01 * dy;
232 up = data [row] [col] + ((y1 - yoff) / dy - row) *
233 (data [row + 1] [col] - data [row] [col]) > border [iBorder];
234 }
235
236 /* Find out whether the point is inside or outside the contour. */
237
238 if (! NUMrotationsPointInPolygon (x1, y1, numberOfPoints,
239 x.asArgumentToFunctionThatExpectsOneBasedArray(), y.asArgumentToFunctionThatExpectsOneBasedArray()))
240 up = ! up;
241
242 double xmin = 1e308, xmax = -1e308, ymin = 1e308, ymax = -1e308;
243 c -> grey = up ? iBorder + 1 : iBorder;
244 for (int i = 1; i <= numberOfPoints; i ++) {
245 c -> x [i] = x [i];
246 c -> y [i] = y [i];
247 if (x [i] < xmin) xmin = x [i];
248 if (x [i] > xmax) xmax = x [i];
249 if (y [i] < ymin) ymin = y [i];
250 if (y [i] > ymax) ymax = y [i];
251 }
252 c -> xmin = xmin;
253 c -> xmax = xmax;
254 c -> ymin = ymin;
255 c -> ymax = ymax;
256 }
257
smallGrey()258 static void smallGrey () {
259 numberOfEdgeContours = 0;
260 numberOfClosedContours = 0;
261 for (iBorder = 1; iBorder <= numberOfBorders; iBorder ++) {
262 for (integer row = 1; row <= MAXGREYSIDE; row ++)
263 for (integer col = 1; col <= MAXGREYSIDE; col ++)
264 right [row] [col] = below [row] [col] = false;
265
266 /* Find all the edge contours of this border value. */
267
268 for (integer col = col1; col < col2; col ++)
269 if (empty (row1, col, 1))
270 makeEdgeContour (row1, col, 1);
271 for (integer row = row1; row < row2; row ++)
272 if (empty (row, col2 - 1, 2))
273 makeEdgeContour (row, col2 - 1, 2);
274 for (integer col = col2 - 1; col >= col1; col --)
275 if (empty (row2 - 1, col, 3))
276 makeEdgeContour (row2 - 1, col, 3);
277 for (integer row = row2 - 1; row >= row1; row --)
278 if (empty (row, col1, 4))
279 makeEdgeContour (row, col1, 4);
280
281 /* Find all the closed contours of this border value. */
282
283 for (integer row = row1 + 1; row < row2; row ++)
284 for (integer col = col1; col < col2; col ++)
285 if (empty (row, col, 1)) makeClosedContour (row, col, 1);
286 for (integer col = col1 + 1; col < col2; col ++)
287 for (integer row = row1; row < row2; row ++)
288 if (empty (row, col, 4)) makeClosedContour (row, col, 4);
289 }
290 numberOfEdgePoints = 2 * numberOfEdgeContours + 4;
291 Melder_assert (numberOfEdgePoints <= MAXGREYEDGEPOINTS * numberOfBorders);
292
293 /* Make a list of all points on the edge. */
294
295 /* The edge points include the four corner points. */
296
297 for (int i = 1; i <= 4; i ++) {
298 EdgePoint p = & edgePoints [i];
299 p -> ori = i;
300 p -> val = 0;
301 p -> iContour = 0;
302 p -> start = 0;
303 p -> usedAsEntry = 0;
304 p -> grey = -1;
305 }
306
307 /* The edge points include the first points of the edge contours. */
308
309 for (int i = 1; i <= numberOfEdgeContours; i ++) {
310 EdgeContour c = edgeContours [i];
311 EdgePoint p = & edgePoints [i + i + 3];
312 switch (p -> ori = c -> beginOri) {
313 case 1: p -> val = c -> x [1] - xoff - col1 * dx; break;
314 case 2: p -> val = c -> y [1] - yoff - row1 * dy; break;
315 case 3: p -> val = xoff + col2 * dx - c -> x [1]; break;
316 case 4: p -> val = yoff + row2 * dy - c -> y [1]; break;
317 }
318 p -> iContour = i;
319 p -> start = 1;
320 p -> usedAsEntry = 0;
321 p -> grey = c -> lowerGrey;
322 }
323
324 /* The edge points include the last points of the edge contours. */
325
326 for (int i = 1; i <= numberOfEdgeContours; i ++) {
327 EdgeContour c = edgeContours [i];
328 EdgePoint p = & edgePoints [i + i + 4];
329 switch (p -> ori = c -> endOri) {
330 case 1: p -> val = c -> x [c -> numberOfPoints] - xoff - col1 * dx; break;
331 case 2: p -> val = c -> y [c -> numberOfPoints] - yoff - row1 * dy; break;
332 case 3: p -> val = xoff + col2 * dx - c -> x [c -> numberOfPoints]; break;
333 case 4: p -> val = yoff + row2 * dy - c -> y [c -> numberOfPoints]; break;
334 }
335 p -> iContour = i;
336 p -> start = 0;
337 p -> usedAsEntry = 0;
338 p -> grey = c -> upperGrey;
339 }
340
341 /* Sort the list of edge points with keys Ori and Val. */
342 for (int i = 1; i < numberOfEdgePoints; i ++) {
343 structEdgePoint p;
344 int min = i, j;
345 for (j = i + 1; j <= numberOfEdgePoints; j ++)
346 if (edgePoints [min]. ori > edgePoints [j]. ori ||
347 (edgePoints [min]. ori == edgePoints [j]. ori && edgePoints [min]. val > edgePoints [j]. val))
348 min = j;
349 p = edgePoints [i];
350 edgePoints [i] = edgePoints [min];
351 edgePoints [min] = p;
352 }
353
354 {
355 for (int edge0 = 1; edge0 <= numberOfEdgePoints; edge0 ++)
356 if (edgePoints [edge0].grey > -1 && ! edgePoints [edge0].usedAsEntry) {
357 int iPoint = 0;
358 int edge1 = edge0;
359 int darkness;
360 do {
361 /*
362 Follow one edge contour.
363 */
364 EdgePoint p = & edgePoints [edge1];
365 integer iContour = p -> iContour;
366 EdgeContour c = edgeContours [iContour];
367 Melder_assert (iContour > 0);
368 darkness = p -> grey;
369 p -> usedAsEntry = 1;
370 if (p -> start) {
371 for (int i = 1; i <= c -> numberOfPoints; i ++) {
372 Melder_assert (iPoint < MAXGREYPATH);
373 x [++ iPoint] = c -> x [i];
374 y [iPoint] = c -> y [i];
375 }
376 for (int i = edge1 + 1; i <= numberOfEdgePoints; i ++)
377 if (edgePoints [i].iContour == iContour)
378 edge1 = i;
379 } else {
380 int edge1dummy = edge1;
381 for (integer i = c -> numberOfPoints; i >= 1; i --) {
382 Melder_assert (iPoint < MAXGREYPATH);
383 x [++ iPoint] = c -> x [i];
384 y [iPoint] = c -> y [i];
385 }
386 for (int i = 1; i <= edge1dummy - 1; i ++)
387 if (edgePoints [i].iContour == iContour)
388 edge1 = i;
389 }
390 edge1 = edge1 % numberOfEdgePoints + 1;
391
392 /*
393 Round some corners.
394 */
395 while (edgePoints [edge1].grey == -1) {
396 ++ iPoint;
397 Melder_assert (iPoint <= MAXGREYPATH);
398 switch (edgePoints [edge1].ori) {
399 case 1: x [iPoint] = xoff + col1 * dx; y [iPoint] = yoff + row1 * dy; break;
400 case 2: x [iPoint] = xoff + col2 * dx; y [iPoint] = yoff + row1 * dy; break;
401 case 3: x [iPoint] = xoff + col2 * dx; y [iPoint] = yoff + row2 * dy; break;
402 case 4: x [iPoint] = xoff + col1 * dx; y [iPoint] = yoff + row2 * dy; break;
403 }
404 edge1 = edge1 % numberOfEdgePoints + 1;
405 }
406 }
407 while (edge1 != edge0);
408 fillGrey (iPoint, x.all(), y.all(), darkness);
409 }
410 }
411 if (numberOfEdgeContours == 0) {
412 int i = 1;
413 while (i <= numberOfBorders && border [i] < data [row1] [col1])
414 i ++;
415 x [1] = x [4] = xoff + col1 * dx;
416 x [2] = x [3] = xoff + col2 * dx;
417 y [1] = y [2] = yoff + row1 * dy;
418 y [3] = y [4] = yoff + row2 * dy;
419 fillGrey (4, x.all(), y.all(), i);
420 }
421
422 /*
423 Iterate over all the closed contours.
424 Those that are not enclosed by any other contour, are filled first.
425 */
426 {
427 bool found;
428 do {
429 found = false;
430 for (integer i = 1; i <= numberOfClosedContours; i ++) {
431 ClosedContour ci = closedContours [i];
432 if (! ci -> drawn) {
433 bool enclosed = false;
434 integer j = 1;
435 while (j <= numberOfClosedContours && ! enclosed) {
436 ClosedContour cj = closedContours [j];
437 if (! cj -> drawn && j != i &&
438 ci -> xmin > cj -> xmin && ci -> xmax < cj -> xmax &&
439 ci -> ymin > cj -> ymin && ci -> ymax < cj -> ymax)
440 enclosed = NUMrotationsPointInPolygon (ci -> x [1], ci -> y [1],
441 cj -> numberOfPoints,
442 cj -> x.asArgumentToFunctionThatExpectsOneBasedArray(),
443 cj -> y.asArgumentToFunctionThatExpectsOneBasedArray());
444 j ++;
445 }
446 if (! enclosed) {
447 found = true;
448 fillGrey (ci -> numberOfPoints, ci -> x.all(), ci -> y.all(), ci -> grey);
449 ci -> drawn = 1;
450 }
451 }
452 }
453 } while (found);
454 }
455 Graphics_setGrey (theGraphics, 0.0);
456 for (int i = 1; i <= numberOfEdgeContours; i ++)
457 EdgeContour_delete (edgeContours [i]);
458 for (int i = 1; i <= numberOfClosedContours; i ++)
459 ClosedContour_delete (closedContours [i]);
460 }
461
Graphics_grey(Graphics me,constMATVU const & z,double x1WC,double x2WC,double y1WC,double y2WC,int _numberOfBorders,double borders[])462 void Graphics_grey (Graphics me, constMATVU const& z,
463 double x1WC, double x2WC, double y1WC, double y2WC,
464 int _numberOfBorders, double borders [])
465 {
466 if (z.nrow <= 1 || z.ncol <= 1)
467 return;
468
469 /* Static variables. */
470
471 theGraphics = me;
472 numberOfBorders = _numberOfBorders;
473 data = z;
474 border = borders;
475 dx = (x2WC - x1WC) / (z.ncol - 1);
476 dy = (y2WC - y1WC) / (z.nrow - 1);
477 xoff = x1WC - dx;
478 yoff = y1WC - dy;
479 if (NUMisEmpty (right.get())) {
480 right = zero_BOOLMAT (MAXGREYSIDE, MAXGREYSIDE); // BUG memory
481 below = zero_BOOLMAT (MAXGREYSIDE, MAXGREYSIDE);
482 x = zero_VEC (MAXGREYPATH);
483 y = zero_VEC (MAXGREYPATH);
484 edgeContours = Melder_calloc (EdgeContour, MAXGREYEDGECONTOURS * numberOfBorders) - 1;
485 closedContours = Melder_calloc (ClosedContour, MAXGREYCLOSEDCONTOURS * numberOfBorders) - 1;
486 edgePoints = Melder_calloc (structEdgePoint, MAXGREYEDGEPOINTS * numberOfBorders);
487 }
488
489 /* The matrix is subdivided into matrices with side MAXGREYSIDE, so that:
490 * 1. All the paths will fit into our memory (we have to remember them all).
491 * 2. The path for filling fits into the PostScript path, which may be max. 1500 points long.
492 */
493 for (row1 = 1; row1 < z.nrow; row1 += MAXGREYSIDE - 1) {
494 row2 = row1 + (MAXGREYSIDE - 1);
495 if (row2 > z.nrow) row2 = z.nrow;
496 for (col1 = 1; col1 < z.ncol; col1 += MAXGREYSIDE - 1) {
497 col2 = col1 + (MAXGREYSIDE - 1);
498 if (col2 > z.ncol) col2 = z.ncol;
499 smallGrey ();
500 }
501 }
502 }
503
504 /* End of file Graphics_grey.cpp */
505