-
Notifications
You must be signed in to change notification settings - Fork 0
/
inc.conrec.php
419 lines (384 loc) · 20.8 KB
/
inc.conrec.php
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
<?php
// ------------------------------------------------------------------
// ------------------------------------------------------------------
// CONREC is a straightforward method of contouring some surface represented
// as a regular triangular mesh.
//
// Ported from C / Fortran code by Paul Borke.
// http://paulbourke.net/papers/conrec/
//
// PHP implementation by Ashley Norris.
// http://norris.org.au/
//
// Contouring aids in visualizing three dimensional surfaces on a two dimensional
// medium (on paper or in this case a computer graphics screen). Two most common
// applications are displaying topological features of an area on a map or the air
// pressure on a weather map. In all cases some parameter is plotted as a function
// of two variables, the longitude and latitude or x and y axis. One problem with
// computer contouring is the process is usually CPU intensive and the algorithms
// often use advanced mathematical techniques making them susceptible to error.
//
// ------------------------------------------------------------------
// ------------------------------------------------------------------
//
// This PHP implementation accepts a 2D scalar field (a matrix) and returns
// a set of line segments that would typically then be rendered in some way by the user.
//
// The form of the returned data has been made as simple as possible, and could
// be easily edited (to minimise memory use) by anyone if they so wish.
//
// The returned data is an array of this structure:
// [top] (array of contours) // An array containing all of the contours created by this algorithm
// -> [contour] (mixed array)
// -> 'value' (double) // The scalar value (height) of the contour
// -> ['segments'] (array of segments) // An array containing all the segments (not necessarily continuous) that make up the contour
// -> 'x1' (double)
// -> 'y1' (double)
// -> 'x2' (double)
// -> 'y2' (double)
//
// Best thing to do is to use the function and then do a "print_r();" on its returned value to explore the structure returned.
//
// ------------------------------------------------------------------
// ------------------------------------------------------------------
//
// This library file produces no output, and contains a definition for one useful function: contour();
//
// There are some other functions defined at the end of the file for internal use, which would be useless for external calls.
//
// ------------------------------------------------------------------
// ------------------------------------------------------------------
// Function contour($d, $x, $y, $z, $SKIP_QC = false);
//
// === Input Parameters ===
//
// $d => 2D matrix of values, typically of type double
// $x => 1D array of the x-coordinate for each column
// $y => 1D array of the y-coordinate for each row
// $z => 1D array of the Z value to use for each contour **OR** an integer number of contours to automatically create
// $SKIP_QC => boolean, set to true to skip the Quality Control stage (speed benefit for large datasets)
//
// Some notes:
// - $d is a 2D jagged array, with X on the first dim and Y on the second.
// I say jagged as PHP does not support a true native multidimensional array
// type in the same way as other languages do. Thus, (C, Fortran, etc...) d[i,j] === (PHP) $d[i][j]
// - $x must be the same length as the width (first dim) of $d
// - $y must be the same length as the height (second dim) of $d
// - The $x and $y arrays are not specifically "needed" to create contours, but are
// used to return contour line segments in that coordinate system, rather than a unit system of col and row number.
// - $z takes advantage of PHP's flexible typing (yes, it really is good for some things!)
// Either, pass in an array of values and the number of contours returned will equal the length of $z
// Or, pass in an integer, and that many contours will be returned.
// ------------------------------------------------------------------
// ------------------------------------------------------------------
function CONREC_contour($d, $x, $y, $z, $SKIP_QC = false)
{
// ==================================================
// QC the inputs (turn off for large datasets)
// If wanting to improve performance, but maintain QC, then feel free to weave the QC checks into the actual algorithm code below
// ==================================================
if (!$SKIP_QC) {
// Basic checks first
if (!is_array($d)) { die('Function contour(); failed QC. Parameter $d is not an array.'); }
if (!is_array($x)) { die('Function contour(); failed QC. Parameter $x is not an array.'); }
if (!is_array($y)) { die('Function contour(); failed QC. Parameter $y is not an array.'); }
// NOTE: the count() function in PHP is NOT recursive by default
if (count($d) <= 1) { die('Function contour(); failed QC. Dimension 1 of $d has length smaller than or equal to 1.'); }
if (count($x) <= 1) { die('Function contour(); failed QC. Array $x has length smaller than or equal to 1.'); }
if (count($y) <= 1) { die('Function contour(); failed QC. Array $y has length smaller than or equal to 1.'); }
// Check first dim of $d against $x
if (count($d) != count($x)) { die('Function contour(); failed QC. Dimension 1 of $d ['.count($d).'] does not match length of $x ['.count($x).'].'); }
// Loop first dim of $d and check second dim for all rows (as a side-effect this ensures $d is regular)
foreach($d as $index=>$row) {
if (!is_int($index)) { die('Function contour(); failed QC. Row "'.$index.'" of $d must have an integer index.'); }
if (count($row) != count($y)) { die('Function contour(); failed QC. Dimension 2 of row "'.$index.'" of $d does not match length of $y.'); }
// Check indexes of row
foreach(array_keys($row) as $col_index) {
if (!is_int($col_index)) { die('Function contour(); failed QC. On row "'.$index.'", col "'.$col_index.'" of $d must have an integer index.'); }
}
}
// $z is either an array, or an integer
if (!is_array($z) AND !is_int($z)) { die('Function contour(); failed QC. Parameter $z must be either an array or an integer.'); }
}
// ==================================================
// ==================================================
// While PHP does not require us to declare our variables, we do so here anyway, as per the original C code
$x1 = 0.0;
$x2 = 0.0;
$y1 = 0.0;
$y2 = 0.0;
// ==================================================
// For these arrays we can not specify a fixed length, we just init them with zeros and start using them
// ==================================================
$h = array(0,0,0,0,0); // Length 5
$sh = array(0,0,0,0,0); // Length 5
$xh = array(0,0,0,0,0); // Length 5
$yh = array(0,0,0,0,0); // Length 5
// ==================================================
// When getting the range of the X and Y values we can use the first row of the matrix thanks to the QC step above
// ==================================================
// NOTE: the count() function in PHP is NOT recursive by default
$ilb = 0; // Can assume this thanks to QC step
$iub = count($d) - 1; // Counts first dim
$jlb = 0; // Can assume this thanks to QC step
$jub = count($d[0]) - 1; // Counts second dim of first row, which we can do thanks to QC step
if (is_array($z)) {
$nc = count($z);
}
else if (is_int($z)) {
$nc = $z;
}
else {
// Sensible hard-coded default value, in case QC was skipped
$nc = 10;
}
// ==================================================
// Automatic contours (in case only the number has been defined)
// ==================================================
if (!is_array($z)) {
$z = CONREC_auto_z($nc,$d);
}
// ==================================================
// Sort $z
sort($z);
// ==================================================
// The indexing of im and jm should be noted as it has to start from zero
// unlike the fortran counter part
$im = array(0, 1, 1, 0 );
$jm = array( 0, 0, 1, 1 );
// ==================================================
// Note that castab is arranged differently from the FORTRAN code because
// Fortran and C/C++ arrays are transposed of each other, in this case
// it is more tricky as castab is in 3 dimensions
$castab = array(
array( array( 0, 0, 8 ), array( 0, 2, 5 ), array( 7, 6, 9 ) ),
array( array( 0, 3, 4 ), array( 1, 3, 1 ), array( 4, 3, 0 ) ),
array( array( 9, 6, 7 ), array( 5, 2, 0 ), array( 8, 0, 0 ) )
);
// ==================================================
// We are now ready to perform the contouring, init the return array
$return_array = array();
foreach($z as $index=>$value) {
$return_array[$index] = array('value'=>$value,'segments'=>array());
}
// ==================================================
// Primary loop: i (row, or $d first dimension) loop begins here
for ($i = $ilb; $i <= $iub - 1; $i++)
{
// ==================================================
// Primary loop: j (col, or $d second dimension) loop begins here
for ($j = $jub - 1; $j >= $jlb; $j--)
{
// ==================================================
// Find the max and min value of the corners
$corners = array($d[$i][$j], $d[$i][$j + 1], $d[$i + 1][$j], $d[$i + 1][$j + 1]);
$dmin = min($corners);
$dmax = max($corners);
// Perform elimination of trivial cases
if ($dmax >= $z[0] && $dmin <= $z[$nc - 1])
{
// Start k (contour level) loop
for ($k = 0; $k < $nc; $k++)
{
// Test against this contour level
if ($z[$k] >= $dmin && $z[$k] <= $dmax)
{
// Loop over $m (corners)
for ($m = 4; $m >= 0; $m--)
{
if ($m > 0)
{
// The indexing of im and jm should be noted as it has to
// start from zero
$h[$m] = $d[$i + $im[$m - 1]][ $j + $jm[$m - 1]] - $z[$k];
$xh[$m] = $x[$i + $im[$m - 1]];
$yh[$m] = $y[$j + $jm[$m - 1]];
}
else
{
// Compute average
$h[0] = 0.25 * ($h[1] + $h[2] + $h[3] + $h[4]);
$xh[0] = 0.5 * ($x[$i] + $x[$i + 1]);
$yh[0] = 0.5 * ($y[$j] + $y[$j + 1]);
}
if ($h[$m] > 0.0)
{
$sh[$m] = 1;
}
else if ($h[$m] < 0.0)
{
$sh[$m] = -1;
}
else
{
$sh[$m] = 0;
}
} // end $m loop
// Note: at this stage the relative heights of the corners and the
// centre are in the h array, and the corresponding coordinates are
// in the xh and yh arrays. The centre of the box is indexed by 0
// and the 4 corners by 1 to 4 as shown below.
// Each triangle is then indexed by the parameter m, and the 3
// vertices of each triangle are indexed by parameters m1,m2,and
// m3.
// It is assumed that the centre of the box is always vertex 2
// though this is important only when all 3 vertices lie exactly on
// the same contour level, in which case only the side of the box
// is drawn.
//
// vertex 4
// +-------------------+ vertex 3
// | \ / |
// | \ m-3 / |
// | \ / |
// | \ / |
// | m=2 X m=2 | the centre is vertex 0
// | / \ |
// | / \ |
// | / m=1 \ |
// | / \ |
// +-------------------+ vertex 2
// vertex 1
//
// Scan each triangle in the box
for ($m = 1; $m <= 4; $m++)
{
$m1 = $m;
$m2 = 0;
$m3 = 0;
if ($m != 4)
{
$m3 = $m + 1;
}
else
{
$m3 = 1;
}
$caseValue = $castab[$sh[$m1] + 1][$sh[$m2] + 1][$sh[$m3] + 1];
if ($caseValue != 0)
{
switch ($caseValue)
{
case 1: // Line between vertices 1 and 2
$x1 = $xh[$m1];
$y1 = $yh[$m1];
$x2 = $xh[$m2];
$y2 = $yh[$m2];
break;
case 2: // Line between vertices 2 and 3
$x1 = $xh[$m2];
$y1 = $yh[$m2];
$x2 = $xh[$m3];
$y2 = $yh[$m3];
break;
case 3: // Line between vertices 3 and 1
$x1 = $xh[$m3];
$y1 = $yh[$m3];
$x2 = $xh[$m1];
$y2 = $yh[$m1];
break;
case 4: // Line between vertex 1 and side 2-3
$x1 = $xh[$m1];
$y1 = $yh[$m1];
$x2 = CONREC_sect($m2, $m3, $h, $xh);
$y2 = CONREC_sect($m2, $m3);
break;
case 5: // Line between vertex 2 and side 3-1
$x1 = $xh[$m2];
$y1 = $yh[$m2];
$x2 = CONREC_sect($m3, $m1, $h, $xh);
$y2 = CONREC_sect($m3, $m1, $h, $yh);
break;
case 6: // Line between vertex 3 and side 1-2
$x1 = $xh[$m3];
$y1 = $yh[$m3];
$x2 = CONREC_sect($m1, $m2, $h, $xh);
$y2 = CONREC_sect($m1, $m2, $h, $yh);
break;
case 7: // Line between sides 1-2 and 2-3
$x1 = CONREC_sect($m1, $m2, $h, $xh);
$y1 = CONREC_sect($m1, $m2, $h, $yh);
$x2 = CONREC_sect($m2, $m3, $h, $xh);
$y2 = CONREC_sect($m2, $m3, $h, $yh);
break;
case 8: // Line between sides 2-3 and 3-1
$x1 = CONREC_sect($m2, $m3, $h, $xh);
$y1 = CONREC_sect($m2, $m3, $h, $yh);
$x2 = CONREC_sect($m3, $m1, $h, $xh);
$y2 = CONREC_sect($m3, $m1, $h, $yh);
break;
case 9: // Line between sides 3-1 and 1-2
$x1 = CONREC_sect($m3, $m1, $h, $xh);
$y1 = CONREC_sect($m3, $m1, $h, $yh);
$x2 = CONREC_sect($m1, $m2, $h, $xh);
$y2 = CONREC_sect($m1, $m2, $h, $yh);
break;
default:
break;
}
// ==================================================
// OUTPUT to return array structure
// ==================================================
$return_array[$k]['segments'][] = array('x1'=>$x1,'y1'=>$y1,'x2'=>$x2,'y2'=>$y2);
// NOTE: If wishing to improve performance, or reduce memory use, just put your
// rendering routine right here.
// ==================================================
} // end if: caseValue != 0
} // end for: m loop
}
} // end k loop
} // end if: trivial case elimination (outside contouring range)
} // end for: primary j loop
} // end for: primary i loop
// Return the computed contours
return $return_array;
} // end function: conrec();
// ----------------------------------------------------------------
// Internal functions
// ----------------------------------------------------------------
// As we don't have access to the $xh and $yh parameters we pass them into this single function to perform the segment computation
function CONREC_sect($p1_index, $p2_index, $h, $dimh) {
return ($h[$p2_index] * $dimh[$p1_index] - $h[$p1_index] * $dimh[$p2_index]) / ($h[$p2_index] - $h[$p1_index]);
}
function CONREC_auto_z($nc, $array) {
// ==================================================
// NOTE: Simply dividing the range by the number of desired contours doesn't
// actually work very well, as the max and min values do not produce
// any contours at all. To correct for this we add 2 to the desired number of
// contours, divide the range by that number and then use that step to return the
// actual contour values to use.
// ==================================================
$z = array();
$d_min = CONREC_arraymin($array);
$d_max = CONREC_arraymax($array);
$d_range = $d_max - $d_min;
$contour_step = $d_range / ($nc + 2 - 1); // See above for an explanation of this algebra :)
for($tmp_i = 1;$tmp_i <= $nc;$tmp_i++) {
$z[] = $d_min + ($contour_step * $tmp_i);
}
// $z should now contain $nc elements
return $z;
}
function CONREC_arraymax($array){
if (is_array($array)){
foreach($array as $key => $value) {
$array[$key] = CONREC_arraymax($value);
}
return max($array);
} else{
return $array;
}
}
function CONREC_arraymin($array){
if (is_array($array)) {
foreach($array as $key => $value) {
$array[$key] = CONREC_arraymin($value);
}
return min($array);
} else{
return $array;
}
}
// ----------------------------------------------------------------
// End of file.
// ----------------------------------------------------------------