Add feature to merge lidar tiles into single dem entry

common.h

@@ -22,6 +22,8 @@
 #define	KM_PER_MILE	1.609344
 #define FOUR_THIRDS	1.3333333333333
 
+#define MAX(x,y)((x)>(y)?(x):(y))
+
 struct dem {
 	float min_north;
 	float max_north;

inputs.cc

@@ -322,58 +322,123 @@ int loadLIDAR(char *filenames)
 
 	}
 
-	/* Iterate through all tiles to find the largest x/y dimension
+	/* Iterate through all of the tiles to find the smallest resolution. We will
-	to use as the global IPPD value. We do this here so that we are
+	 * need to rescale every tile from here on out to this value */
-	sure to allocate enough memory for the remainder of the calculations */
+	int smallest_res = 0;
-	int dimension_max = -1;
+	int pix_per_deg = 0;
 	for (size_t i = 0; i < fc; i++) {
-		if( tiles[i].width > tiles[i].height && tiles[i].width > dimension_max )
+		if ( smallest_res == 0 || tiles[i].resolution < smallest_res ){
-			dimension_max = tiles[i].width;
+			smallest_res = tiles[i].resolution;
-		else if( tiles[i].height > dimension_max )
+			pix_per_deg = MAX(tiles[i].width,tiles[i].height) / MAX(tiles[i].max_north - tiles[i].min_north, tiles[i].max_west - tiles[i].min_west);
-			dimension_max = tiles[i].height;
+		}
 	}
 
-	MAXPAGES = fc;
+	/* Now we need to rescale all tiles the the lowest resolution. ie if we have
-	IPPD = dimension_max;
+	 * one 1m lidar and one 2m lidar, resize the 2m to fake 1m */
+	for (size_t i = 0; i< fc; i++) {
+		float rescale = tiles[i].resolution / smallest_res;
+		if (tiles[i].resolution != smallest_res)
+			tile_rescale(&tiles[i],rescale);
+	}
+
+	/* Now we work out the size of the giant lidar tile. */
+	double total_width = max_west - min_west;
+	double total_height = max_north - min_north;
+	/* This is how we should _theoretically_ work this out, but due to
+	 * the nature of floating point arithmetic and rounding errors, we need to
+	 * crunch the numbers the hard way */
+	// size_t new_width = total_width * pix_per_deg;
+	// size_t new_height = total_height * pix_per_deg;
+
+	size_t new_height = 0;
+	size_t new_width = 0;
+	for ( size_t i = 0; i < fc; i++ ) {
+		double north_offset = max_north - tiles[i].max_north;
+		double west_offset = max_west - tiles[i].max_west;
+		size_t north_pixel_offset = north_offset * pix_per_deg;
+		size_t west_pixel_offset = west_offset * pix_per_deg;
+
+		if ( west_pixel_offset + tiles[i].width > new_width )
+			new_width = west_pixel_offset + tiles[i].width;
+		if ( north_pixel_offset + tiles[i].height > new_height )
+			new_height = north_pixel_offset + tiles[i].height;
+	}
+	size_t new_tile_alloc = new_width * new_height;
+
+	int * new_tile = (int*) calloc( new_tile_alloc, sizeof(int) );
+	if (debug)
+		fprintf(stderr,"Lidar tile dimensions w:%lf(%zu) h:%lf(%zu)\n", total_width, new_width, total_height, new_height);
+
+	/* ...If we wanted a value other than sea level here, we would need to initialize the array... */
+
+	/* Fill out the array one tile at a time */
+	for (size_t i = 0; i< fc; i++) {
+		double north_offset = max_north - tiles[i].max_north;
+		double west_offset = max_west - tiles[i].max_west;
+		size_t north_pixel_offset = north_offset * pix_per_deg;
+		size_t west_pixel_offset = west_offset * pix_per_deg;
+
+		if (debug) {
+			fprintf(stderr,"mn: %lf mw:%lf globals: %lf %lf\n", tiles[i].max_north, tiles[i].max_west, max_north, max_west);
+			fprintf(stderr,"Offset n:%zu w:%zu\n", north_pixel_offset, west_pixel_offset);
+			fprintf(stderr,"Height: %d\n", tiles[i].height);
+		}
+
+		/* Copy it row-by-row from the tile */
+		for (size_t h = 0; h < tiles[i].height; h++) {
+			register int *dest_addr = &new_tile[ (north_pixel_offset+h)*new_width + west_pixel_offset];
+			register int *src_addr = &tiles[i].data[h*tiles[i].width];
+			// Check if we might overflow
+			if ( dest_addr + tiles[i].width > new_tile + new_tile_alloc ){
+				if (debug)
+					fprintf(stderr, "Overflow %zu\n",i);
+				continue;
+			}
+			//fprintf(stderr,"dest:%p src:%p\n", dest_addr, src_addr);
+			memcpy( dest_addr, src_addr, tiles[i].width * sizeof(int) );
+		}
+	}
+
+	MAXPAGES = 1;
+	IPPD = MAX(new_width,new_height);
 	if(debug){
 		fprintf(stderr,"Setting IPPD to %d\n",IPPD);
 		fflush(stderr);
 	}
 	// add fudge as reprojected tiles sometimes vary by a pixel or ten
