1 /***************************************************************************
2 *
3 * Project: libLAS -- C/C++ read/write library for LAS LIDAR data
4 * Purpose: LAS translation to MonetDB binary format with optional configuration
5 * Author: Romulo Goncalves r.goncalves@esciencecenter.nl
6 Oscar Martinez Rubi o.rubi@esciencecenter.nl
7 ***************************************************************************
8 * This tool has been developed by the Netherlands eScience Center
9 * (https://www.esciencecenter.nl/)
10 *
11 * Copyright (c) 2016, Romulo Goncalves r.goncalves@esciencecenter.nl
12 *
13 * See LICENSE.txt in this source distribution for more information.
14 **************************************************************************/
15 #include <sys/stat.h>
16 #include <stdio.h>
17 #include <stdlib.h>
18 #include <string.h>
19 #include <pthread.h>
20 #include <assert.h>
21
22 #include "liblas.h"
23 #include "lascommon.h"
24
25 #include <limits.h>
26 #include <inttypes.h>
27 #if defined(__linux__) || defined(__CYGWIN__) || defined(__FreeBSD_kernel__) || defined(__GNU__)
28
29 #include <endian.h>
30 #include <unistd.h>
31 #elif defined(__APPLE__)
32
33 #include <unistd.h>
34 # include <libkern/OSByteOrder.h>
35
36 # define htobe16(x) OSSwapHostToBigInt16(x)
37 # define htole16(x) OSSwapHostToLittleInt16(x)
38 # define be16toh(x) OSSwapBigToHostInt16(x)
39 # define le16toh(x) OSSwapLittleToHostInt16(x)
40
41 # define htobe32(x) OSSwapHostToBigInt32(x)
42 # define htole32(x) OSSwapHostToLittleInt32(x)
43 # define be32toh(x) OSSwapBigToHostInt32(x)
44 # define le32toh(x) OSSwapLittleToHostInt32(x)
45
46 # define htobe64(x) OSSwapHostToBigInt64(x)
47 # define htole64(x) OSSwapHostToLittleInt64(x)
48 # define be64toh(x) OSSwapBigToHostInt64(x)
49 # define le64toh(x) OSSwapLittleToHostInt64(x)
50
51 # define __BYTE_ORDER BYTE_ORDER
52 # define __BIG_ENDIAN BIG_ENDIAN
53 # define __LITTLE_ENDIAN LITTLE_ENDIAN
54 # define __PDP_ENDIAN PDP_ENDIAN
55
56 #elif defined(__OpenBSD__) || defined(__FreeBSD__) || defined(__DragonFly__)
57
58 #include <unistd.h>
59 # include <sys/endian.h>
60
61 #elif defined(__NetBSD__)
62
63 # include <sys/endian.h>
64
65 # define be16toh(x) betoh16(x)
66 # define le16toh(x) letoh16(x)
67
68 # define be32toh(x) betoh32(x)
69 # define le32toh(x) letoh32(x)
70
71 # define be64toh(x) betoh64(x)
72 # define le64toh(x) letoh64(x)
73
74 #elif defined(__WINDOWS__)
75
76 # include <winsock2.h>
77 # include <sys/param.h>
78
79 # if BYTE_ORDER == LITTLE_ENDIAN
80
81 # define htobe16(x) htons(x)
82 # define htole16(x) (x)
83 # define be16toh(x) ntohs(x)
84 # define le16toh(x) (x)
85
86 # define htobe32(x) htonl(x)
87 # define htole32(x) (x)
88 # define be32toh(x) ntohl(x)
89 # define le32toh(x) (x)
90
91 # define htobe64(x) htonll(x)
92 # define htole64(x) (x)
93 # define be64toh(x) ntohll(x)
94 # define le64toh(x) (x)
95
96 # elif BYTE_ORDER == BIG_ENDIAN
97
98 /* that would be xbox 360 */
99 # define htobe16(x) (x)
100 # define htole16(x) __builtin_bswap16(x)
101 # define be16toh(x) (x)
102 # define le16toh(x) __builtin_bswap16(x)
103
104 # define htobe32(x) (x)
105 # define htole32(x) __builtin_bswap32(x)
106 # define be32toh(x) (x)
107 # define le32toh(x) __builtin_bswap32(x)
108
109 # define htobe64(x) (x)
110 # define htole64(x) __builtin_bswap64(x)
111 # define be64toh(x) (x)
112 # define le64toh(x) __builtin_bswap64(x)
113
114 # else
115
116 # error byte order not supported
117
118 # endif
119
120 # define __BYTE_ORDER BYTE_ORDER
121 # define __BIG_ENDIAN BIG_ENDIAN
122 # define __LITTLE_ENDIAN LITTLE_ENDIAN
123 # define __PDP_ENDIAN PDP_ENDIAN
124
125 #else
126
127 # error platform not supported
128
129 #endif
130
131 #include <math.h>
132
133 #ifndef _BSD_SOURCE
134 #define _BSD_SOURCE
135 #endif
136
137 #define boolean short
138 #define int64_t long long int
139
140 #define NUM_OF_ENTRIES 21
141 #define DEFAULT_NUM_READ_THREADS 1
142 #define DEFAULT_NUM_INPUT_FILES 2000
143 #define TOLERANCE 0.0000001
144 #define MAX_INT_31 2147483648.0
145
146 #define set_lock(lock, s) \
147 { MT_lock_set(&lock, s);}
148 #define unset_lock(node, lock, s) \
149 { MT_lock_unset(&lock, s);}
150
151
152
153
154 void print_header(FILE *file, LASHeaderH header, const char* file_name);
155
usage()156 void usage()
157 {
158 fprintf(stderr,"----------------------------------------------------------\n");
159 fprintf(stderr," las2col (version %s) usage:\n", LAS_GetVersion());
160 fprintf(stderr,"----------------------------------------------------------\n");
161 fprintf(stderr,"\n");
162
163 fprintf(stderr,"Convert a las/laz file into columnar format (binary) of MonetDB, outputs for each entry a file <output_prefix>_col_<entry_name>.dat:\n");
164 fprintf(stderr," las2col -i <input_file> -o <output_prefix>\n");
165 fprintf(stderr,"\n");
166
167 fprintf(stderr,"Convert a list of las/laz files (still outputs for each entry a file <output_prefix>_col_<entry_name>.dat):\n");
168 fprintf(stderr," las2col -i <las_file_1> -i <las_file_2> -o <output_prefix>\n");
169 fprintf(stderr,"Alternatively:\n");
170 fprintf(stderr," las2col -f <file_with_the_list_las/laz_files> -o <output_prefix>\n");
171 fprintf(stderr,"\n");
172
173 fprintf(stderr,"Convert a list of las/laz files using <num_read_threads> threads (default is 1):\n");
174 fprintf(stderr," las2col -f <file_with_the_list_las/laz_files> -o <output_prefix> --num_read_threads <number_of_threads>\n");
175 fprintf(stderr,"\n\n");
176
177 fprintf(stderr,"----------------------------------------------------------\n");
178 fprintf(stderr," The '--parse txyz' flag specifies which entries of the LAS/LAZ\n");
179 fprintf(stderr," will be extracted (default is --parse xyz). For example, 'txyzia'\n");
180 fprintf(stderr," means that six columnar (binary) MonetDB files will be generated,\n");
181 fprintf(stderr," the first one containing all gpstime values, \n");
182 fprintf(stderr," the next three containing values for x, y, and\n");
183 fprintf(stderr," z coordinates, the next one with intensity values\n");
184 fprintf(stderr," and the last one with scan angle values.\n");
185 fprintf(stderr," The supported entries are:\n");
186 fprintf(stderr," t - gpstime as double\n");
187 fprintf(stderr," x - x coordinate as double\n");
188 fprintf(stderr," y - y coordinate as double\n");
189 fprintf(stderr," z - z coordinate as double\n");
190 fprintf(stderr," X - x coordinate as decimal(<num_digits_unscaled_max_x>,<num_digits_scale_x>)\n");
191 fprintf(stderr," Y - y coordinate as decimal(<num_digits_unscaled_max_y>,<num_digits_scale_y>)\n");
192 fprintf(stderr," Z - z coordinate as decimal(<num_digits_unscaled_max_z>,<num_digits_scale_z>)\n");
193 fprintf(stderr," a - scan angle as tinyint\n");
194 fprintf(stderr," i - intensity as smallint\n");
195 fprintf(stderr," n - number of returns for given pulse as smallint\n");
196 fprintf(stderr," r - number of this return as smallint\n");
197 fprintf(stderr," c - classification number as tinyint\n");
198 fprintf(stderr," u - user data as tinyint\n");
199 fprintf(stderr," p - point source ID as smallint\n");
200 fprintf(stderr," e - edge of flight line as smallint\n");
201 fprintf(stderr," d - direction of scan flag as smallint\n");
202 fprintf(stderr," R - red channel of RGB color as smallint\n");
203 fprintf(stderr," G - green channel of RGB color as smallint\n");
204 fprintf(stderr," B - blue channel of RGB color as smallint\n");
205 fprintf(stderr," M - vertex index number as integer\n");
206 fprintf(stderr," k - Morton 2D code using X and Y (unscaled and no offset) as bigint\n\n");
207
208 fprintf(stderr," The '--moffset 8600000,40000000' flag specifies a global offset in X and Y \n");
209 fprintf(stderr," to be used when computing the Morton 2D code. Values must be unscaled \n\n");
210
211 fprintf(stderr," The '--check 0.01,0.01' flag checks suitability to compute Morton 2D codes \n");
212 fprintf(stderr," It checks specified scale matches the one in input file. \n");
213 fprintf(stderr," If moffset is provided it also checks that obtained Morton 2D codes \n");
214 fprintf(stderr," will be consistent, i.e. global X,Y within [0,2^31] \n\n");
215
216 fprintf(stderr,"----------------------------------------------------------\n");
217
218 fprintf(stderr," After generating the columnar files, import them in MonetDB. Example: \n");
219 fprintf(stderr," mclient <db_name> -s \"COPY BINARY INTO flat FROM ('<full_parent_path>/out_col_x.dat','<full_parent_path>/out_col_y.dat','<full_parent_path>/out_col_z.dat')\"\n");
220 fprintf(stderr," Note that full paths of the columnar files MUST be used. Also note that a table called flat has to be created in a MonetDB DB beforehand. The table must have \n");
221 fprintf(stderr," the columns in the same order as specified by the --parse option, and the column types must be the ones specified above. Example: \n");
222 fprintf(stderr," mclient <db_name> -s \"create table flat (x double, y double, z double)\"\n");
223 fprintf(stderr," Note that for decimal entries (XYZ) the column definition at table-creation time must be decimal(<num_digits_unscaled_max>,<num_digits_scale>)\n");
224 fprintf(stderr," For example, if the maximum X value of a file (or a list of files) is 638982.55, then the X definition when creating the table is decimal(8,2). Example:\n");
225 fprintf(stderr," mclient <db_name> -s \"create table flat (x decimal(8,2), y decimal(8,2), z decimal(8,2))\"\n");
226 }
227
228 /*Global structures*/
229 #define MT_Lock pthread_mutex_t
230 #define MT_set_lock(p) pthread_mutex_lock(p)
231 #define MT_unset_lock(p) pthread_mutex_unlock(p)
232 #define MT_lock_init(p) pthread_mutex_init(p,NULL)
233 #define MT_lock_destroy(p) pthread_mutex_destroy(p)
234
235 #define MT_Cond pthread_cond_t
236 #define MT_cond_wait(p,t) pthread_cond_wait(p,t)
237 #define MT_cond_init(p) pthread_cond_init(p,NULL)
238 #define MT_cond_destroy(p) pthread_cond_destroy(p)
239
240 typedef void (*f_ptr)( void );
241
242 MT_Lock dataLock;
243 MT_Cond mainCond, writeTCond, readCond;
244 int entries[NUM_OF_ENTRIES];
245 double (*entriesFuncD[NUM_OF_ENTRIES])();
246 int (*entriesFuncI[NUM_OF_ENTRIES])();
247 short (*entriesFuncS[NUM_OF_ENTRIES])();
248 char (*entriesFuncC[NUM_OF_ENTRIES])();
249 int entriesType[NUM_OF_ENTRIES];
250 char **files_name_in = NULL;
251 int files_in_index = 0 ;
252 int skip_invalid = FALSE;
253 int verbose = TRUE;
254 struct writeT **data = NULL;
255 struct writeT *dataWriteT = NULL;
256 int stop;
257
258 typedef enum {
259 ENTRY_x,
260 ENTRY_y,
261 ENTRY_z,
262 ENTRY_X,
263 ENTRY_Y,
264 ENTRY_Z,
265 ENTRY_t,
266 ENTRY_i,
267 ENTRY_a,
268 ENTRY_r,
269 ENTRY_c,
270 ENTRY_u,
271 ENTRY_n,
272 ENTRY_R,
273 ENTRY_G,
274 ENTRY_B,
275 ENTRY_M,
276 ENTRY_p,
277 ENTRY_e,
278 ENTRY_d,
279 ENTRY_k
280 } ENTRIES;
281
282 struct writeThreadArgs {
283 int id;
284 FILE *out;
285 };
286
287 struct writeT {
288 long num_points;
289 char* values;
290 int type;
291 };
292
293 struct readThreadArgs {
294 int id;
295 int num_read_threads;
296 int num_of_entries;
297 int check;
298 int64_t global_offset_x;
299 int64_t global_offset_y;
300 double scale_x;
301 double scale_y;
302 };
303
writeFile(void * arg)304 void* writeFile(void *arg) {
305 int i = 0;
306 struct writeThreadArgs *wTA = (struct writeThreadArgs*) arg;
307
308 /*Obtain lock over data to get the pointer*/
309 while (stop == 0) {
310 MT_set_lock(&dataLock);
311 while ((stop == 0) && (dataWriteT == NULL || (dataWriteT && dataWriteT[wTA->id].values == NULL))) {
312 /*Sleep and wait for data to be read*/
313 MT_cond_wait(&writeTCond,&dataLock);
314 }
315 //Release the lock
316 MT_unset_lock(&dataLock);
317
318 if (stop) {
319 return NULL;
320 }
321
322 fwrite(dataWriteT[wTA->id].values, dataWriteT[wTA->id].type, dataWriteT[wTA->id].num_points, wTA->out);
323 //for (i = 0; i < 1000; i++)
324 // printf("%d\n", dataWriteT[wTA->id].values[i]);
325 MT_set_lock(&dataLock);
326 free(dataWriteT[wTA->id].values);
327 dataWriteT[wTA->id].values = NULL;
328 MT_unset_lock(&dataLock);
329 fflush(wTA->out);
330
331 /*Wake up the main*/
332 pthread_cond_broadcast(&mainCond);
333 }
334 return NULL;
335 }
336
readFile(void * arg)337 void* readFile(void *arg) {
338 struct readThreadArgs *rTA = (struct readThreadArgs*) arg;
339 LASReaderH reader = NULL;
340 LASHeaderH header = NULL;
341 LASPointH p = NULL;
342 unsigned int index = 0;
343 int read_index = 0;
344 char *file_name_in = NULL;
345 int i, j;
346
347 while(1) {
348 file_name_in = NULL;
349 /*Get next file to read*/
350 MT_set_lock(&dataLock);
351 file_name_in = files_name_in[files_in_index];
352 if (file_name_in == NULL) {
353 MT_unset_lock(&dataLock);
354 return NULL;
355 }
356 read_index = (files_in_index % rTA->num_read_threads);
357 files_in_index++;
358
359 struct writeT *dataWriteTT = (struct writeT*) malloc(sizeof(struct writeT)*rTA->num_of_entries);
360 /*Lets read the data*/
361 reader = LASReader_Create(file_name_in);
362 if (!reader) {
363 LASError_Print("Unable to read file");
364 MT_unset_lock(&dataLock);
365 exit(1);
366 }
367 MT_unset_lock(&dataLock);
368
369 header = LASReader_GetHeader(reader);
370 if (!header) {
371 LASError_Print("Unable to fetch header for file");
372 exit(1);
373 }
374
375 if (verbose)
376 {
377 print_header(stderr, header, file_name_in);
378 }
379
380 /*Allocate arrays for the columns*/
381 long num_points = LASHeader_GetPointRecordsCount(header);
382 for (i = 0; i < rTA->num_of_entries; i++) {
383 dataWriteTT[i].num_points = num_points;
384 dataWriteTT[i].values = malloc(entriesType[i]*num_points);
385 dataWriteTT[i].type = entriesType[i];
386 }
387
388 /*Changes for Oscar's new Morton code function*/
389 //unsigned int factorX = (unsigned int) (LASHeader_GetOffsetX(header) / LASHeader_GetScaleX(header));
390 //unsigned int factorY = (unsigned int) (LASHeader_GetOffsetY(header) / LASHeader_GetScaleY(header));
391
392 /*Compute factors to add to X and Y and cehck sanity of generated codes*/
393 double file_scale_x = LASHeader_GetScaleX(header);
394 double file_scale_y = LASHeader_GetScaleY(header);
395 double file_scale_z = LASHeader_GetScaleZ(header);
396 //printf("The scales are x:%lf y:%lf z:%lf\n", file_scale_x, file_scale_y, file_scale_z);
397
398 /* scaled offsets to add for the morton encoding */
399 int64_t factorX = ((int64_t) (LASHeader_GetOffsetX(header) / file_scale_x)) - rTA->global_offset_x;
400 int64_t factorY = ((int64_t) (LASHeader_GetOffsetY(header) / file_scale_y)) - rTA->global_offset_y;
401
402 if (rTA->check)
403 {
404 // Check specified scales are like in the LAS file
405 if (fabs(rTA->scale_x - file_scale_x) > TOLERANCE){
406 fprintf(stderr, "ERROR: x scale in input file (%lf) does not match specified x scale (%lf)\n",file_scale_x, rTA->scale_x);
407 exit(1);
408 }
409 if (fabs(rTA->scale_y - file_scale_y) > TOLERANCE){
410 fprintf(stderr, "ERROR: y scale in input file (%lf) does not match specified y scale (%lf)\n",file_scale_y, rTA->scale_y);
411 exit(1);
412 }
413 /* Check that the extent of the file (taking into account the global offset)
414 * is within 0,2^31 */
415 double check_min_x = 1.0 + LASHeader_GetMinX(header) - (((double) rTA->global_offset_x) * rTA->scale_x);
416 if (check_min_x < TOLERANCE) {
417 fprintf(stderr, "ERROR: Specied X global offset is too large. (MinX - (GlobalX*ScaleX)) < 0\n");
418 exit(1);
419 }
420 double check_min_y = 1.0 + LASHeader_GetMinY(header) - (((double) rTA->global_offset_y) * rTA->scale_y);
421 if (check_min_y < TOLERANCE) {
422 fprintf(stderr, "ERROR: Specied Y global offset is too large. (MinY - (GlobalY*ScaleY)) < 0\n");
423 exit(1);
424 }
425 double check_max_x = LASHeader_GetMaxX(header) - (((double) rTA->global_offset_x) * rTA->scale_x);
426 if (check_max_x > (MAX_INT_31 * rTA->scale_x)) {
427 fprintf(stderr, "ERROR: Specied X global offset is too small. (MaxX - (GlobalX*ScaleX)) > (2^31)*ScaleX\n");
428 exit(1);
429 }
430 double check_max_y = LASHeader_GetMaxY(header) - (((double) rTA->global_offset_y) * rTA->scale_y);
431 if (check_max_y > (MAX_INT_31 * rTA->scale_y)) {
432 fprintf(stderr, "ERROR: Specied Y global offset is too small. (MaxY - (GlobalY*ScaleY)) > (2^31)*ScaleY\n");
433 exit(1);
434 }
435 }
436
437 p = LASReader_GetNextPoint(reader);
438 index = 0;
439 while (p)
440 {
441 if (skip_invalid && !LASPoint_IsValid(p)) {
442 if (verbose) {
443 LASError_Print("Skipping writing invalid point...");
444 }
445 p = LASReader_GetNextPoint(reader);
446 index -=1;
447 continue;
448 }
449
450 LASColorH color = NULL;
451 for (j = 0; j < rTA->num_of_entries; j++) {
452 uint64_t res;
453 switch (entries[j]) {
454 case ENTRY_x:
455 case ENTRY_y:
456 case ENTRY_z:
457 case ENTRY_t:
458 ((double*) dataWriteTT[j].values)[index] = entriesFuncD[j](p);
459 //printf(" Point is:%lf\n", ((double*) dataWriteTT[j].values)[index]);
460 break;
461 case ENTRY_X:
462 ((int*) dataWriteTT[j].values)[index] = entriesFuncD[j](p) / file_scale_x;
463 break;
464 case ENTRY_Y:
465 ((int*) dataWriteTT[j].values)[index] = entriesFuncD[j](p) / file_scale_y;
466 break;
467 case ENTRY_Z:
468 ((int*) dataWriteTT[j].values)[index] = entriesFuncD[j](p) / file_scale_z;
469 break;
470 case ENTRY_i:
471 case ENTRY_r:
472 case ENTRY_n:
473 case ENTRY_p:
474 case ENTRY_e:
475 case ENTRY_d:
476 ((short*) dataWriteTT[j].values)[index] = entriesFuncS[j](p);
477 break;
478 case ENTRY_a:
479 case ENTRY_c:
480 case ENTRY_u:
481 ((char*) dataWriteTT[j].values)[index] = entriesFuncC[j](p);
482 break;
483 case ENTRY_k:
484 entriesFuncD[j](&res, p, factorX, factorY);
485 ((int64_t*)dataWriteTT[j].values)[index] = res;
486 break;
487 case ENTRY_R:
488 case ENTRY_G:
489 case ENTRY_B:
490 color = (color == NULL) ? LASPoint_GetColor(p) : color;
491 ((unsigned short*) dataWriteTT[j].values)[index] = entriesFuncS[j](color);;
492 break;
493 case ENTRY_M:
494 ((unsigned int*)dataWriteTT[j].values)[index] = index;
495 break;
496 default:
497 LASError_Print("las2col:readFile: Invalid Entry.");
498 }
499 }
500 if (color != NULL)
501 LASColor_Destroy(color);
502
503 p = LASReader_GetNextPoint(reader);
504 index +=1;
505 }
506 if (verbose)
507 printf("Num of points:%d %ld for file:%s \n", index, num_points, file_name_in);
508
509 /*Give the data to the writer threads*/
510 MT_set_lock(&dataLock);
511 LASHeader_Destroy(header);
512 header = NULL;
513 LASReader_Destroy(reader);
514 reader = NULL;
515
516 /*TODO: make sure you are not overtaking other reading threads*/
517 while (data[read_index] != NULL) {
518 MT_cond_wait(&readCond, &dataLock);
519 }
520 data[read_index] = dataWriteTT;
521 /*Wake up the main*/
522 pthread_cond_broadcast(&mainCond);
523 MT_unset_lock(&dataLock);
524
525 }
526 return NULL;
527 }
528
doesFileExist(const char * filename)529 int doesFileExist(const char *filename) {
530 struct stat st;
531 int result = stat(filename, &st);
532 return result == 0;
533 }
534
EncodeMorton2D_1(unsigned int rawx,unsigned int rawy)535 int64_t EncodeMorton2D_1(unsigned int rawx, unsigned int rawy){
536 int64_t answer = 0;
537 int64_t i;
538 for (i = 0; i < (sizeof(int64_t)* CHAR_BIT)/2; ++i) {
539 answer |= ((rawy & ((int64_t)1 << i)) << i) | ((rawx & ((int64_t)1 << i)) << (i + 1));
540 }
541 return answer;
542 }
543
Expand1(uint32_t a)544 uint64_t Expand1(uint32_t a)
545 {
546 uint64_t b = a & 0x7fffffff; // b = ---- ---- ---- ---- ---- ---- ---- ---- 0edc ba98 7654 3210 fedc ba98 7654 3210
547 b = (b ^ (b << 16)) & 0x0000ffff0000ffff; // b = ---- ---- ---- ---- 0edc ba98 7654 3210 ---- ---- ---- ---- fedc ba98 7654 3210
548 b = (b ^ (b << 8)) & 0x00ff00ff00ff00ff; // b = ---- ---- 0edc ba98 ---- ---- 7654 3210 ---- ---- fedc ba98 ---- ---- 7654 3210
549 b = (b ^ (b << 4)) & 0x0f0f0f0f0f0f0f0f; // b = ---- 0edc ---- ba98 ---- 7654 ---- 3210 ---- fedc ---- ba98 ---- 7654 ---- 3210
550 b = (b ^ (b << 2)) & 0x3333333333333333; // b = --0e --dc --ba --98 --76 --54 --32 --10 --fe --dc --ba --98 --76 --54 --32 --10
551 b = (b ^ (b << 1)) & 0x5555555555555555; // b = -0-e -d-c -b-a -9-8 -7-6 -5-4 -3-2 -1-0 -f-e -d-c -b-a -9-8 -7-6 -5-4 -3-2 -1-0
552 return b;
553 }
554
morton2D_encode(int64_t * answer,LASPointH p,unsigned int factorX,unsigned int factorY)555 int64_t morton2D_encode(int64_t *answer, LASPointH p, unsigned int factorX, unsigned int factorY){
556 unsigned int rawx = ((unsigned int) LASPoint_GetRawX(p)) + factorX;
557 unsigned int rawy = ((unsigned int) LASPoint_GetRawY(p)) + factorY;
558 *answer = EncodeMorton2D_1(rawx, rawy);
559 return *answer;
560 }
561
562 /*Changes for Oscar's new Morton code function*/
morton2D_encodeOscar(uint64_t * answer,LASPointH p,unsigned int factorX,unsigned int factorY)563 uint64_t morton2D_encodeOscar(uint64_t *answer, LASPointH p, unsigned int factorX, unsigned int factorY){
564 uint32_t x = (uint32_t) (((int64_t) LASPoint_GetRawX(p)) + factorX);
565 uint32_t y = (uint32_t) (((int64_t) LASPoint_GetRawY(p)) + factorY);
566 *answer = (Expand1(x) << 1) + Expand1(y);
567
568 return *answer;
569 }
570
S64(const char * s)571 int64_t S64(const char *s) {
572 int64_t i;
573 char c ;
574 int scanned = sscanf(s, "lld%" SCNd64 "%c", &i, &c);
575 if (scanned == 1) return i;
576 fprintf(stderr, "ERROR: parsing string to int64_t.\n");
577 exit(1);
578 }
579
main(int argc,char * argv[])580 int main(int argc, char *argv[])
581 {
582 /*Initialize the catalog locks*/
583 MT_lock_init(&dataLock);
584 MT_cond_init(&mainCond);
585 MT_cond_init(&writeTCond);
586 MT_cond_init(&readCond);
587
588 char* file_name_in = 0;
589 char* file_name_out = 0;
590 char separator_sign = ' ';
591 char* parse_string = "xyz";
592 char* buffer;
593 char printstring[256];
594 LASReaderH reader = NULL;
595 LASHeaderH header = NULL;
596 LASPointH p = NULL;
597 FILE** files_out = NULL;
598 int len, j;
599 int64_t mortonkey = 0;
600 int num_files_in = 0, num_files, num_of_entries=0, check = 0, num_read_threads = DEFAULT_NUM_READ_THREADS;
601 int i;
602 pthread_t *writeThreads = NULL;
603 pthread_t *readThreads = NULL;
604 struct readThreadArgs *dataRead = NULL;
605 boolean input_file = FALSE;
606 int64_t global_offset_x = 0;
607 int64_t global_offset_y = 0;
608 double scale_x;
609 double scale_y;
610
611
612 if (argc == 1) {
613 usage();
614 exit(0);
615 }
616
617 /*Allocate space for input files*/
618 files_name_in = (char**) malloc(sizeof(char*)*DEFAULT_NUM_INPUT_FILES);
619
620 for (i = 1; i < argc; i++)
621 {
622 if ( strcmp(argv[i],"-h") == 0 ||
623 strcmp(argv[i],"-help") == 0 ||
624 strcmp(argv[i],"--help") == 0
625 )
626 {
627 usage();
628 exit(0);
629 }
630 else if ( strcmp(argv[i],"-v") == 0 ||
631 strcmp(argv[i],"--verbose") == 0
632 )
633 {
634 verbose = TRUE;
635 }
636 else if ( strcmp(argv[i],"--num_read_threads") == 0)
637 {
638 num_read_threads = atoi(argv[++i]);
639 }
640 else if ( strcmp(argv[i],"-s") == 0 ||
641 strcmp(argv[i],"--skip_invalid") == 0
642 )
643 {
644 skip_invalid = TRUE;
645 }
646 else if ( strcmp(argv[i], "--parse") == 0 ||
647 strcmp(argv[i], "-parse") == 0
648 )
649 {
650 i++;
651 if ( (parse_string = argv[i]) == NULL) {
652 usage();
653 exit(0);
654 }
655 }
656 else if ( strcmp(argv[i], "--moffset") == 0 ||
657 strcmp(argv[i], "-moffset") == 0
658 )
659 {
660 i++;
661 buffer = strtok (argv[i], ",");
662 j = 0;
663 while (buffer) {
664 if (j == 0) {
665 global_offset_x = S64(buffer);
666 }
667 else if (j == 1) {
668 global_offset_y = S64(buffer);
669 }
670 j++;
671 buffer = strtok (NULL, ",");
672 while (buffer && *buffer == '\040')
673 buffer++;
674 }
675 if (j != 2){
676 fprintf(stderr, "Only two int64_t are required in moffset option!\n");
677 exit(1);
678 }
679
680 }
681 else if ( strcmp(argv[i], "--check") == 0 ||
682 strcmp(argv[i], "-check") == 0
683 )
684 {
685 i++;
686 check = 1;
687 buffer = strtok (argv[i], ",");
688 j = 0;
689 while (buffer) {
690 if (j == 0) {
691 sscanf(buffer, "%lf", &scale_x);
692 }
693 else if (j == 1) {
694 sscanf(buffer, "%lf", &scale_y);
695 }
696 j++;
697 buffer = strtok (NULL, ",");
698 while (buffer && *buffer == '\040')
699 buffer++;
700 }
701 if (j != 2){
702 fprintf(stderr, "Only two doubles are required in moffset option!\n");
703 exit(1);
704 }
705 }
706 else if ( strcmp(argv[i],"--input") == 0 ||
707 strcmp(argv[i],"-input") == 0 ||
708 strcmp(argv[i],"-i") == 0 ||
709 strcmp(argv[i],"-in") == 0
710 )
711 {
712 i++;
713 files_name_in[num_files_in++] = argv[i];
714
715 if (num_files_in % DEFAULT_NUM_INPUT_FILES)
716 files_name_in = (char**) realloc(files_name_in, (num_files_in*2)*sizeof(char*));
717 }
718 else if (strcmp(argv[i],"--file") == 0 ||
719 strcmp(argv[i],"-file") == 0 ||
720 strcmp(argv[i],"-f") == 0
721 )
722 {
723 i++;
724 int read;
725 char line_buffer[BUFSIZ];
726 FILE* in = NULL;
727
728 in = fopen(argv[i], "r");
729 if (!in) {
730 fprintf(stderr, "ERROR: the path for file containing the input files is invalid %s\n", argv[i]);
731 exit(1);
732 }
733 while (fgets(line_buffer, sizeof(line_buffer), in)) {
734 line_buffer[strlen(line_buffer)-1]='\0';
735 files_name_in[num_files_in++] = strdup(line_buffer);
736
737 if (num_files_in % DEFAULT_NUM_INPUT_FILES)
738 files_name_in = (char**) realloc(files_name_in, (num_files_in*2)*sizeof(char*));
739 }
740 fclose(in);
741 input_file = TRUE;
742 }
743 else if ( strcmp(argv[i],"--output") == 0 ||
744 strcmp(argv[i],"--out") == 0 ||
745 strcmp(argv[i],"-out") == 0 ||
746 strcmp(argv[i],"-o") == 0
747 )
748 {
749 i++;
750 file_name_out = argv[i];
751 }
752 else
753 {
754 fprintf(stderr, "ERROR: unknown argument '%s'\n",argv[i]);
755 usage();
756 exit(1);
757 }
758 } /* end looping through argc/argv */
759 num_of_entries = strlen(parse_string);
760
761 if (num_files_in == 0)
762 {
763 LASError_Print("No input filename was specified");
764 usage();
765 exit(1);
766 }
767 files_name_in[num_files_in] = NULL;
768 num_files = num_files_in;
769
770 if (file_name_out == 0){
771 LASError_Print("No output prefix was specified");
772 usage();
773 exit(1);
774 }
775
776 /*Entries metadata*/
777 i = 0;
778 for (;;)
779 {
780 switch (parse_string[i])
781 {
782 /* // the morton code on xy */
783 case 'k':
784 entries[i] = ENTRY_k;
785 entriesType[i] = sizeof(int64_t);
786 /*Changes for Oscar's new Morton code function*/
787 //entriesFunc[i] = (void*)morton2D_encode;
788 entriesFuncD[i] = (void*)morton2D_encodeOscar;
789 break;
790 /* // the x coordinate double*/
791 case 'x':
792 entries[i] = ENTRY_x;
793 entriesType[i] = sizeof(double);
794 entriesFuncD[i] = (void*)LASPoint_GetX;
795 break;
796 /* // the y coordinate double*/
797 case 'y':
798 entries[i] = ENTRY_y;
799 entriesType[i] = sizeof(double);
800 entriesFuncD[i] = (void*)LASPoint_GetY;
801 break;
802 /* // the z coordinate double*/
803 case 'z':
804 entries[i] = ENTRY_z;
805 entriesType[i] = sizeof(double);
806 entriesFuncD[i] = (void*)LASPoint_GetZ;
807 break;
808 /* // the X coordinate decimal*/
809 case 'X':
810 entries[i] = ENTRY_X;
811 entriesType[i] = sizeof(int);
812 entriesFuncD[i] = (void*)LASPoint_GetX;
813 break;
814 /* // the y coordinate decimal*/
815 case 'Y':
816 entries[i] = ENTRY_Y;
817 entriesType[i] = sizeof(int);
818 entriesFuncD[i] = (void*)LASPoint_GetY;
819 break;
820 /* // the z coordinate decimal*/
821 case 'Z':
822 entries[i] = ENTRY_Z;
823 entriesType[i] = sizeof(int);
824 entriesFuncD[i] = (void*)LASPoint_GetZ;
825 break;
826 /* // the gps-time */
827 case 't':
828 entries[i] = ENTRY_t;
829 entriesType[i] = sizeof(double);
830 entriesFuncD[i] = (void*)LASPoint_GetTime;
831 break;
832 /* // the intensity */
833 case 'i':
834 entries[i] = ENTRY_i;
835 entriesType[i] = sizeof(unsigned short);
836 entriesFuncS[i] = (void*)LASPoint_GetIntensity;
837 break;
838 /* the scan angle */
839 case 'a':
840 entries[i] = ENTRY_a;
841 entriesType[i] = sizeof(char);
842 entriesFuncC[i] = (void*)LASPoint_GetScanAngleRank;
843 break;
844 /* the number of the return */
845 case 'r':
846 entries[i] = ENTRY_r;
847 entriesType[i] = sizeof(short);
848 entriesFuncS[i] = (void*)LASPoint_GetReturnNumber;
849 break;
850 /* the classification */
851 case 'c':
852 entries[i] = ENTRY_c;
853 entriesType[i] = sizeof(char);
854 entriesFuncC[i] = (void*)LASPoint_GetClassification;
855 break;
856 /* the user data */
857 case 'u':
858 entries[i] = ENTRY_u;
859 entriesType[i] = sizeof(char);
860 entriesFuncC[i] = (void*)LASPoint_GetUserData;
861 break;
862 /* the number of returns of given pulse */
863 case 'n':
864 entries[i] = ENTRY_n;
865 entriesType[i] = sizeof(short);
866 entriesFuncS[i] = (void*)LASPoint_GetNumberOfReturns;
867 break;
868 /* the red channel color */
869 case 'R':
870 entries[i] = ENTRY_R;
871 entriesType[i] = sizeof(unsigned short);
872 entriesFuncS[i] = (void*)LASColor_GetRed;
873 break;
874 /* the green channel color */
875 case 'G':
876 entries[i] = ENTRY_G;
877 entriesType[i] = sizeof(unsigned short);
878 entriesFuncS[i] = (void*)LASColor_GetGreen;
879 break;
880 /* the blue channel color */
881 case 'B':
882 entries[i] = ENTRY_B;
883 entriesType[i] = sizeof(unsigned short);
884 entriesFuncS[i] = (void*)LASColor_GetBlue;
885 break;
886 case 'M':
887 entries[i] = ENTRY_M;
888 entriesType[i] = sizeof(int);
889 break;
890 case 'p':
891 entries[i] = ENTRY_p;
892 entriesType[i] = sizeof(short);
893 entriesFuncS[i] = (void*)LASPoint_GetPointSourceId;
894 break;
895 /* the edge of flight line flag */
896 case 'e':
897 entries[i] = ENTRY_e;
898 entriesType[i] = sizeof(short);
899 entriesFuncS[i] = (void*)LASPoint_GetFlightLineEdge;
900 break;
901 /* the direction of scan flag */
902 case 'd':
903 entries[i] = ENTRY_d;
904 entriesType[i] = sizeof(short);
905 entriesFuncS[i] = (void*)LASPoint_GetScanDirection;
906 break;
907 }
908 i++;
909 if (parse_string[i] == 0)
910 {
911 break;
912 }
913 }
914
915 /*Prepare the output files*/
916 if (file_name_out == NULL)
917 {
918 len = (int)strlen(file_name_in);
919 file_name_out = LASCopyString(file_name_in);
920 if (file_name_out[len-3] == '.' && file_name_out[len-2] == 'g' && file_name_out[len-1] == 'z')
921 {
922 len = len - 4;
923 }
924 while (len > 0 && file_name_out[len] != '.')
925 {
926 len--;
927 }
928 file_name_out[len] = '\0';
929 }
930
931 char *str = malloc(sizeof(char)*(strlen(file_name_out)+12));
932 files_out = (FILE**) malloc(sizeof(FILE*)*num_of_entries);
933 for (i = 0; i < num_of_entries; i++) {
934 sprintf(str, "%s_col_%c.dat", file_name_out, parse_string[i]);
935 if(doesFileExist(str)) {
936 remove(str);
937 }
938
939 files_out[i] = fopen(str, "wb");
940
941 if (files_out[i] == 0) {
942 LASError_Print("Could not open file for write");
943 usage();
944 exit(1);
945 }
946 }
947 free(str);
948
949 /*Initialize structures for the reading threads*/
950 //data = (struct writeT**) malloc(num_read_threads*sizeof(struct writeT*)); //Malloc is more efficient than calloc
951 data = (struct writeT**) calloc(num_read_threads, sizeof(struct writeT*));
952
953 dataRead = (struct readThreadArgs*) malloc(sizeof(struct readThreadArgs)*num_read_threads);
954 /* Launch read Threads */
955 stop = 0;
956 readThreads = (pthread_t*) malloc(sizeof(pthread_t)*num_read_threads);
957 for (i=0; i < num_read_threads; i++) {
958 dataRead[i].id = i;
959 dataRead[i].num_read_threads = num_read_threads;
960 dataRead[i].num_of_entries = num_of_entries;
961 dataRead[i].check = check;
962 dataRead[i].global_offset_x = global_offset_x;
963 dataRead[i].global_offset_y = global_offset_y;
964 dataRead[i].scale_x = scale_x;
965 dataRead[i].scale_y = scale_y;
966 pthread_create(&readThreads[i], NULL, readFile, (void*)dataRead);
967 }
968
969 int writeIndex = 0;
970 writeThreads = (pthread_t*) malloc(sizeof(pthread_t)*num_of_entries);
971
972 /* Launch Threads */
973 struct writeThreadArgs *dataWrite = (struct writeThreadArgs *) malloc(sizeof(struct writeThreadArgs) *num_of_entries);
974 for (i = 0; i < num_of_entries; i++) {
975 dataWrite[i].id = i;
976 dataWrite[i].out = files_out[i];
977 pthread_create(&writeThreads[i], NULL, writeFile, (void*)(&dataWrite[i]));
978 }
979 sleep(1);
980 //Do we need to comment this one out!?
981 int done = 0;
982 while (num_files) {
983 /*Obtain lock over data to get the pointer*/
984 MT_set_lock(&dataLock);
985 dataWriteT = data[writeIndex];
986 while (dataWriteT == NULL) {
987 /*Sleep and wait for data to be read*/
988 MT_cond_wait(&mainCond,&dataLock);
989 dataWriteT = data[writeIndex];
990 }
991 data[writeIndex] = NULL;
992 //Release the lock
993
994 /*Tell the write threads there is new data*/
995 pthread_cond_broadcast(&writeTCond);
996
997 /*Tell the read threads there is a new buf empty*/
998 pthread_cond_broadcast(&readCond);
999 MT_unset_lock(&dataLock);
1000
1001 /*Keep looping*/
1002 writeIndex++;
1003 writeIndex = (writeIndex % num_read_threads);
1004
1005 MT_set_lock(&dataLock);
1006 while (done == 0) {
1007 /*Sleep and wait for data to be read*/
1008 MT_cond_wait(&mainCond,&dataLock);
1009 done = 1;
1010 for (i = 0; i < num_of_entries; i++) {
1011 if (dataWriteT[i].values != NULL) {
1012 done = 0;
1013 break;
1014 }
1015 }
1016 }
1017 num_files--;
1018 if (verbose)
1019 printf("Files to go %d\n", num_files);
1020 free(dataWriteT);
1021 dataWriteT = NULL;
1022 done = 0;
1023 MT_unset_lock(&dataLock);
1024 }
1025
1026 /*Tell the write threads to exit*/
1027 MT_set_lock(&dataLock);
1028 stop = 1;
1029 pthread_cond_broadcast(&writeTCond);
1030 MT_unset_lock(&dataLock);
1031
1032 /* Wait for Threads to Finish */
1033 for (i=0; i<num_of_entries; i++) {
1034 pthread_join(writeThreads[i], NULL);
1035 }
1036 free(dataWrite);
1037 free(writeThreads);
1038
1039 MT_cond_destroy(&readCond);
1040 MT_cond_destroy(&writeTCond);
1041 MT_cond_destroy(&mainCond);
1042 MT_lock_destroy(&dataLock);
1043
1044 for (i = 0; i < num_of_entries; i++) {
1045 fflush(files_out[i]);
1046 if (verbose)
1047 printf("close file %d\n", i);
1048 fsync(files_out[i]);
1049 fclose(files_out[i]);
1050 }
1051 free(files_out);
1052 if (input_file) {
1053 for (i=0 ; i < num_files_in; i++)
1054 free(files_name_in[i]);
1055
1056 free(files_name_in);
1057 }
1058
1059 free(dataRead);
1060
1061 if (readThreads)
1062 free(readThreads);
1063
1064 return 0;
1065 }
1066