Bus Stop density heatmap
Bus Stop density heatmap

--- a/include/common-session.inc.php
+++ b/include/common-session.inc.php
@@ -51,7 +51,7 @@
 	session_destroy();
 	session_start();
 }
-debug(print_r($_SESSION, true) , "session");
+//debug(print_r($_SESSION, true) , "session");
 
 function current_time() {
 	return ($_SESSION['time']? $_SESSION['time'] : date("H:i:s"));

--- /dev/null
+++ b/labs/busstopdensity.php
@@ -1,1 +1,77 @@
+<?php
+include ('../include/common.inc.php');
+//include_header("Bus Stop Density", "busstopdensity")
+?>
 
+<style type="text/css">
+
+			#map_container{
+				width:100%;
+				height:100%;
+			}
+			#status, #legend {
+				background-color: #eeeeee;
+				padding: 3px 5px 3px 5px;
+				opacity: .9;
+			}
+		</style>
+	
+		<div id="map_container"></div>
+			<div id="status">
+			Status:&nbsp; <span id="log"> <img alt="progess bar" src="progress_bar.gif" width="150" height="16"/></span>
+		</div>
+		<script type="text/javascript" src="http://www.google.com/jsapi?autoload={%22modules%22:[{%22name%22:%22maps%22,version:3,other_params:%22sensor=false%22},{%22name%22:%22jquery%22,%22version%22:%221.4.2%22}]}"></script>
+		<script type="text/javascript">
+		//<![CDATA[
+		//Google Map API v3
+		
+		var googleMap = null;
+		var previousPos = null;
+
+		$(function($){//Called when page is loaded
+			googleMap = new google.maps.Map(document.getElementById("map_container"), {
+				zoom: 17, 
+				center: new google.maps.LatLng(-35.25,149.125), 
+				mapTypeId: google.maps.MapTypeId.SATELLITE});
+			//Set status bar
+			googleMap.controls[google.maps.ControlPosition.TOP_LEFT].push($("#status").get(0));
+			//Set legend
+			googleMap.controls[google.maps.ControlPosition.BOTTOM_LEFT].push($("#legend").get(0));
+
+			google.maps.event.addListener(googleMap, "zoom_changed", function(){
+				google.maps.event.trigger(googleMap, "mousemove", previousPos);
+			});//onzoomend
+				
+			//Add a listener when mouse moves
+			google.maps.event.addListener(googleMap, "mousemove", function(event){
+				var latLng = event.latLng;
+				var xy = googleMap.getProjection().fromLatLngToPoint(latLng);
+				var ratio = Math.pow(2,googleMap.getZoom());
+				$("#log").html("Zoom:" + googleMap.getZoom() + " WGS84:(" + latLng.lat().toFixed(5) + ", " + latLng.lng().toFixed(5) + ") Px:(" + Math.floor(xy.x * ratio)  + "," + Math.floor(xy.y *ratio) + ")");
+				previousPos = event;
+			});//onmouseover
+
+			//Add a listener when mouse leaves the map area
+			google.maps.event.addListener(googleMap, "mouseout", function(event){
+				$("#log").html("");
+			});//onmouseout
+
+			//Add tile overlay
+			var myOverlay = new google.maps.ImageMapType({
+				getTileUrl: function(coord, zoom) {
+					return 'busstopdensity.tile.php?x=' + coord.x + '&y=' + coord.y + '&zoom=' +zoom;
+				},
+				tileSize: new google.maps.Size(256, 256),
+				isPng: true,
+				opacity:1.0
+			});
+			googleMap.overlayMapTypes.insertAt(0, myOverlay);
+			$("#log").html("Map loaded!");
+		});//onload
+		//]]>
+		</script>
+<?php
+include_footer()
+?>
+        
+

--- /dev/null
+++ b/labs/busstopdensity.tile.php
@@ -1,1 +1,121 @@
+<?php
+include ('../include/common.inc.php');
+$debugOkay = Array();
 
