source: trunk/xgraph/jpgraph/jpgraph_contour.php@ 45

Last change on this file since 45 was 42, checked in by marrucho, 11 years ago
File size: 21.2 KB
Line 
1<?php
2/*=======================================================================
3// File: JPGRAPH_CONTOUR.PHP
4// Description: Contour plot
5// Created: 2009-03-08
6// Ver: $Id: jpgraph_contour.php 1870 2009-09-29 04:24:18Z ljp $
7//
8// Copyright (c) Asial Corporation. All rights reserved.
9//========================================================================
10*/
11require_once('jpgraph_meshinterpolate.inc.php');
12define('HORIZ_EDGE',0);
13define('VERT_EDGE',1);
14
15/**
16 * This class encapsulates the core contour plot algorithm. It will find the path
17 * of the specified isobars in the data matrix specified. It is assumed that the
18 * data matrix models an equspaced X-Y mesh of datavalues corresponding to the Z
19 * values.
20 *
21 */
22class Contour {
23
24 private $dataPoints = array();
25 private $nbrCols=0,$nbrRows=0;
26 private $horizEdges = array(), $vertEdges=array();
27 private $isobarValues = array();
28 private $stack = null;
29 private $isobarCoord = array();
30 private $nbrIsobars = 10, $isobarColors = array();
31 private $invert = true;
32 private $highcontrast = false, $highcontrastbw = false;
33
34 /**
35 * Create a new contour level "algorithm machine".
36 * @param $aMatrix The values to find the contour from
37 * @param $aIsobars Mixed. If integer it determines the number of isobars to be used. The levels are determined
38 * automatically as equdistance between the min and max value of the matrice.
39 * If $aIsobars is an array then this is interpretated as an array of values to be used as isobars in the
40 * contour plot.
41 * @return an instance of the contour algorithm
42 */
43 function __construct($aMatrix,$aIsobars=10, $aColors=null) {
44
45 $this->nbrRows = count($aMatrix);
46 $this->nbrCols = count($aMatrix[0]);
47 $this->dataPoints = $aMatrix;
48
49 if( is_array($aIsobars) ) {
50 // use the isobar values supplied
51 $this->nbrIsobars = count($aIsobars);
52 $this->isobarValues = $aIsobars;
53 }
54 else {
55 // Determine the isobar values automatically
56 $this->nbrIsobars = $aIsobars;
57 list($min,$max) = $this->getMinMaxVal();
58 $stepSize = ($max-$min) / $aIsobars ;
59 $isobar = $min+$stepSize/2;
60 for ($i = 0; $i < $aIsobars; $i++) {
61 $this->isobarValues[$i] = $isobar;
62 $isobar += $stepSize;
63 }
64 }
65
66 if( $aColors !== null && count($aColors) > 0 ) {
67
68 if( !is_array($aColors) ) {
69 JpGraphError::RaiseL(28001);
70 //'Third argument to Contour must be an array of colors.'
71 }
72
73 if( count($aColors) != count($this->isobarValues) ) {
74 JpGraphError::RaiseL(28002);
75 //'Number of colors must equal the number of isobar lines specified';
76 }
77
78 $this->isobarColors = $aColors;
79 }
80 }
81
82 /**
83 * Flip the plot around the Y-coordinate. This has the same affect as flipping the input
84 * data matrice
85 *
86 * @param $aFlg If true the the vertice in input data matrice position (0,0) corresponds to the top left
87 * corner of teh plot otherwise it will correspond to the bottom left corner (a horizontal flip)
88 */
89 function SetInvert($aFlg=true) {
90 $this->invert = $aFlg;
91 }
92
93 /**
94 * Find the min and max values in the data matrice
95 *
96 * @return array(min_value,max_value)
97 */
98 function getMinMaxVal() {
99 $min = $this->dataPoints[0][0];
100 $max = $this->dataPoints[0][0];
101 for ($i = 0; $i < $this->nbrRows; $i++) {
102 if( ($mi=min($this->dataPoints[$i])) < $min ) $min = $mi;
103 if( ($ma=max($this->dataPoints[$i])) > $max ) $max = $ma;
104 }
105 return array($min,$max);
106 }
107
108 /**
109 * Reset the two matrices that keeps track on where the isobars crosses the
110 * horizontal and vertical edges
111 */
112 function resetEdgeMatrices() {
113 for ($k = 0; $k < 2; $k++) {
114 for ($i = 0; $i <= $this->nbrRows; $i++) {
115 for ($j = 0; $j <= $this->nbrCols; $j++) {
116 $this->edges[$k][$i][$j] = false;
117 }
118 }
119 }
120 }
121
122 /**
123 * Determine if the specified isobar crosses the horizontal edge specified by its row and column
124 *
125 * @param $aRow Row index of edge to be checked
126 * @param $aCol Col index of edge to be checked
127 * @param $aIsobar Isobar value
128 * @return true if the isobar is crossing this edge
129 */
130 function isobarHCrossing($aRow,$aCol,$aIsobar) {
131
132 if( $aCol >= $this->nbrCols-1 ) {
133 JpGraphError::RaiseL(28003,$aCol);
134 //'ContourPlot Internal Error: isobarHCrossing: Coloumn index too large (%d)'
135 }
136 if( $aRow >= $this->nbrRows ) {
137 JpGraphError::RaiseL(28004,$aRow);
138 //'ContourPlot Internal Error: isobarHCrossing: Row index too large (%d)'
139 }
140
141 $v1 = $this->dataPoints[$aRow][$aCol];
142 $v2 = $this->dataPoints[$aRow][$aCol+1];
143
144 return ($aIsobar-$v1)*($aIsobar-$v2) < 0 ;
145
146 }
147
148 /**
149 * Determine if the specified isobar crosses the vertical edge specified by its row and column
150 *
151 * @param $aRow Row index of edge to be checked
152 * @param $aCol Col index of edge to be checked
153 * @param $aIsobar Isobar value
154 * @return true if the isobar is crossing this edge
155 */
156 function isobarVCrossing($aRow,$aCol,$aIsobar) {
157
158 if( $aRow >= $this->nbrRows-1) {
159 JpGraphError::RaiseL(28005,$aRow);
160 //'isobarVCrossing: Row index too large
161 }
162 if( $aCol >= $this->nbrCols ) {
163 JpGraphError::RaiseL(28006,$aCol);
164 //'isobarVCrossing: Col index too large
165 }
166
167 $v1 = $this->dataPoints[$aRow][$aCol];
168 $v2 = $this->dataPoints[$aRow+1][$aCol];
169
170 return ($aIsobar-$v1)*($aIsobar-$v2) < 0 ;
171
172 }
173
174 /**
175 * Determine all edges, horizontal and vertical that the specified isobar crosses. The crossings
176 * are recorded in the two edge matrices.
177 *
178 * @param $aIsobar The value of the isobar to be checked
179 */
180 function determineIsobarEdgeCrossings($aIsobar) {
181
182 $ib = $this->isobarValues[$aIsobar];
183
184 for ($i = 0; $i < $this->nbrRows-1; $i++) {
185 for ($j = 0; $j < $this->nbrCols-1; $j++) {
186 $this->edges[HORIZ_EDGE][$i][$j] = $this->isobarHCrossing($i,$j,$ib);
187 $this->edges[VERT_EDGE][$i][$j] = $this->isobarVCrossing($i,$j,$ib);
188 }
189 }
190
191 // We now have the bottom and rightmost edges unsearched
192 for ($i = 0; $i < $this->nbrRows-1; $i++) {
193 $this->edges[VERT_EDGE][$i][$j] = $this->isobarVCrossing($i,$this->nbrCols-1,$ib);
194 }
195 for ($j = 0; $j < $this->nbrCols-1; $j++) {
196 $this->edges[HORIZ_EDGE][$i][$j] = $this->isobarHCrossing($this->nbrRows-1,$j,$ib);
197 }
198
199 }
200
201 /**
202 * Return the normalized coordinates for the crossing of the specified edge with the specified
203 * isobar- The crossing is simpy detrmined with a linear interpolation between the two vertices
204 * on each side of the edge and the value of the isobar
205 *
206 * @param $aRow Row of edge
207 * @param $aCol Column of edge
208 * @param $aEdgeDir Determine if this is a horizontal or vertical edge
209 * @param $ib The isobar value
210 * @return unknown_type
211 */
212 function getCrossingCoord($aRow,$aCol,$aEdgeDir,$aIsobarVal) {
213
214 // In order to avoid numerical problem when two vertices are very close
215 // we have to check and avoid dividing by close to zero denumerator.
216 if( $aEdgeDir == HORIZ_EDGE ) {
217 $d = abs($this->dataPoints[$aRow][$aCol] - $this->dataPoints[$aRow][$aCol+1]);
218 if( $d > 0.001 ) {
219 $xcoord = $aCol + abs($aIsobarVal - $this->dataPoints[$aRow][$aCol]) / $d;
220 }
221 else {
222 $xcoord = $aCol;
223 }
224 $ycoord = $aRow;
225 }
226 else {
227 $d = abs($this->dataPoints[$aRow][$aCol] - $this->dataPoints[$aRow+1][$aCol]);
228 if( $d > 0.001 ) {
229 $ycoord = $aRow + abs($aIsobarVal - $this->dataPoints[$aRow][$aCol]) / $d;
230 }
231 else {
232 $ycoord = $aRow;
233 }
234 $xcoord = $aCol;
235 }
236 if( $this->invert ) {
237 $ycoord = $this->nbrRows-1 - $ycoord;
238 }
239 return array($xcoord,$ycoord);
240
241 }
242
243 /**
244 * In order to avoid all kinds of unpleasent extra checks and complex boundary
245 * controls for the degenerated case where the contour levels exactly crosses
246 * one of the vertices we add a very small delta (0.1%) to the data point value.
247 * This has no visible affect but it makes the code sooooo much cleaner.
248 *
249 */
250 function adjustDataPointValues() {
251
252 $ni = count($this->isobarValues);
253 for ($k = 0; $k < $ni; $k++) {
254 $ib = $this->isobarValues[$k];
255 for ($row = 0 ; $row < $this->nbrRows-1; ++$row) {
256 for ($col = 0 ; $col < $this->nbrCols-1; ++$col ) {
257 if( abs($this->dataPoints[$row][$col] - $ib) < 0.0001 ) {
258 $this->dataPoints[$row][$col] += $this->dataPoints[$row][$col]*0.001;
259 }
260 }
261 }
262 }
263
264 }
265
266 /**
267 * @param $aFlg
268 * @param $aBW
269 * @return unknown_type
270 */
271 function UseHighContrastColor($aFlg=true,$aBW=false) {
272 $this->highcontrast = $aFlg;
273 $this->highcontrastbw = $aBW;
274 }
275
276 /**
277 * Calculate suitable colors for each defined isobar
278 *
279 */
280 function CalculateColors() {
281 if ( $this->highcontrast ) {
282 if ( $this->highcontrastbw ) {
283 for ($ib = 0; $ib < $this->nbrIsobars; $ib++) {
284 $this->isobarColors[$ib] = 'black';
285 }
286 }
287 else {
288 // Use only blue/red scale
289 $step = round(255/($this->nbrIsobars-1));
290 for ($ib = 0; $ib < $this->nbrIsobars; $ib++) {
291 $this->isobarColors[$ib] = array($ib*$step, 50, 255-$ib*$step);
292 }
293 }
294 }
295 else {
296 $n = $this->nbrIsobars;
297 $v = 0; $step = 1 / ($this->nbrIsobars-1);
298 for ($ib = 0; $ib < $this->nbrIsobars; $ib++) {
299 $this->isobarColors[$ib] = RGB::GetSpectrum($v);
300 $v += $step;
301 }
302 }
303 }
304
305 /**
306 * This is where the main work is done. For each isobar the crossing of the edges are determined
307 * and then each cell is analyzed to find the 0, 2 or 4 crossings. Then the normalized coordinate
308 * for the crossings are determined and pushed on to the isobar stack. When the method is finished
309 * the $isobarCoord will hold one arrayfor each isobar where all the line segments that makes
310 * up the contour plot are stored.
311 *
312 * @return array( $isobarCoord, $isobarValues, $isobarColors )
313 */
314 function getIsobars() {
315
316 $this->adjustDataPointValues();
317
318 for ($isobar = 0; $isobar < $this->nbrIsobars; $isobar++) {
319
320 $ib = $this->isobarValues[$isobar];
321 $this->resetEdgeMatrices();
322 $this->determineIsobarEdgeCrossings($isobar);
323 $this->isobarCoord[$isobar] = array();
324
325 $ncoord = 0;
326
327 for ($row = 0 ; $row < $this->nbrRows-1; ++$row) {
328 for ($col = 0 ; $col < $this->nbrCols-1; ++$col ) {
329
330 // Find out how many crossings around the edges
331 $n = 0;
332 if ( $this->edges[HORIZ_EDGE][$row][$col] ) $neigh[$n++] = array($row, $col, HORIZ_EDGE);
333 if ( $this->edges[HORIZ_EDGE][$row+1][$col] ) $neigh[$n++] = array($row+1,$col, HORIZ_EDGE);
334 if ( $this->edges[VERT_EDGE][$row][$col] ) $neigh[$n++] = array($row, $col, VERT_EDGE);
335 if ( $this->edges[VERT_EDGE][$row][$col+1] ) $neigh[$n++] = array($row, $col+1,VERT_EDGE);
336
337 if ( $n == 2 ) {
338 $n1=0; $n2=1;
339 $this->isobarCoord[$isobar][$ncoord++] = array(
340 $this->getCrossingCoord($neigh[$n1][0],$neigh[$n1][1],$neigh[$n1][2],$ib),
341 $this->getCrossingCoord($neigh[$n2][0],$neigh[$n2][1],$neigh[$n2][2],$ib) );
342 }
343 elseif ( $n == 4 ) {
344 // We must determine how to connect the edges either northwest->southeast or
345 // northeast->southwest. We do that by calculating the imaginary middle value of
346 // the cell by averaging the for corners. This will compared with the value of the
347 // top left corner will help determine the orientation of the ridge/creek
348 $midval = ($this->dataPoints[$row][$col]+$this->dataPoints[$row][$col+1]+$this->dataPoints[$row+1][$col]+$this->dataPoints[$row+1][$col+1])/4;
349 $v = $this->dataPoints[$row][$col];
350 if( $midval == $ib ) {
351 // Orientation "+"
352 $n1=0; $n2=1; $n3=2; $n4=3;
353 } elseif ( ($midval > $ib && $v > $ib) || ($midval < $ib && $v < $ib) ) {
354 // Orientation of ridge/valley = "\"
355 $n1=0; $n2=3; $n3=2; $n4=1;
356 } elseif ( ($midval > $ib && $v < $ib) || ($midval < $ib && $v > $ib) ) {
357 // Orientation of ridge/valley = "/"
358 $n1=0; $n2=2; $n3=3; $n4=1;
359 }
360
361 $this->isobarCoord[$isobar][$ncoord++] = array(
362 $this->getCrossingCoord($neigh[$n1][0],$neigh[$n1][1],$neigh[$n1][2],$ib),
363 $this->getCrossingCoord($neigh[$n2][0],$neigh[$n2][1],$neigh[$n2][2],$ib) );
364
365 $this->isobarCoord[$isobar][$ncoord++] = array(
366 $this->getCrossingCoord($neigh[$n3][0],$neigh[$n3][1],$neigh[$n3][2],$ib),
367 $this->getCrossingCoord($neigh[$n4][0],$neigh[$n4][1],$neigh[$n4][2],$ib) );
368
369 }
370 }
371 }
372 }
373
374 if( count($this->isobarColors) == 0 ) {
375 // No manually specified colors. Calculate them automatically.
376 $this->CalculateColors();
377 }
378 return array( $this->isobarCoord, $this->isobarValues, $this->isobarColors );
379 }
380}
381
382
383/**
384 * This class represent a plotting of a contour outline of data given as a X-Y matrice
385 *
386 */
387class ContourPlot extends Plot {
388
389 private $contour, $contourCoord, $contourVal, $contourColor;
390 private $nbrCountours = 0 ;
391 private $dataMatrix = array();
392 private $invertLegend = false;
393 private $interpFactor = 1;
394 private $flipData = false;
395 private $isobar = 10;
396 private $showLegend = false;
397 private $highcontrast = false, $highcontrastbw = false;
398 private $manualIsobarColors = array();
399
400 /**
401 * Construct a contour plotting algorithm. The end result of the algorithm is a sequence of
402 * line segments for each isobar given as two vertices.
403 *
404 * @param $aDataMatrix The Z-data to be used
405 * @param $aIsobar A mixed variable, if it is an integer then this specified the number of isobars to use.
406 * The values of the isobars are automatically detrmined to be equ-spaced between the min/max value of the
407 * data. If it is an array then it explicetely gives the isobar values
408 * @param $aInvert By default the matrice with row index 0 corresponds to Y-value 0, i.e. in the bottom of
409 * the plot. If this argument is true then the row with the highest index in the matrice corresponds to
410 * Y-value 0. In affect flipping the matrice around an imaginary horizontal axis.
411 * @param $aHighContrast Use high contrast colors (blue/red:ish)
412 * @param $aHighContrastBW Use only black colors for contours
413 * @return an instance of the contour plot algorithm
414 */
415 function __construct($aDataMatrix, $aIsobar=10, $aFactor=1, $aInvert=false, $aIsobarColors=array()) {
416
417 $this->dataMatrix = $aDataMatrix;
418 $this->flipData = $aInvert;
419 $this->isobar = $aIsobar;
420 $this->interpFactor = $aFactor;
421
422 if ( $this->interpFactor > 1 ) {
423
424 if( $this->interpFactor > 5 ) {
425 JpGraphError::RaiseL(28007);// ContourPlot interpolation factor is too large (>5)
426 }
427
428 $ip = new MeshInterpolate();
429 $this->dataMatrix = $ip->Linear($this->dataMatrix, $this->interpFactor);
430 }
431
432 $this->contour = new Contour($this->dataMatrix,$this->isobar,$aIsobarColors);
433
434 if( is_array($aIsobar) )
435 $this->nbrContours = count($aIsobar);
436 else
437 $this->nbrContours = $aIsobar;
438 }
439
440
441 /**
442 * Flipe the data around the center
443 *
444 * @param $aFlg
445 *
446 */
447 function SetInvert($aFlg=true) {
448 $this->flipData = $aFlg;
449 }
450
451 /**
452 * Set the colors for the isobar lines
453 *
454 * @param $aColorArray
455 *
456 */
457 function SetIsobarColors($aColorArray) {
458 $this->manualIsobarColors = $aColorArray;
459 }
460
461 /**
462 * Show the legend
463 *
464 * @param $aFlg true if the legend should be shown
465 *
466 */
467 function ShowLegend($aFlg=true) {
468 $this->showLegend = $aFlg;
469 }
470
471
472 /**
473 * @param $aFlg true if the legend should start with the lowest isobar on top
474 * @return unknown_type
475 */
476 function Invertlegend($aFlg=true) {
477 $this->invertLegend = $aFlg;
478 }
479
480 /* Internal method. Give the min value to be used for the scaling
481 *
482 */
483 function Min() {
484 return array(0,0);
485 }
486
487 /* Internal method. Give the max value to be used for the scaling
488 *
489 */
490 function Max() {
491 return array(count($this->dataMatrix[0])-1,count($this->dataMatrix)-1);
492 }
493
494 /**
495 * Internal ramewrok method to setup the legend to be used for this plot.
496 * @param $aGraph The parent graph class
497 */
498 function Legend($aGraph) {
499
500 if( ! $this->showLegend )
501 return;
502
503 if( $this->invertLegend ) {
504 for ($i = 0; $i < $this->nbrContours; $i++) {
505 $aGraph->legend->Add(sprintf('%.1f',$this->contourVal[$i]), $this->contourColor[$i]);
506 }
507 }
508 else {
509 for ($i = $this->nbrContours-1; $i >= 0 ; $i--) {
510 $aGraph->legend->Add(sprintf('%.1f',$this->contourVal[$i]), $this->contourColor[$i]);
511 }
512 }
513 }
514
515
516 /**
517 * Framework function which gets called before the Stroke() method is called
518 *
519 * @see Plot#PreScaleSetup($aGraph)
520 *
521 */
522 function PreScaleSetup($aGraph) {
523 $xn = count($this->dataMatrix[0])-1;
524 $yn = count($this->dataMatrix)-1;
525
526 $aGraph->xaxis->scale->Update($aGraph->img,0,$xn);
527 $aGraph->yaxis->scale->Update($aGraph->img,0,$yn);
528
529 $this->contour->SetInvert($this->flipData);
530 list($this->contourCoord,$this->contourVal,$this->contourColor) = $this->contour->getIsobars();
531 }
532
533 /**
534 * Use high contrast color schema
535 *
536 * @param $aFlg True, to use high contrast color
537 * @param $aBW True, Use only black and white color schema
538 */
539 function UseHighContrastColor($aFlg=true,$aBW=false) {
540 $this->highcontrast = $aFlg;
541 $this->highcontrastbw = $aBW;
542 $this->contour->UseHighContrastColor($this->highcontrast,$this->highcontrastbw);
543 }
544
545 /**
546 * Internal method. Stroke the contour plot to the graph
547 *
548 * @param $img Image handler
549 * @param $xscale Instance of the xscale to use
550 * @param $yscale Instance of the yscale to use
551 */
552 function Stroke($img,$xscale,$yscale) {
553
554 if( count($this->manualIsobarColors) > 0 ) {
555 $this->contourColor = $this->manualIsobarColors;
556 if( count($this->manualIsobarColors) != $this->nbrContours ) {
557 JpGraphError::RaiseL(28002);
558 }
559 }
560
561 $img->SetLineWeight($this->line_weight);
562
563 for ($c = 0; $c < $this->nbrContours; $c++) {
564
565 $img->SetColor( $this->contourColor[$c] );
566
567 $n = count($this->contourCoord[$c]);
568 $i = 0;
569 while ( $i < $n ) {
570 list($x1,$y1) = $this->contourCoord[$c][$i][0];
571 $x1t = $xscale->Translate($x1);
572 $y1t = $yscale->Translate($y1);
573
574 list($x2,$y2) = $this->contourCoord[$c][$i++][1];
575 $x2t = $xscale->Translate($x2);
576 $y2t = $yscale->Translate($y2);
577
578 $img->Line($x1t,$y1t,$x2t,$y2t);
579 }
580
581 }
582 }
583
584}
585
586// EOF
587?>
Note: See TracBrowser for help on using the repository browser.