-	IPPD += 50;
+	// IPPD += 50;
 	ARRAYSIZE = (MAXPAGES * IPPD) + 50;
 	do_allocs();
 	// reset the IPPD after allocations
-	IPPD -= 50;
+	// IPPD -= 50;
 
 	/* Load the data into the global dem array */
-	for (size_t indx = 0; indx < fc; indx++) {
+	dem[0].max_north = max_north;
-		dem[indx].max_north = tiles[indx].yur;
+	dem[0].min_west = min_west;
-		dem[indx].min_west = tiles[indx].xur;
+	dem[0].min_north = min_north;
-		dem[indx].min_north = tiles[indx].yll;
+	dem[0].max_west = max_west;
-		dem[indx].max_west = tiles[indx].xll;
+	dem[0].max_el = max_elevation;
-		dem[indx].max_el = tiles[indx].max_el;
+	dem[0].min_el = min_elevation;
-		dem[indx].min_el = tiles[indx].min_el;
 
 	/* 
 	 * Copy the lidar tile data into the dem array. The dem array is rotated
 	 * 90 degrees (christ knows why...)
 	 */
-		int y = tiles[indx].height-1;
+	int y = new_height-1;
-		for (size_t h = 0; h < tiles[indx].height; h++, y--) {
+	for (size_t h = 0; h < new_height; h++, y--) {
-			int x = tiles[indx].width-1;
+		int x = new_width-1;
-			for (size_t w = 0; w < tiles[indx].width; w++, x--) {
+		for (size_t w = 0; w < new_width; w++, x--) {
-				dem[indx].data[y][x] = tiles[indx].data[h*tiles[indx].width + w];
+			dem[0].data[y][x] = new_tile[h*new_width + w];
-				dem[indx].signal[y][x] = 0;
+			dem[0].signal[y][x] = 0;
-				dem[indx].mask[y][x] = 0;
+			dem[0].mask[y][x] = 0;
 		}
 	}
 
-	}
-
 	ippd=IPPD;
-	height = (unsigned)((max_north-min_north) / smCellsize);
+	// height = (unsigned)((max_north-min_north) / smCellsize);
-	width = (unsigned)((max_west-min_west) / smCellsize);
+	// width = (unsigned)((max_west-min_west) / smCellsize);
+	height = (unsigned)((max_north-min_north) * pix_per_deg);
+	width = (unsigned)((max_west-min_west) * pix_per_deg);
 
 	if (debug)
 		fprintf(stderr, "LIDAR LOADED %d x %d\n", width, height);

tiles.cc

@@ -2,11 +2,27 @@
 #include <stdio.h>
 #include <errno.h>
 #include <string.h>
+#include <math.h>
 #include "tiles.hh"
 #include "common.h"
 
 #define MAX_LINE 25000
 
+/* Computes the distance between two long/lat points */
+double haversine_formulaz(double th1, double ph1, double th2, double ph2)
+{
+	#define TO_RAD (3.1415926536 / 180)
+	int R = 6371;
+	double dx, dy, dz;
+	ph1 -= ph2;
+	ph1 *= TO_RAD, th1 *= TO_RAD, th2 *= TO_RAD;
+ 
+	dz = sin(th1) - sin(th2);
+	dx = cos(ph1) * cos(th1) - cos(th2);
+	dy = sin(ph1) * cos(th1);
+	return asin(sqrt(dx * dx + dy * dy + dz * dz) / 2) * 2 * R;
+}
+
 int tile_load_lidar(tile_t *tile, char *filename){
 	FILE *fd;
 	char line[MAX_LINE];
@@ -100,6 +116,9 @@ int tile_load_lidar(tile_t *tile, char *filename){
 		}//if
 	}
 
+	double current_res_km = haversine_formulaz(tile->max_north, tile->max_west, tile->max_north, tile->min_west);
+	tile->resolution = (int) ceil((current_res_km/MAX(tile->width,tile->height))*1000);
+
 	if (debug)
 		fprintf(stderr,"Pixels loaded: %zu/%d\n",loaded,tile->width*tile->height);
 
@@ -108,3 +127,78 @@ int tile_load_lidar(tile_t *tile, char *filename){
 
 	return 0;
 }
+
+/*
+ *	A positive scale will _increase_ the size of the data
+ */
+int tile_rescale(tile_t *tile, float scale){
+	int *new_data;
+	size_t skip_count = 1;
+	size_t copy_count = 1;
+
+	size_t new_height = tile->height * scale;
+	size_t new_width = tile->width * scale;
+
+	/* Allocate the array for the lidar data */
+	if ( (new_data = (int*) calloc(new_height * new_width, sizeof(int))) == NULL ) {
+		return ENOMEM;
+	}
+
+	tile->max_el = -32768;
+	tile->min_el = 32768;
+
+	/* Making the tile data smaller */
+	if (scale < 0) {
+		skip_count = 1 / scale;
+	} else {
+		copy_count = (size_t) scale;
+	}
+
+	fprintf(stderr,"Skip: %zu Copy: %zu\n", skip_count, copy_count);
+
+	
+
+	/* Nearest neighbour normalization. For each subsample of the original, simply
+	 * assign the value in the top left to the new pixel */
+	if (scale < 0){
+		for (size_t x = 0, i = 0; i < new_width; x += skip_count, i++) {
+			for (size_t y = 0, j = 0; j < new_height; y += skip_count, j++){
+				new_data[i*new_width+j] = tile->data[x*tile->width+y];
+				/* Update local min / max values */
+				if (tile->data[x * tile->width + y] > tile->max_el)
+					tile->max_el = tile->data[x * tile->width + y];
+				if (tile->data[x * tile->width + y] < tile->min_el)
+					tile->min_el = tile->data[x * tile->width + y];
+			}
+		}
+	}else{
+		for (size_t x = 0; x < tile->width; x++) {
+			for (size_t y = 0; y < tile->height; y++){
+				/* These are for scaling up the data */
+				for (size_t copy_x = 0; copy_x < copy_count; copy_x++) {
+					for (size_t copy_y = 0; copy_y < copy_count; copy_y++) {
+						size_t new_x = (x * skip_count) + copy_x;
+						size_t new_y = (y * skip_count) + copy_y;
+						new_data[ new_x * new_width + new_y ] = tile->data[x * tile->width + y];
+						// new_data[(x + copy_x) * new_width + (y + copy_y)] = tile->data[x * tile->width + y];
+					}
+				}
+				/* Update local min / max values */
+				if (tile->data[x * tile->width + y] > tile->max_el)
+					tile->max_el = tile->data[x * tile->width + y];
+				if (tile->data[x * tile->width + y] < tile->min_el)
+					tile->min_el = tile->data[x * tile->width + y];
+			}
+		}
+	}
+
+	/* Update the date in the tile */
+	free(tile->data);
+	tile->data = new_data;
+
+	/* Update the height and width values */
+	tile->height = new_height;
+	tile->width = new_width;
+
+	return 0;
+}

tiles.hh

@@ -33,8 +33,10 @@ typedef struct _tile_t{
 	int 	max_el;
 	int		min_el;
 	int		*data;
+	int		resolution;
 } tile_t, *ptile_t;
 
 int tile_load_lidar(tile_t*, char *);
+int tile_rescale(tile_t *, float);
 
 #endif