+/*
+*DISCLAIMER
+*  http://blog.gmapify.fr/create-beautiful-tiled-heat-maps-with-php-and-gd
+*THIS SOFTWARE IS PROVIDED BY THE AUTHOR 'AS IS' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES *OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, *INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF *USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT *(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*
+*	@author: Olivier G. <olbibigo_AT_gmail_DOT_com>
+*	@version: 1.0
+*	@history:
+*		1.0	creation
+*/
+	set_time_limit(120);//2mn
+	ini_set('memory_limit', '256M');
+error_reporting(E_ALL ^ E_DEPRECATED);
+	require_once ('lib/GoogleMapUtility.php');
+	require_once ('lib/HeatMap.php');
+
+	//Root folder to store generated tiles
+	define('TILE_DIR', 'tiles/');
+	//Covered geographic areas
+	define('MIN_LAT', -35.48);
+	define('MAX_LAT', -35.15);
+	define('MIN_LNG', 148.98);
+	define('MAX_LNG', 149.25);
+	define('TILE_SIZE_FACTOR', 0.5);
+	define('SPOT_RADIUS', 30);
+	define('SPOT_DIMMING_LEVEL', 50);
+	
+	//Input parameters
+	if(isset($_GET['x']))
+		$X = (int)$_GET['x'];
+	else
+		exit("x missing");
+	if(isset($_GET['y']))
+		$Y = (int)$_GET['y'];
+	else
+		exit("y missing");
+	if(isset($_GET['zoom']))
+		$zoom = (int)$_GET['zoom'];
+	else
+		exit("zoom missing");
+
+	$dir = TILE_DIR.$zoom;
+	$tilename = $dir.'/'.$X.'_'.$Y.'.png';
+	//HTTP headers  (data type and caching rule)
+	header("Cache-Control: must-revalidate");
+	header("Expires: " . gmdate("D, d M Y H:i:s", time() + 86400) . " GMT");
+	if(!file_exists($tilename)){
+		$rect = GoogleMapUtility::getTileRect($X, $Y, $zoom);
+		//A tile can contain part of a spot with center in an adjacent tile (overlaps).
+		//Knowing the spot radius (in pixels) and zoom level, a smart way to process tiles would be to compute the box (in decimal degrees) containing only spots that can be drawn on current tile. We choose a simpler solution by increeasing  geo bounds by 2*TILE_SIZE_FACTOR whatever the zoom level and spot radius.
+		$extend_X = $rect->width * TILE_SIZE_FACTOR;//in decimal degrees
+		$extend_Y = $rect->height * TILE_SIZE_FACTOR;//in decimal degrees
+		$swlat = $rect->y - $extend_Y;
+		$swlng = $rect->x - $extend_X;
+		$nelat = $swlat + $rect->height + 2 * $extend_Y;
+		$nelng = $swlng + $rect->width + 2 * $extend_X;
+
+		if( ($nelat <= MIN_LAT) || ($swlat >= MAX_LAT) || ($nelng <= MIN_LNG) || ($swlng >= MAX_LNG)){
+			//No geodata so return generic empty tile
+			echo file_get_contents(TILE_DIR.'empty.png');
+			exit();
+		}
+
+		//Get McDonald's spots
+		$spots = fGetPOI('Select * from stops where
+				(stop_lon > '.$swlng.' AND stop_lon < '.$nelng.')
+			AND (stop_lat < '.$nelat.' AND stop_lat > '.$swlat.')', $im, $X, $Y, $zoom, SPOT_RADIUS);
+
+		
+		if(empty($spots)){
+			//No geodata so return generic empty tile
+			header('Content-type: image/png');
+			echo file_get_contents(TILE_DIR.'empty.png');
+		}else{
+			if(!file_exists($dir)){
+				mkdir($dir, 0705);
+			}
+			//All the magics is in HeatMap class :)
+			$im = HeatMap::createImage($spots, GoogleMapUtility::TILE_SIZE, GoogleMapUtility::TILE_SIZE, heatMap::$WITH_ALPHA, SPOT_RADIUS, SPOT_DIMMING_LEVEL, HeatMap::$GRADIENT_FIRE);
+			//Store tile for reuse and output it
+			header('content-type:image/png;');
+			imagepng($im, $tilename);
+			echo file_get_contents($tilename);
+			imagedestroy($im);
+			unset($im);
+		}
+	}else{
+		//Output stored tile
+		header('content-type:image/png;');
+		echo file_get_contents($tilename);
+	}
+	/////////////
+	//Functions//
+	/////////////
+	function fGetPOI($query, &$im, $X, $Y, $zoom, $offset){
+            global $conn;
+		$nbPOIInsideTile = 0;
+
+	$result = pg_query($conn, $query);
+        $spots = Array();
+	if (!$result) {
+		databaseError(pg_result_error($result));
+		return Array();
+	}
+	foreach( pg_fetch_all($result) as $row){
+				$point = GoogleMapUtility::getOffsetPixelCoords($row['stop_lat'], $row['stop_lon'], $zoom, $X, $Y);
+				//Count result only in the tile
+				if( ($point->x > -$offset) && ($point->x < (GoogleMapUtility::TILE_SIZE+$offset)) && ($point->y > -$offset) && ($point->y < (GoogleMapUtility::TILE_SIZE+$offset))){
+					$spots[] = new HeatMapPoint($point->x, $point->y);
+				}
+				
+			}//while
+		return $spots;
+	}//fAddPOI
+?>
+
+

 Binary files /dev/null and b/labs/gradient/classic.png differ
 Binary files /dev/null and b/labs/gradient/fire.png differ
 Binary files /dev/null and b/labs/gradient/pgaitch.png differ
--- a/labs/index.php
+++ b/labs/index.php
@@ -6,8 +6,10 @@
 		<li data-role="list-divider" > Experimental Features </li>
 		<li><a href="mywaybalance.php"><h3>MyWay Balance for mobile</h3>
 		<p>Mobile viewer for MyWay balance. Warning! No HTTPS security.</p></a></li>
-		<li><a href="networkstats.php"><h3>Network/Route Statistics</h3>
-		<p>Statistical analysis of network/routes</p></a></li>
+		<li><a href="networkstats.php"><h3>Route Statistics</h3>
+		<p>Analysis of route timing points</p></a></li>
+		<li><a href="busstopdensity.php"><h3>Bus Stop Density Map</h3>
+		<p>Analysis of bus stop coverage</p></a></li>
 		<li>More coming soon!</li>
             </ul>
 	    </div>

--- a/labs/networkstats.php
+++ b/labs/networkstats.php
@@ -1,6 +1,6 @@
 <?php
 include ('../include/common.inc.php');
-include_header("Network/Route Statistics", "networkstats")
+include_header("Route Statistics", "networkstats")
 ?>
 <script type="text/javascript" src="js/flotr/lib/prototype-1.6.0.2.js"></script>
 

 Binary files /dev/null and b/labs/tiles/empty.png differ
--- /dev/null
+++ b/lib/GoogleMapUtility.php
@@ -1,1 +1,97 @@
+<?php
+/*
+*DISCLAIMER
+* 
+*THIS SOFTWARE IS PROVIDED BY THE AUTHOR 'AS IS' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES *OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, *INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF *USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT *(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*
+*	@author: Olivier G. <olbibigo_AT_gmail_DOT_com>
+*	@version: 1.1
+*	@history:
+*		1.0	creation
+		1.1	disclaimer added
+*/
+class GoogleMapUtility {
+	const TILE_SIZE = 256;
+	
+	//(lat, lng, z) -> parent tile (X,Y)
+	public static function getTileXY($lat, $lng, $zoom) {
+		$normalised = GoogleMapUtility::_toNormalisedMercatorCoords(GoogleMapUtility::_toMercatorCoords($lat, $lng));
+		$scale = 1 << ($zoom);
+		return new Point(
+			(int)($normalised->x * $scale), 
+			(int)($normalised->y * $scale)
+		);
+	}//toTileXY
+	
+	//(lat, lng, z) -> (x,y) with (0,0) in the upper left corner of the MAP
+	public static function getPixelCoords($lat, $lng, $zoom) {
+		$normalised = GoogleMapUtility::_toNormalisedMercatorCoords(GoogleMapUtility::_toMercatorCoords($lat, $lng));
+		$scale = (1 << ($zoom)) * GoogleMapUtility::TILE_SIZE;
+		return new Point(
+			(int)($normalised->x * $scale), 
+			(int)($normalised->y * $scale)
+		);
+	}//getPixelCoords
 
+	//(lat, lng, z) -> (x,y) in the upper left corner of the TILE ($X, $Y)
+	public static function getOffsetPixelCoords($lat,$lng,$zoom, $X, $Y) {
+		$pixelCoords = GoogleMapUtility::getPixelCoords($lat, $lng, $zoom);
+		return new Point(
+			$pixelCoords->x - $X * GoogleMapUtility::TILE_SIZE, 
+			$pixelCoords->y - $Y * GoogleMapUtility::TILE_SIZE
+		);
+	}//getPixelOffsetInTile
+	
+	public static function getTileRect($X,$Y,$zoom) {
+		$tilesAtThisZoom = 1 << $zoom;
+		$lngWidth = 360.0 / $tilesAtThisZoom;
+		$lng = -180 + ($X * $lngWidth);	
+		$latHeightMerc = 1.0 / $tilesAtThisZoom;
+		$topLatMerc = $Y * $latHeightMerc;
+		$bottomLatMerc = $topLatMerc + $latHeightMerc;
+		$bottomLat = (180 / M_PI) * ((2 * atan(exp(M_PI * (1 - (2 * $bottomLatMerc))))) - (M_PI / 2));
+		$topLat = (180 / M_PI) * ((2 * atan(exp(M_PI * (1 - (2 * $topLatMerc))))) - (M_PI / 2));
+		$latHeight = $topLat - $bottomLat;
+		return new Boundary($lng, $bottomLat, $lngWidth, $latHeight);
+	}//getTileRect	
+
+	private static function _toMercatorCoords($lat, $lng) {
+		if ($lng > 180) {
+			$lng -= 360;
+		}
+		$lng /= 360;
+		$lat = asinh(tan(deg2rad($lat)))/M_PI/2;
+		return new Point($lng, $lat);
+	}//_toMercatorCoords
+
+	private static function _toNormalisedMercatorCoords($point) {
+		$point->x += 0.5;
+		$point->y = abs($point->y-0.5);
+		return $point;
+	}//_toNormalisedMercatorCoords
+}//GoogleMapUtility
+
+class Point {
+	public $x,$y;
+	function __construct($x,$y) {
+		$this->x = $x;
+		$this->y = $y;
+	}
+	function __toString() {
+		return "({$this->x},{$this->y})";
+	}
+}//Point
+
+class Boundary {
+	public $x,$y,$width,$height;
+	function __construct($x,$y,$width,$height) {
+		$this->x = $x;
+		$this->y = $y;
+		$this->width = $width;
+		$this->height = $height;
+	}
+	function __toString() {
+		return "({$this->x} x {$this->y},{$this->width},{$this->height})";
+	}
+}//Boundary
+?>

file:b/lib/HeatMap.php (new)
--- /dev/null
+++ b/lib/HeatMap.php
@@ -1,1 +1,275 @@
-
+<?php
+/*
+*DISCLAIMER
+* http://blog.gmapify.fr/create-beautiful-tiled-heat-maps-with-php-and-gd
+*THIS SOFTWARE IS PROVIDED BY THE AUTHOR 'AS IS' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES *OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, *INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF *USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT *(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*
+*	@author: Olivier G. <olbibigo_AT_gmail_DOT_com>
+*	@version: 1.0
+*	@history:
+*		1.0	creation
+*/
+	define('PI2', 2*M_PI);
+
+	class HeatMapPoint{
+		public $x,$y;
+		function __construct($x,$y) {
+			$this->x = $x;
+			$this->y = $y;
+		}
+		function __toString() {
+			return "({$this->x},{$this->y})";
+		}
+	}//Point
+
+	class HeatMap{
+		//TRANSPARENCY
+		public static $WITH_ALPHA = 0;
+		public static $WITH_TRANSPARENCY = 1;
+		//GRADIENT STYLE
+		public static $GRADIENT_CLASSIC = 'classic';
+		public static $GRADIENT_FIRE = 'fire';
+		public static $GRADIENT_PGAITCH = 'pgaitch';
+		//GRADIENT MODE (for heatImage)
+		public static $GRADIENT_NO_NEGATE_NO_INTERPOLATE = 0;
+		public static $GRADIENT_NO_NEGATE_INTERPOLATE = 1;
+		public static $GRADIENT_NEGATE_NO_INTERPOLATE = 2;
+		public static $GRADIENT_NEGATE_INTERPOLATE = 3;
+		//NOT PROCESSED PIXEL (for heatImage)
+		public static $KEEP_VALUE = 0;
+		public static $NO_KEEP_VALUE = 1;
+		//CONSTRAINTS
+		private static $MIN_RADIUS = 2;//in px
+		private static $MAX_RADIUS = 50;//in px
+		private static $MAX_IMAGE_SIZE = 10000;//in px
+		
+		//generate an $image_width by $image_height pixels heatmap image of $points
+		public static function createImage($data, $image_width, $image_height, $mode=0, $spot_radius = 30, $dimming = 75, $gradient_name = 'classic'){
+			$_gradient_name = $gradient_name;
+			if(($_gradient_name != self::$GRADIENT_CLASSIC) && ($_gradient_name != self::$GRADIENT_FIRE) && ($_gradient_name != self::$GRADIENT_PGAITCH)){
+				$_gradient_name = self::$GRADIENT_CLASSIC;
+			}
+			$_image_width = min(self::$MAX_IMAGE_SIZE, max(0, intval($image_width)));
+			$_image_height = min(self::$MAX_IMAGE_SIZE, max(0, intval($image_height)));
+			$_spot_radius = min(self::$MAX_RADIUS, max(self::$MIN_RADIUS, intval($spot_radius)));
+			$_dimming = min(255, max(0, intval($dimming)));
+			if(!is_array($data)){
+				return false;
+			}
+			$im = imagecreatetruecolor($_image_width, $_image_height);
+			$white = imagecolorallocate($im, 255, 255, 255);
+			imagefill($im, 0, 0, $white);
+			if(self::$WITH_ALPHA == $mode){
+				imagealphablending($im, false);
+				imagesavealpha($im,true);
+			}
+			//Step 1: create grayscale image
+			foreach($data as $datum){
+				if( (is_array($datum) && (count($datum)==1)) || (!is_array($datum) && ('HeatMapPoint' == get_class($datum)))){//Plot points
+					if('HeatMapPoint' != get_class($datum)){
+						$datum = $datum[0];
+					}
+					self::_drawCircularGradient($im, $datum->x, $datum->y, $_spot_radius, $_dimming);
+				}else if(is_array($datum)){//Draw lines
+					$length = count($datum)-1;
+					for($i=0; $i < $length; ++$i){//Loop through points
+						//Bresenham's algorithm to plot from from $datum[$i] to $datum[$i+1];
+						self::_drawBilinearGradient($im, $datum[$i], $datum[$i+1], $_spot_radius, $_dimming);
+					}
+				}
+			}
+			//Gaussian filter
+			if($_spot_radius >= 30){
+				imagefilter($im, IMG_FILTER_GAUSSIAN_BLUR);
+			}
+			//Step 2: create colored image
+			if(FALSE === ($grad_rgba = self::_createGradient($im, $mode, $_gradient_name))){
+				return FALSE;
+			}
+			$grad_size = count($grad_rgba);
+			for($x=0; $x <$_image_width; ++$x){
+				for($y=0; $y <$_image_height; ++$y){
+					$level = imagecolorat($im, $x, $y) & 0xFF;
+					if( ($level >= 0) && ($level < $grad_size) ){
+						imagesetpixel($im, $x, $y, $grad_rgba[imagecolorat($im, $x, $y) & 0xFF]);
+					}
+				}
+			}
+			if(self::$WITH_TRANSPARENCY == $mode){
+				imagecolortransparent($im, $grad_rgba[count($grad_rgba)-1]);
+			}
+			return $im;
+		}//createImage
+
+		//Heat an image
+		public static function heatImage($filepath, $gradient_name = 'classic', $mode= 0, $min_level=0, $max_level=255, $gradient_interpolate=0, $keep_value=0){
+			$_gradient_name = $gradient_name;
+			if(($_gradient_name != self::$GRADIENT_CLASSIC) && ($_gradient_name != self::$GRADIENT_FIRE) && ($_gradient_name != self::$GRADIENT_PGAITCH)){
+				$_gradient_name = self::$GRADIENT_CLASSIC;
+			}
+			$_min_level = min(255, max(0, intval($min_level)));
+			$_max_level = min(255, max(0, intval($max_level)));
+
+			//try opening jpg first then png then gif format
+			if(FALSE === ($im = @imagecreatefromjpeg($filepath))){
+				if(FALSE === ($im = @imagecreatefrompng($filepath))){
+					if(FALSE === ($im = @imagecreatefromgif($filepath))){
+						return FALSE;
+					}
+				}
+			}
+			if(self::$WITH_ALPHA == $mode){
+				imagealphablending($im, false);
+				imagesavealpha($im,true);
+			}
+			$width = imagesx($im);
+			$height = imagesy($im);	
+			if(FALSE === ($grad_rgba = self::_createGradient($im, $mode, $_gradient_name))){
+				return FALSE;
+			}
+			//Convert to grayscale
+			$grad_size = count($grad_rgba);
+			$level_range = $_max_level - $_min_level;
+			for($x=0; $x <$width; ++$x){
+				for($y=0; $y <$height; ++$y){
+					$rgb = imagecolorat($im, $x, $y);
+					$r = ($rgb >> 16) & 0xFF;
+					$g = ($rgb >> 8) & 0xFF;
+					$b = $rgb & 0xFF;
+					$gray_level = Min(255, Max(0, floor(0.33 * $r + 0.5 * $g + 0.16 * $b)));//between 0 and 255				
+					if( ($gray_level >= $_min_level) && ($gray_level <= $_max_level) ){
+						switch($gradient_interpolate){
+							case self::$GRADIENT_NO_NEGATE_NO_INTERPOLATE:
+								//$_max_level takes related lowest gradient color
+								//$_min_level takes related highest gradient color
+								$value = 255 - $gray_level;
+								break;
+							case self::$GRADIENT_NEGATE_NO_INTERPOLATE:
+								//$_max_level takes related highest gradient color
+								//$_min_level takes related lowest gradient color
+								$value = $gray_level;
+								break;
+							case self::$GRADIENT_NO_NEGATE_INTERPOLATE:
+								//$_max_level takes lowest gradient color
+								//$_min_level takes highest gradient color
+								$value = 255- floor(($gray_level - $_min_level) * $grad_size / $level_range);
+								break;
+							case self::$GRADIENT_NEGATE_INTERPOLATE:
+								//$_max_level takes highest gradient color
+								//$_min_level takes lowest gradient color
+								$value = floor(($gray_level - $_min_level) * $grad_size / $level_range);
+								break;
+							default:
+						}
+						imagesetpixel($im, $x, $y, $grad_rgba[$value]);
+					}else{
+						if(self::$KEEP_VALUE == $keep_value){
+							//Do nothing
+						}else{//self::$NO_KEEP_VALUE
+							imagesetpixel($im, $x, $y, imagecolorallocatealpha($im,0,0,0,0));
+						}
+					}
+				}
+			}			
+			if(self::$WITH_TRANSPARENCY == $mode){
+				imagecolortransparent($im, $grad_rgba[count($grad_rgba)-1]);
+			}
+			return $im;
+		}//heatImage
+		
+		private static function _drawCircularGradient(&$im, $center_x, $center_y, $spot_radius, $dimming){
+			$dirty = array();
+			$ratio = (255 - $dimming) / $spot_radius;
+			for($r=$spot_radius; $r > 0; --$r){
+				$channel = $dimming + $r * $ratio;
+				$angle_step = 0.45/$r; //0.01;
+				//Process pixel by pixel to draw a radial grayscale radient
+				for($angle=0; $angle <= PI2; $angle += $angle_step){
+					$x = floor($center_x + $r*cos($angle));
+					$y = floor($center_y + $r*sin($angle));
+					if(!isset($dirty[$x][$y])){
+						$previous_channel = @imagecolorat($im, $x, $y) & 0xFF;//grayscale so same value
+						$new_channel = Max(0, Min(255,($previous_channel * $channel)/255));
+						imagesetpixel($im, $x, $y, imagecolorallocate($im, $new_channel, $new_channel, $new_channel));
+						$dirty[$x][$y] = 0;
+					}
+				}
+			}
+		}//_drawCircularGradient
+		
+		private static function _drawBilinearGradient(&$im, $point0, $point1, $spot_radius, $dimming){
+			if($point0->x < $point1->x){
+				$x0 = $point0->x;
+				$y0 = $point0->y;
+				$x1 = $point1->x;
+				$y1 = $point1->y;
+			}else{
+				$x0 = $point1->x;
+				$y0 = $point1->y;
+				$x1 = $point0->x;
+				$y1 = $point0->y;
+			}
+
+			if( ($x0==$x1) && ($y0==$y1)){//check if same coordinates
+				return false;
+			}
+			$steep = (abs($y1 - $y0) > abs($x1 - $x0))? true: false;
+			if($steep){
+				list($x0, $y0) = array($y0, $x0);//swap
+				list($x1, $y1) = array($y1, $x1);//swap
+			}
+			if($x0>$x1){
+				list($x0, $x1) = array($x1, $x0);//swap
+				list($y0, $y1) = array($y1, $y0);//swap
+			}
+			$deltax = $x1 - $x0;
+			$deltay = abs($y1 - $y0);
+			$error = $deltax / 2;
+			$y = $y0;
+			if( $y0 < $y1){
+				$ystep = 1; 
+			}else{
+				$ystep = -1;
+			}
+			$step = max(1, floor($spot_radius/ 3));
+			for($x=$x0; $x<=$x1; ++$x){//Loop through x value
+				if(0==(($x-$x0) % $step)){
+					if($steep){
+						self::_drawCircularGradient(&$im, $y, $x, $spot_radius, $dimming);
+					}else{ 
+						self::_drawCircularGradient(&$im, $x, $y, $spot_radius, $dimming);
+					}
+				}
+				$error -= $deltay;
+				if($error<0){
+						$y = $y + $ystep;
+						$error = $error + $deltax;
+				}
+			}		
+		}//_drawBilinearGradient
+		
+		private static function _createGradient($im, $mode, $gradient_name){
+			//create the gradient from an image
+			if(FALSE === ($grad_im = imagecreatefrompng('gradient/'.$gradient_name.'.png'))){
+				return FALSE;
+			}
+			$width_g = imagesx($grad_im);
+			$height_g = imagesy($grad_im);
+			//Get colors along the longest dimension
+			//Max density is for lower channel value
+			for($y=$height_g-1; $y >= 0 ; --$y){
+					$rgb = imagecolorat($grad_im, 1, $y);
+					//Linear function
+					$alpha = Min(127, Max(0, floor(127 - $y/2)));
+					if(self::$WITH_ALPHA == $mode){
+						$grad_rgba[] = imagecolorallocatealpha($im, ($rgb >> 16) & 0xFF, ($rgb >> 8) & 0xFF, $rgb & 0xFF, $alpha);
+					}else{
+						$grad_rgba[] = imagecolorallocate($im, ($rgb >> 16) & 0xFF, ($rgb >> 8) & 0xFF, $rgb & 0xFF);
+					}
+			}
+			imagedestroy($grad_im);
+			unset($grad_im);
+			return($grad_rgba);
+		}//_createGradient
+	}//Heatmap
+?>