Context Navigation

rt_api.c

Visit:

Revision 5865, 85.2 KB (checked in by pracine, 21 months ago)
-Fix segfault introduced in precedent commit. Better copy of old band values.
Property svn:keywords set to `Id`

Line
1	/*
2	* $Id$
3	*
4	* WKTRaster - Raster Types for PostGIS
5	* http://www.postgis.org/support/wiki/index.php?WKTRasterHomePage
6	*
7	* Copyright (C) 2009 Sandro Santilli <strk@keybit.net>, Pierre Racine <pierre.racine@sbf.ulaval.ca>,
8	* Jorge Arevalo <jorge.arevalo@deimos-space.com>
9	*
10	* This program is free software; you can redistribute it and/or modify
11	* it under the terms of the GNU General Public License as published by
12	* the Free Software Foundation; either version 3 of the License, or
13	* (at your option) any later version.
14	*
15	* This program is distributed in the hope that it will be useful,
16	* but WITHOUT ANY WARRANTY; without even the implied warranty of
17	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18	* GNU General Public License for more details.
19	*
20	* You should have received a copy of the GNU General Public License
21	* along with this program; if not, write to the Free Software
22	* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
23	*
24	*/
25
26	#include <math.h>
27	#include <stdio.h> /* for printf (default message handler) */
28	#include <stdarg.h> /* for va_list, va_start etc */
29	#include <string.h> /* for memcpy */
30	#include <assert.h>
31	#include <float.h> /* for FLT_EPSILON */
32
33	#include "rt_api.h"
34
35	/* Define this to debug low-level RASTER API activity */
36	/* Alternative, use ./configure --enable-rtapi-debug */
37	/#define POSTGIS_RASTER_API_DEBUG/
38
39	#define POSTGIS_RASTER_WARN_ON_TRUNCATION
40
41	/--- Utilities -------------------------------------------------/
42
43	static void
44	swap_char(uint8_t a, uint8_t b)
45	{
46	uint8_t c = 0;
47
48	assert(NULL != a && NULL != b);
49
50	c = *a;
51	a=b;
52	*b=c;
53	}
54
55	static void
56	flip_endian_16(uint8_t *d)
57	{
58	assert(NULL != d);
59
60	swap_char(d, d+1);
61	}
62
63	static void
64	flip_endian_32(uint8_t *d)
65	{
66	assert(NULL != d);
67
68	swap_char(d, d+3);
69	swap_char(d+1, d+2);
70	}
71
72	static void
73	flip_endian_64(uint8_t *d)
74	{
75	assert(NULL != d);
76
77	swap_char(d+7, d);
78	swap_char(d+6, d+1);
79	swap_char(d+5, d+2);
80	swap_char(d+4, d+3);
81	}
82
83	/- rt_context -------------------------------------------------------/
84
85	static void
86	default_error_handler(const char *fmt, ...)
87	{
88	va_list ap;
89
90	static const char *label = "ERROR: ";
91	char newfmt[1024] = { 0 };
92	snprintf(newfmt, 1024, "%s%s\n", label, fmt);
93	newfmt[1023] = '\0';
94
95	va_start(ap, fmt);
96
97	vprintf(newfmt, ap);
98
99	va_end(ap);
100	}
101
102	static void
103	default_warning_handler(const char *fmt, ...)
104	{
105	va_list ap;
106
107	static const char *label = "WARNING: ";
108	char newfmt[1024] = { 0 };
109	snprintf(newfmt, 1024, "%s%s\n", label, fmt);
110	newfmt[1023] = '\0';
111
112	va_start(ap, fmt);
113
114	vprintf(newfmt, ap);
115
116	va_end(ap);
117	}
118
119	static void
120	default_info_handler(const char *fmt, ...)
121	{
122	va_list ap;
123
124	static const char *label = "INFO: ";
125	char newfmt[1024] = { 0 };
126	snprintf(newfmt, 1024, "%s%s\n", label, fmt);
127	newfmt[1023] = '\0';
128
129	va_start(ap, fmt);
130
131	vprintf(newfmt, ap);
132
133	va_end(ap);
134	}
135
136	struct rt_context_t {
137	rt_allocator alloc;
138	rt_reallocator realloc;
139	rt_deallocator dealloc;
140	rt_message_handler err;
141	rt_message_handler warn;
142	rt_message_handler info;
143	};
144
145	rt_context
146	rt_context_new(rt_allocator allocator, rt_reallocator reallocator,
147	rt_deallocator deallocator)
148	{
149	rt_context ret;
150
151	if ( ! allocator ) allocator = malloc;
152	if ( ! reallocator ) reallocator = realloc;
153	if ( ! deallocator ) deallocator = free;
154
155	ret = (rt_context)allocator(sizeof(struct rt_context_t));
156	if ( ! ret ) {
157	default_error_handler("Out of virtual memory creating an rt_context");
158	return 0;
159	}
160	#ifdef POSTGIS_RASTER_API_DEBUG
161	default_info_handler("Created rt_context @ %p", ret);
162	#endif
163	ret->alloc = allocator;
164	ret->realloc = reallocator;
165	ret->dealloc = deallocator;
166	ret->err = default_error_handler;
167	ret->warn = default_warning_handler;
168	ret->info = default_info_handler;
169
170	assert(NULL != ret->alloc);
171	assert(NULL != ret->realloc);
172	assert(NULL != ret->dealloc);
173	assert(NULL != ret->err);
174	assert(NULL != ret->warn);
175	assert(NULL != ret->info);
176	return ret;
177	}
178
179	void
180	rt_context_set_message_handlers(rt_context ctx,
181	rt_message_handler error_handler,
182	rt_message_handler warning_handler,
183	rt_message_handler info_handler)
184	{
185	ctx->err = error_handler;
186	ctx->warn = warning_handler;
187	ctx->info = info_handler;
188
189	assert(NULL != ctx->err);
190	assert(NULL != ctx->warn);
191	assert(NULL != ctx->info);
192	}
193
194	void
195	rt_context_destroy(rt_context ctx)
196	{
197	#ifdef POSTGIS_RASTER_API_DEBUG
198	ctx->info("Destroying rt_context @ %p", ctx);
199	#endif
200	ctx->dealloc(ctx);
201	}
202
203	/--- Debug and Testing Utilities --------------------------------------------/
204
205	#ifdef POSTGIS_RASTER_API_DEBUG
206
207	static char*
208	d_binary_to_hex(rt_context ctx, const uint8_t* const raw, uint32_t size, uint32_t *hexsize)
209	{
210	char* hex = NULL;
211	uint32_t i = 0;
212
213	assert(NULL != ctx);
214	assert(NULL != raw);
215	assert(NULL != hexsize);
216
217	hexsize = size 2; /* hex is 2 times bytes */
218	hex = (char)ctx->alloc( (hexsize) + 1);
219	if ( ! hex )
220	{
221	ctx->err("Out of memory hexifying raw binary\n");
222	return NULL;
223	}
224	hex[hexsize] = '\0'; / Null-terminate */
225
226	for (i = 0; i < size; ++i)
227	{
228	deparse_hex(raw[i], &(hex[2 * i]));
229	}
230
231	assert(NULL != hex);
232	assert(0 == strlen(hex) % 2);
233	return hex;
234	}
235
236	static void
237	d_print_binary_hex(rt_context ctx, const char* msg, const uint8_t* const raw, uint32_t size)
238	{
239	char* hex = NULL;
240	uint32_t hexsize = 0;
241
242	assert(NULL != ctx);
243	assert(NULL != msg);
244	assert(NULL != raw);
245
246	hex = d_binary_to_hex(ctx, raw, size, &hexsize);
247	if (NULL != hex)
248	{
249	ctx->info("%s\t%s", msg, hex);
250	ctx->dealloc(hex);
251	}
252	}
253
254	static size_t
255	d_binptr_to_pos(const uint8_t* const ptr, const uint8_t* const end, size_t size)
256	{
257	assert(NULL != ptr && NULL != end);
258
259	return (size - (end - ptr));
260	}
261
262	#define CHECK_BINPTR_POSITION(ptr, end, size, pos) \
263	{ if (pos != d_binptr_to_pos(ptr, end, size)) { \
264	fprintf(stderr, "Check of binary stream pointer position failed on line %d\n", __LINE__); \
265	fprintf(stderr, "\tactual = %lu, expected = %lu\n", (long unsigned)d_binptr_to_pos(ptr, end, size), (long unsigned)pos); \
266	} }
267
268	#else
269
270	#define CHECK_BINPTR_POSITION(ptr, end, size, pos) ((void)0);
271
272	#endif /* ifndef POSTGIS_RASTER_API_DEBUG */
273
274	/- rt_pixeltype -----------------------------------------------------/
275
276	int
277	rt_pixtype_size(rt_context ctx, rt_pixtype pixtype)
278	{
279	int pixbytes = -1;
280
281	assert(NULL != ctx);
282
283	switch (pixtype)
284	{
285	case PT_1BB:
286	case PT_2BUI:
287	case PT_4BUI:
288	case PT_8BSI:
289	case PT_8BUI:
290	pixbytes = 1;
291	break;
292	case PT_16BSI:
293	case PT_16BUI:
294	pixbytes = 2;
295	break;
296	case PT_32BSI:
297	case PT_32BUI:
298	case PT_32BF:
299	pixbytes = 4;
300	break;
301	case PT_64BF:
302	pixbytes = 8;
303	break;
304	default:
305	ctx->err("Unknown pixeltype %d", pixtype);
306	pixbytes = -1;
307	break;
308	}
309
310	#ifdef POSTGIS_RASTER_API_DEBUG
311	ctx->info("Pixel type = %s and size = %d bytes",
312	rt_pixtype_name(ctx, pixtype), pixbytes);
313	#endif
314
315	return pixbytes;
316	}
317
318	int
319	rt_pixtype_alignment(rt_context ctx, rt_pixtype pixtype)
320	{
321	return rt_pixtype_size(ctx, pixtype);
322	}
323
324	rt_pixtype
325	rt_pixtype_index_from_name(rt_context ctx, const char* pixname)
326	{
327	assert(strlen(pixname) > 0);
328
329	if (strcmp(pixname, "1BB") == 0)
330	return PT_1BB;
331	if (strcmp(pixname, "2BUI") == 0)
332	return PT_2BUI;
333	if (strcmp(pixname, "4BUI") == 0)
334	return PT_4BUI;
335	if (strcmp(pixname, "8BSI") == 0)
336	return PT_8BSI;
337	if (strcmp(pixname, "8BUI") == 0)
338	return PT_8BUI;
339	if (strcmp(pixname, "16BSI") == 0)
340	return PT_16BSI;
341	if (strcmp(pixname, "16BUI") == 0)
342	return PT_16BUI;
343	if (strcmp(pixname, "32BSI") == 0)
344	return PT_32BSI;
345	if (strcmp(pixname, "32BUI") == 0)
346	return PT_32BUI;
347	if (strcmp(pixname, "32BF") == 0)
348	return PT_32BF;
349	if (strcmp(pixname, "64BF") == 0)
350	return PT_64BF;
351	return PT_END;
352	}
353
354	const char*
355	rt_pixtype_name(rt_context ctx, rt_pixtype pixtype)
356	{
357	assert(NULL != ctx);
358
359	switch (pixtype)
360	{
361	case PT_1BB:
362	return "1BB";
363	case PT_2BUI:
364	return "2BUI";
365	case PT_4BUI:
366	return "4BUI";
367	case PT_8BSI:
368	return "8BSI";
369	case PT_8BUI:
370	return "8BUI";
371	case PT_16BSI:
372	return "16BSI";
373	case PT_16BUI:
374	return "16BUI";
375	case PT_32BSI:
376	return "32BSI";
377	case PT_32BUI:
378	return "32BUI";
379	case PT_32BF:
380	return "32BF";
381	case PT_64BF:
382	return "64BF";
383	default:
384	ctx->err("Unknown pixeltype %d", pixtype);
385	return "Unknown";
386	}
387	}
388
389	/- rt_band ----------------------------------------------------------/
390
391	struct rt_extband_t {
392	uint8_t bandNum;
393	char* path; /* externally owned ? */
394	};
395
396	struct rt_band_t {
397
398	rt_pixtype pixtype;
399	int32_t offline;
400	uint16_t width;
401	uint16_t height;
402	int32_t hasnodata; /* a flag indicating if this band contains nodata values */
403	double nodataval; /* int will be converted ... */
404	int32_t ownsData; /* XXX mloskot: its behaviour needs to be documented */
405	union {
406	void* mem; /* actual data, externally owned */
407	struct rt_extband_t offline;
408	} data;
409
410	};
411
412	rt_band
413	rt_band_new_inline(rt_context ctx, uint16_t width, uint16_t height,
414	rt_pixtype pixtype, uint32_t hasnodata, double nodataval,
415	uint8_t* data)
416	{
417	rt_band band = NULL;
418
419	assert(NULL != ctx);
420	assert(NULL != data);
421
422	band = ctx->alloc(sizeof(struct rt_band_t));
423	if ( ! band )
424	{
425	ctx->err("Out of memory allocating rt_band");
426	return 0;
427	}
428
429	#ifdef POSTGIS_RASTER_API_DEBUG
430	ctx->info("Created rt_band @ %p with pixtype %s",
431	band, rt_pixtype_name(ctx, pixtype));
432	#endif
433
434	band->pixtype = pixtype;
435	band->offline = 0;
436	band->width = width;
437	band->height = height;
438	band->hasnodata = hasnodata;
439	band->nodataval = nodataval;
440	band->data.mem = data;
441	band->ownsData = 0;
442
443	return band;
444	}
445
446	rt_band
447	rt_band_new_offline(rt_context ctx, uint16_t width, uint16_t height,
448	rt_pixtype pixtype, uint32_t hasnodata, double nodataval,
449	uint8_t bandNum, const char* path)
450	{
451	rt_band band = NULL;
452
453	assert(NULL != ctx);
454	assert(NULL != path);
455
456	band = ctx->alloc(sizeof(struct rt_band_t));
457	if ( ! band )
458	{
459	ctx->err("Out of memory allocating rt_band");
460	return 0;
461	}
462
463	#ifdef POSTGIS_RASTER_API_DEBUG
464	ctx->info("Created rt_band @ %p with pixtype %s",
465	band, rt_pixtype_name(ctx, pixtype));
466	#endif
467
468	band->pixtype = pixtype;
469	band->offline = 1;
470	band->width = width;
471	band->height = height;
472	band->hasnodata = hasnodata;
473	band->nodataval = nodataval;
474	band->data.offline.bandNum = bandNum;
475
476	/* memory for data.offline.path should be managed externally */
477	band->data.offline.path = (char *)path;
478
479	/* XXX QUESTION (jorgearevalo): What does exactly ownsData mean?? I think that
480	* ownsData = 0 ==> the memory for band->data is externally owned
481	* ownsData = 1 ==> the memory for band->data is internally owned
482	*/
483	band->ownsData = 0;
484
485	return band;
486	}
487
488
489
490	int
491	rt_band_is_offline(rt_context ctx, rt_band band)
492	{
493	assert(NULL != ctx);
494	assert(NULL != band);
495
496	return band->offline;
497	}
498
499	void
500	rt_band_destroy(rt_context ctx, rt_band band)
501	{
502	#ifdef POSTGIS_RASTER_API_DEBUG
503	ctx->info("Destroying rt_band @ %p", band);
504	#endif
505
506	/* band->data content is externally owned */
507	/* XXX jorgearevalo: not really... rt_band_from_wkb allocates memory for
508	* data.me
509	*/
510	ctx->dealloc(band);
511	}
512
513	const char*
514	rt_band_get_ext_path(rt_context ctx, rt_band band)
515	{
516	assert(NULL != ctx);
517	assert(NULL != band);
518
519	if ( ! band->offline ) {
520	ctx->warn("rt_band_get_ext_path: non-offline band doesn't have "
521	"an associated path");
522	return 0;
523	}
524	return band->data.offline.path;
525	}
526
527	uint8_t
528	rt_band_get_ext_band_num(rt_context ctx, rt_band band)
529	{
530	assert(NULL != ctx);
531	assert(NULL != band);
532
533	if ( ! band->offline ) {
534	ctx->warn("rt_band_get_ext_path: non-offline band doesn't have "
535	"an associated band number");
536	return 0;
537	}
538	return band->data.offline.bandNum;
539	}
540
541	void *
542	rt_band_get_data(rt_context ctx, rt_band band)
543	{
544	assert(NULL != ctx);
545	assert(NULL != band);
546
547	if ( band->offline )
548	{
549	ctx->warn("rt_band_get_data: "
550	"offline band doesn't have associated data");
551	return 0;
552	}
553	return band->data.mem;
554	}
555
556	rt_pixtype
557	rt_band_get_pixtype(rt_context ctx, rt_band band)
558	{
559	assert(NULL != ctx);
560	assert(NULL != band);
561
562	return band->pixtype;
563	}
564
565	uint16_t
566	rt_band_get_width(rt_context ctx, rt_band band)
567	{
568	assert(NULL != ctx);
569	assert(NULL != band);
570
571	return band->width;
572	}
573
574	uint16_t
575	rt_band_get_height(rt_context ctx, rt_band band)
576	{
577	assert(NULL != ctx);
578	assert(NULL != band);
579
580	return band->height;
581	}
582
583	#ifdef OPTIMIZE_SPACE
584	/*
585	* Set given number of bits of the given byte,
586	* starting from given bitOffset (from the first)
587	* to the given value.
588	*
589	* Examples:
590	* char ch;
591	* ch=0; setBits(&ch, 1, 1, 0) -> ch==8
592	* ch=0; setBits(&ch, 3, 2, 1) -> ch==96 (0x60)
593	*
594	* Note that number of bits set must be <= 8-bitOffset
595	*
596	*/
597	static void
598	setBits(char* ch, double val, int bits, int bitOffset)
599	{
600	char mask = 0xFF >> (8 - bits);
601	char ival = val;
602
603	assert(8-bitOffset >= bits);
604
605	#ifdef POSTGIS_RASTER_API_DEBUG
606	printf("ival:%d bits:%d mask:%hhx bitoffset:%d\n",
607	ival, bits, mask, bitOffset);
608	#endif
609
610	/* clear all but significant bits from ival */
611	ival &= mask;
612	#ifdef POSTGIS_RASTER_WARN_ON_TRUNCATION
613	if ( ival != val ) {
614	ctx->warn("Pixel value for %d-bits band got truncated"
615	" from %g to %hhu\n", bits, val, ival);
616	}
617	#endif /* POSTGIS_RASTER_WARN_ON_TRUNCATION */
618
619	#ifdef POSTGIS_RASTER_API_DEBUG
620	printf(" cleared ival:%hhx\n", ival);
621	#endif
622
623	/* Shift ival so the significant bits start at
624	* the first bit */
625	ival <<= (8-bitOffset-bits);
626
627	#ifdef POSTGIS_RASTER_API_DEBUG
628	printf(" ival shifted:%hhx\n", ival);
629	printf(" ch:%hhx\n", *ch);
630	#endif
631
632	/* clear first bits of target */
633	*ch &= ~(mask << (8-bits-bitOffset));
634
635	#ifdef POSTGIS_RASTER_API_DEBUG
636	printf(" ch cleared:%hhx\n", *ch);
637	#endif
638
639	/* Set the first bit of target */
640	*ch \|= ival;
641
642	#ifdef POSTGIS_RASTER_API_DEBUG
643	printf(" ch ored:%hhx\n", *ch);
644	#endif
645
646	}
647	#endif /* OPTIMIZE_SPACE */
648
649	int
650	rt_band_get_hasnodata_flag(rt_context ctx, rt_band band)
651	{
652	assert(NULL != ctx);
653	assert(NULL != band);
654
655	return band->hasnodata;
656	}
657
658	void
659	rt_band_set_hasnodata_flag(rt_context ctx, rt_band band, int flag)
660	{
661	assert(NULL != ctx);
662	assert(NULL != band);
663
664	band->hasnodata = (flag) ? 1 : 0;
665	}
666
667	int
668	rt_band_set_nodata(rt_context ctx, rt_band band, double val)
669	{
670	rt_pixtype pixtype = PT_END;
671
672	assert(NULL != ctx);
673	assert(NULL != band);
674
675	pixtype = band->pixtype;
676
677	#ifdef POSTGIS_RASTER_API_DEBUG
678	ctx->info("rt_band_set_nodata: setting NODATA %g with band type %s", val, rt_pixtype_name(ctx, pixtype));
679	#endif
680
681	/* return -1 on out of range */
682	switch (pixtype)
683	{
684	case PT_1BB:
685	{
686	uint8_t v = val;
687	v &= 0x01;
688	band->nodataval = v;
689	break;
690	}
691	case PT_2BUI:
692	{
693	uint8_t v = val;
694	v &= 0x03;
695	band->nodataval = v;
696	break;
697	}
698	case PT_4BUI:
699	{
700	uint8_t v = val;
701	v &= 0x0F;
702	band->nodataval = v;
703	break;
704	}
705	case PT_8BSI:
706	{
707	int8_t v = val;
708	band->nodataval = v;
709	break;
710	}
711	case PT_8BUI:
712	{
713	uint8_t v = val;
714	band->nodataval = v;
715	break;
716	}
717	case PT_16BSI:
718	{
719	int16_t v = val;
720	band->nodataval = v;
721	break;
722	}
723	case PT_16BUI:
724	{
725	uint16_t v = val;
726	band->nodataval = v;
727	break;
728	}
729	case PT_32BSI:
730	{
731	int32_t v = val;
732	band->nodataval = v;
733	break;
734	}
735	case PT_32BUI:
736	{
737	uint32_t v = val;
738	band->nodataval = v;
739	break;
740	}
741	case PT_32BF:
742	{
743	float v = val;
744	band->nodataval = v;
745	break;
746	}
747	case PT_64BF:
748	{
749	band->nodataval = val;
750	break;
751	}
752	default:
753	{
754	ctx->err("Unknown pixeltype %d", pixtype);
755	band->hasnodata = 0;
756	return -1;
757	}
758	}
759
760	#ifdef POSTGIS_RASTER_API_DEBUG
761	ctx->info("rt_band_set_nodata: band->hasnodata = %d", band->hasnodata);
762	ctx->info("rt_band_set_nodata: band->nodataval = %g", band->nodataval);
763	#endif
764
765
766	// the NODATA value was just set, so this band has NODATA
767	rt_band_set_hasnodata_flag(ctx, band, 1);
768
769	if ( fabs(band->nodataval - val) > 0.0001 )
770	{
771	#ifdef POSTGIS_RASTER_WARN_ON_TRUNCATION
772	ctx->warn("rt_band_set_nodata: NODATA value for %s band got truncated"
773	" from %g to %g",
774	rt_pixtype_name(ctx, pixtype),
775	val, band->nodataval);
776	#endif
777	return -1;
778	}
779
780	return 0;
781	}
782
783	int
784	rt_band_set_pixel(rt_context ctx, rt_band band, uint16_t x, uint16_t y,
785	double val)
786	{
787	rt_pixtype pixtype = PT_END;
788	unsigned char* data = NULL;
789	uint32_t offset = 0;
790	double checkval = 0;
791
792	assert(NULL != ctx);
793	assert(NULL != band);
794
795	pixtype = band->pixtype;
796
797	if ( x >= band->width \|\| y >= band->height )
798	{
799	ctx->err("Coordinates ouf of range");
800	return -1;
801	}
802
803	if ( band->offline )
804	{
805	ctx->err("rt_band_set_pixel not implemented yet for OFFDB bands");
806	return -1;
807	}
808
809	data = rt_band_get_data(ctx, band);
810	offset = x + (y * band->width);
811
812	switch (pixtype)
813	{
814	case PT_1BB:
815	{
816	data[offset] = (int)val&0x01;
817	checkval = data[offset];
818	break;
819	}
820	case PT_2BUI:
821	{
822	data[offset] = (int)val&0x03;
823	checkval = data[offset];
824	break;
825	}
826	case PT_4BUI:
827	{
828	data[offset] = (int)val&0x0F;
829	checkval = data[offset];
830	break;
831	}
832	case PT_8BSI:
833	{
834	data[offset] = (int8_t)val;
835	checkval = (int8_t)data[offset];
836	break;
837	}
838	case PT_8BUI:
839	{
840	data[offset] = val;
841	checkval = data[offset];
842	break;
843	}
844	case PT_16BSI:
845	{
846	int16_t ptr = (int16_t)data; /* we assume correct alignment */
847	ptr[offset] = val;
848	checkval = (int16_t)ptr[offset];
849	break;
850	}
851	case PT_16BUI:
852	{
853	uint16_t ptr = (uint16_t)data; /* we assume correct alignment */
854	ptr[offset] = val;
855	checkval = ptr[offset];
856	break;
857	}
858	case PT_32BSI:
859	{
860	int32_t ptr = (int32_t)data; /* we assume correct alignment */
861	ptr[offset] = val;
862	checkval = (int32_t)ptr[offset];
863	break;
864	}
865	case PT_32BUI:
866	{
867	uint32_t ptr = (uint32_t)data; /* we assume correct alignment */
868	ptr[offset] = val;
869	checkval = ptr[offset];
870	break;
871	}
872	case PT_32BF:
873	{
874	float ptr = (float)data; /* we assume correct alignment */
875	ptr[offset] = val;
876	checkval = ptr[offset];
877	break;
878	}
879	case PT_64BF:
880	{
881	double ptr = (double)data; /* we assume correct alignment */
882	ptr[offset] = val;
883	checkval = ptr[offset];
884	break;
885	}
886	default:
887	{
888	ctx->err("Unknown pixeltype %d", pixtype);
889	return -1;
890	}
891	}
892
893	/* Overflow checking */
894	if ( fabs(checkval-val) > 0.0001 )
895	{
896	#ifdef POSTGIS_RASTER_WARN_ON_TRUNCATION
897	ctx->warn("Pixel value for %s band got truncated"
898	" from %g to %g",
899	rt_pixtype_name(ctx, band->pixtype),
900	val, checkval);
901	#endif /* POSTGIS_RASTER_WARN_ON_TRUNCATION */
902	return -1;
903	}
904
905	return 0;
906	}
907
908	int
909	rt_band_get_pixel(rt_context ctx, rt_band band, uint16_t x, uint16_t y, double *result)
910	{
911	rt_pixtype pixtype = PT_END;
912	uint8_t* data = NULL;
913	uint32_t offset = 0;
914
915	assert(NULL != ctx);
916	assert(NULL != band);
917
918	pixtype = band->pixtype;
919
920	if ( x >= band->width \|\| y >= band->height )
921	{
922	ctx->warn("Attempting to get pixel value with out of range raster coordinates");
923	return -1;
924	}
925
926	if ( band->offline )
927	{
928	ctx->err("rt_band_get_pixel not implemented yet for OFFDB bands");
929	return -1;
930	}
931
932	data = rt_band_get_data(ctx, band);
933	offset = x + (y * band->width); /* +1 for the nodata value */
934
935	switch (pixtype)
936	{
937	case PT_1BB:
938	#ifdef OPTIMIZE_SPACE
939	{
940	int byteOffset = offset/8;
941	int bitOffset = offset%8;
942	data += byteOffset;
943
944	/* Bit to set is bitOffset into data */
945	*result = getBits(data, val, 1, bitOffset);
946	return 0;
947	}
948	#endif
949
950	case PT_2BUI:
951	#ifdef OPTIMIZE_SPACE
952	{
953	int byteOffset = offset/4;
954	int bitOffset = offset%4;
955	data += byteOffset;
956
957	/* Bits to set start at bitOffset into data */
958	*result = getBits(data, val, 2, bitOffset);
959	return 0;
960	}
961	#endif
962
963	case PT_4BUI:
964	#ifdef OPTIMIZE_SPACE
965	{
966	int byteOffset = offset/2;
967	int bitOffset = offset%2;
968	data += byteOffset;
969
970	/* Bits to set start at bitOffset into data */
971	*result = getBits(data, val, 2, bitOffset);
972	return 0;
973	}
974	#endif
975
976	case PT_8BSI:
977	{
978	int8_t val = data[offset];
979	*result = val;
980	return 0;
981	}
982	case PT_8BUI:
983	{
984	uint8_t val = data[offset];
985	*result = val;
986	return 0;
987	}
988	case PT_16BSI:
989	{
990	int16_t ptr = (int16_t)data; /* we assume correct alignment */
991	*result = ptr[offset];
992	return 0;
993	}
994	case PT_16BUI:
995	{
996	uint16_t ptr = (uint16_t)data; /* we assume correct alignment */
997	*result = ptr[offset];
998	return 0;
999	}
1000	case PT_32BSI:
1001	{
1002	int32_t ptr = (int32_t)data; /* we assume correct alignment */
1003	*result = ptr[offset];
1004	return 0;
1005	}
1006	case PT_32BUI:
1007	{
1008	uint32_t ptr = (uint32_t)data; /* we assume correct alignment */
1009	*result = ptr[offset];
1010	return 0;
1011	}
1012	case PT_32BF:
1013	{
1014	float ptr = (float)data; /* we assume correct alignment */
1015	*result = ptr[offset];
1016	return 0;
1017	}
1018	case PT_64BF:
1019	{
1020	double ptr = (double)data; /* we assume correct alignment */
1021	*result = ptr[offset];
1022	return 0;
1023	}
1024	default:
1025	{
1026	ctx->err("Unknown pixeltype %d", pixtype);
1027	return -1;
1028	}
1029	}
1030	}
1031
1032	double
1033	rt_band_get_nodata(rt_context ctx, rt_band band)
1034	{
1035	assert(NULL != ctx);
1036	assert(NULL != band);
1037
1038	if (!band->hasnodata)
1039	ctx->warn("Getting NODATA value for a band without NODATA values. Using %g", band->nodataval);
1040
1041	return band->nodataval;
1042	}
1043
1044
1045	/- rt_raster --------------------------------------------------------/
1046
1047	struct rt_raster_serialized_t {
1048
1049	/---[ 8 byte boundary ]---{ /
1050	uint32_t size; /* required by postgresql: 4 bytes */
1051	uint16_t version; /* format version (this is version 0): 2 bytes */
1052	uint16_t numBands; /* Number of bands: 2 bytes */
1053
1054	/* }---[ 8 byte boundary ]---{ */
1055	double scaleX; /* pixel width: 8 bytes */
1056
1057	/* }---[ 8 byte boundary ]---{ */
1058	double scaleY; /* pixel height: 8 bytes */
1059
1060	/* }---[ 8 byte boundary ]---{ */
1061	double ipX; /* insertion point X: 8 bytes */
1062
1063	/* }---[ 8 byte boundary ]---{ */
1064	double ipY; /* insertion point Y: 8 bytes */
1065
1066	/* }---[ 8 byte boundary ]---{ */
1067	double skewX; /* skew about the X axis: 8 bytes */
1068
1069	/* }---[ 8 byte boundary ]---{ */
1070	double skewY; /* skew about the Y axis: 8 bytes */
1071
1072	/* }---[ 8 byte boundary ]--- */
1073	int32_t srid; /* Spatial reference id: 4 bytes */
1074	uint16_t width; /* pixel columns: 2 bytes */
1075	uint16_t height; /* pixel rows: 2 bytes */
1076	};
1077
1078	/* NOTE: the initial part of this structure matches the layout
1079	* of data in the serialized form version 0, starting
1080	* from the numBands element
1081	*/
1082	struct rt_raster_t {
1083
1084	uint32_t size;
1085	uint16_t version;
1086
1087	/* Number of bands, all share the same dimension
1088	* and georeference */
1089	uint16_t numBands;
1090
1091	/* Georeference (in projection units) */
1092	double scaleX; /* pixel width */
1093	double scaleY; /* pixel height */
1094	double ipX; /* geo x ordinate of the corner of upper-left pixel */
1095	double ipY; /* geo y ordinate of the corner of bottom-right pixel */
1096	double skewX; /* skew about the X axis*/
1097	double skewY; /* skew about the Y axis */
1098
1099	int32_t srid; /* spatial reference id */
1100	uint16_t width; /* pixel columns - max 65535 */
1101	uint16_t height; /* pixel rows - max 65535 */
1102	rt_band bands; / actual bands */
1103
1104	};
1105
1106	rt_raster
1107	rt_raster_new(rt_context ctx, uint16_t width, uint16_t height)
1108	{
1109	rt_raster ret = NULL;
1110
1111	assert(NULL != ctx);
1112	assert(NULL != ctx->alloc);
1113
1114	ret = (rt_raster)ctx->alloc(sizeof(struct rt_raster_t));
1115	if ( ! ret ) {
1116	ctx->err("Out of virtual memory creating an rt_raster");
1117	return 0;
1118	}
1119
1120	#ifdef POSTGIS_RASTER_API_DEBUG
1121	ctx->info("Created rt_raster @ %p", ret);
1122	#endif
1123
1124	assert(NULL != ret);
1125
1126	ret->width = width;
1127
1128	ret->height = height;
1129	ret->scaleX = 1;
1130	ret->scaleY = 1;
1131	ret->ipX = 0.0;
1132	ret->ipY = 0.0;
1133	ret->skewX = 0.0;
1134	ret->skewY = 0.0;
1135	ret->srid = -1;
1136
1137	ret->numBands = 0;
1138	ret->bands = 0;
1139
1140	return ret;
1141	}
1142
1143	void
1144	rt_raster_destroy(rt_context ctx, rt_raster raster)
1145	{
1146	#ifdef POSTGIS_RASTER_API_DEBUG
1147	ctx->info("Destroying rt_raster @ %p", raster);
1148	#endif
1149	if ( raster->bands )
1150	{
1151	ctx->dealloc(raster->bands);
1152	}
1153	ctx->dealloc(raster);
1154	}
1155
1156	uint16_t
1157	rt_raster_get_width(rt_context ctx, rt_raster raster)
1158	{
1159	assert(NULL != ctx);
1160	assert(NULL != raster);
1161
1162	return raster->width;
1163	}
1164
1165	uint16_t
1166	rt_raster_get_height(rt_context ctx, rt_raster raster)
1167	{
1168	assert(NULL != ctx);
1169	assert(NULL != raster);
1170
1171	return raster->height;
1172	}
1173
1174	void
1175	rt_raster_set_pixel_sizes(rt_context ctx, rt_raster raster,
1176	double scaleX, double scaleY)
1177	{
1178	assert(NULL != ctx);
1179	assert(NULL != raster);
1180
1181	raster->scaleX = scaleX;
1182	raster->scaleY = scaleY;
1183	}
1184
1185
1186	double
1187	rt_raster_get_pixel_width(rt_context ctx, rt_raster raster)
1188	{
1189	assert(NULL != ctx);
1190	assert(NULL != raster);
1191
1192	return raster->scaleX;
1193	}
1194
1195	double
1196	rt_raster_get_pixel_height(rt_context ctx, rt_raster raster)
1197	{
1198	assert(NULL != ctx);
1199	assert(NULL != raster);
1200
1201	return raster->scaleY;
1202	}
1203
1204	void
1205	rt_raster_set_skews(rt_context ctx, rt_raster raster,
1206	double skewX, double skewY)
1207	{
1208	assert(NULL != ctx);
1209	assert(NULL != raster);
1210
1211	raster->skewX = skewX;
1212	raster->skewY = skewY;
1213	}
1214
1215	double
1216	rt_raster_get_x_skew(rt_context ctx, rt_raster raster)
1217	{
1218	assert(NULL != ctx);
1219	assert(NULL != raster);
1220
1221	return raster->skewX;
1222	}
1223
1224	double
1225	rt_raster_get_y_skew(rt_context ctx, rt_raster raster)
1226	{
1227	assert(NULL != ctx);
1228	assert(NULL != raster);
1229
1230	return raster->skewY;
1231	}
1232
1233	void
1234	rt_raster_set_offsets(rt_context ctx, rt_raster raster, double x, double y)
1235	{
1236	assert(NULL != ctx);
1237	assert(NULL != raster);
1238
1239	raster->ipX = x;
1240	raster->ipY = y;
1241	}
1242
1243	double
1244	rt_raster_get_x_offset(rt_context ctx, rt_raster raster)
1245	{
1246	assert(NULL != ctx);
1247	assert(NULL != raster);
1248
1249	return raster->ipX;
1250	}
1251
1252	double
1253	rt_raster_get_y_offset(rt_context ctx, rt_raster raster)
1254	{
1255	assert(NULL != ctx);
1256	assert(NULL != raster);
1257
1258	return raster->ipY;
1259	}
1260
1261	int32_t
1262	rt_raster_get_srid(rt_context ctx, rt_raster raster)
1263	{
1264	assert(NULL != ctx);
1265	assert(NULL != raster);
1266
1267	return raster->srid;
1268	}
1269
1270	void
1271	rt_raster_set_srid(rt_context ctx, rt_raster raster, int32_t srid)
1272	{
1273	assert(NULL != ctx);
1274	assert(NULL != raster);
1275
1276	raster->srid = srid;
1277	}
1278
1279	int
1280	rt_raster_get_num_bands(rt_context ctx, rt_raster raster)
1281	{
1282	assert(NULL != ctx);
1283	assert(NULL != raster);
1284
1285	return raster->numBands;
1286	}
1287
1288	rt_band
1289	rt_raster_get_band(rt_context ctx, rt_raster raster, int n)
1290	{
1291	assert(NULL != ctx);
1292	assert(NULL != raster);
1293
1294	if ( n >= raster->numBands ) return 0;
1295	return raster->bands[n];
1296	}
1297
1298	int32_t
1299	rt_raster_add_band(rt_context ctx, rt_raster raster, rt_band band, int index)
1300	{
1301	rt_band *oldbands = NULL;
1302	rt_band oldband = NULL;
1303	rt_band tmpband = NULL;
1304	uint16_t i = 0;
1305
1306	assert(NULL != ctx);
1307	assert(NULL != raster);
1308
1309	#ifdef POSTGIS_RASTER_API_DEBUG
1310	ctx->info("Adding band %p to raster %p", band, raster);
1311	#endif
1312
1313	if ( band->width != raster->width )
1314	{
1315	ctx->err("Can't add a %dx%d band to a %dx%d raster",
1316	band->width, band->height, raster->width, raster->height);
1317	return -1;
1318	}
1319
1320	if (index > raster->numBands)
1321	index = raster->numBands;
1322
1323	if (index < 0)
1324	index = 0;
1325
1326	oldbands = raster->bands;
1327
1328	#ifdef POSTGIS_RASTER_API_DEBUG
1329	ctx->info("Oldbands at %p", oldbands);
1330	#endif
1331
1332	raster->bands = (rt_band*)ctx->realloc(raster->bands,
1333	sizeof(rt_band)*(raster->numBands + 1)
1334	);
1335
1336	#ifdef POSTGIS_RASTER_API_DEBUG
1337	ctx->info("realloc returned %p", raster->bands);
1338	#endif
1339
1340	if ( ! raster->bands ) {
1341	ctx->err("Out of virtual memory "
1342	"reallocating band pointers");
1343	raster->bands = oldbands;
1344	return -1;
1345	}
1346
1347	for (i = 0; i <= raster->numBands; ++i)
1348	{
1349	if (i == index)
1350	{
1351	oldband = raster->bands[i];
1352	raster->bands[i] = band;
1353	}
1354	else if (i > index)
1355	{
1356	tmpband = raster->bands[i];
1357	raster->bands[i] = oldband;
1358	oldband = tmpband;
1359	}
1360	}
1361
1362	raster->numBands++;
1363	return index;
1364	}
1365
1366	void
1367	rt_raster_cell_to_geopoint(rt_context ctx, rt_raster raster,
1368	double x, double y,
1369	double* x1, double* y1)
1370	{
1371	assert(NULL != ctx);
1372	assert(NULL != raster);
1373	assert(NULL != x1);
1374	assert(NULL != y1);
1375
1376	/* Six parameters affine transformation */
1377	x1 = raster->scaleXx + raster->skewX*y + raster->ipX;
1378	y1 = raster->scaleYy + raster->skewY*x + raster->ipY;
1379
1380	#ifdef POSTGIS_RASTER_API_DEBUG
1381	ctx->info("rt_raster_cell_to_geopoint(%g,%g)", x, y);
1382	ctx->info(" ipx/y:%g/%g", raster->ipX, raster->ipY);
1383	ctx->info("cell_to_geopoint: ipX:%g, ipY:%g, %g,%g -> %g,%g",
1384	raster->ipX, raster->ipY, x, y, x1, y1);
1385	#endif
1386	}
1387
1388
1389	/* WKT string representing each polygon in WKT format acompagned by its
1390	correspoding value */
1391	struct rt_wktgeomval_t {
1392	int srid;
1393	double val;
1394	char * wktGeom;
1395	};
1396
1397
1398	rt_wktgeomval
1399	rt_raster_dump_as_wktpolygons(rt_context ctx, rt_raster raster, int nband,
1400	int * pnElements)
1401	{
1402
1403	char * pszDataPointer;
1404	char szGdalOption[50];
1405	long j;
1406	GDALDataType nPixelType = GDT_Unknown;
1407	char * apszOptions[4];
1408	OGRSFDriverH ogr_drv = NULL;
1409	GDALDriverH gdal_drv = NULL;
1410	GDALDatasetH memdataset = NULL;
1411	GDALRasterBandH gdal_band = NULL;
1412	OGRDataSourceH memdatasource = NULL;
1413	rt_pixtype pt = PT_END;
1414	rt_wktgeomval pols = NULL;
1415	OGRLayerH hLayer = NULL;
1416	OGRFeatureH hFeature = NULL;
1417	OGRGeometryH hGeom = NULL;
1418	OGRFieldDefnH hFldDfn = NULL;
1419	char * pszSrcText = NULL;
1420	int nFeatureCount = 0;
1421	rt_band band = NULL;
1422	int iPixVal = -1;
1423	int nValidPols = 0;
1424	double dValueToCompare = 0.0;
1425	int iBandHasNodataValue = FALSE;
1426	double dBandNoData = 0.0;
1427	double dEpsilon = 0.0;
1428	int nCont = 0;
1429
1430	/* Checkings */
1431	assert(NULL != ctx);
1432	assert(NULL != raster);
1433	assert(nband > 0 && nband <= rt_raster_get_num_bands(ctx, raster));
1434
1435	#ifdef POSTGIS_RASTER_API_DEBUG
1436	ctx->info("In rt_raster_dump_as_wktpolygons");
1437	#endif
1438
1439	/*******************************
1440	* Get band
1441	*******************************/
1442	band = rt_raster_get_band(ctx, raster, nband - 1);
1443	if (NULL == band)
1444	{
1445	ctx->err("Error getting band %d from raster", nband);
1446	return 0;
1447	}
1448
1449
1450	/*****************************
1451	* Register ogr mem driver
1452	*****************************/
1453	OGRRegisterAll();
1454
1455	#ifdef POSTGIS_RASTER_API_DEBUG
1456	ctx->info("rt_raster_dump_as_wktpolygons: creating OGR MEM vector");
1457	#endif
1458
1459
1460	/*****************************************************
1461	* Create an OGR in-memory vector for layers
1462	*****************************************************/
1463	ogr_drv = OGRGetDriverByName("Memory");
1464	memdatasource = OGR_Dr_CreateDataSource(ogr_drv, "", NULL);
1465
1466	if (NULL == memdatasource)
1467	{
1468	ctx->err("Couldn't create a OGR Datasource to store pols\n");
1469	return 0;
1470	}
1471
1472	/**
1473	* Can MEM driver create new layers?
1474	**/
1475	if (!OGR_DS_TestCapability(memdatasource, ODsCCreateLayer))
1476	{
1477	ctx->err("MEM driver can't create new layers, aborting\n");
1478	/* xxx jorgearevalo: what should we do now? */
1479	OGRReleaseDataSource(memdatasource);
1480	return 0;
1481	}
1482
1483	#ifdef POSTGIS_RASTER_API_DEBUG
1484	ctx->info("rt_raster_dump_as_wktpolygons: creating GDAL MEM raster");
1485	#endif
1486
1487	/****************************************************************
1488	* Create a GDAL MEM raster with one band, from rt_band object
1489	****************************************************************/
1490	GDALRegister_MEM();
1491
1492
1493	/**
1494	* First, create a Dataset with no bands using MEM driver
1495	**/
1496	gdal_drv = GDALGetDriverByName("MEM");
1497	memdataset = GDALCreate(gdal_drv, "",rt_band_get_width(ctx, band),
1498	rt_band_get_height(ctx, band), 0, GDT_Byte, NULL);
1499
1500	if (NULL == memdataset)
1501	{
1502	ctx->err("Couldn't create a GDALDataset to polygonize it\n");
1503	GDALDeregisterDriver(gdal_drv);
1504	GDALDestroyDriver(gdal_drv);
1505	OGRReleaseDataSource(memdatasource);
1506	return 0;
1507	}
1508
1509	/**
1510	* Add geotransform
1511	*/
1512	double adfGeoTransform[6] = {0.0, 1.0, 0.0, 0.0, 0.0, 1.0 };
1513	adfGeoTransform[0] = rt_raster_get_x_offset(ctx, raster);
1514	adfGeoTransform[1] = rt_raster_get_pixel_width(ctx, raster);
1515	adfGeoTransform[2] = rt_raster_get_x_skew(ctx, raster);
1516	adfGeoTransform[3] = rt_raster_get_y_offset(ctx, raster);
1517	adfGeoTransform[4] = rt_raster_get_y_skew(ctx, raster);
1518	adfGeoTransform[5] = rt_raster_get_pixel_height(ctx, raster);
1519	GDALSetGeoTransform(memdataset, adfGeoTransform);
1520
1521	#ifdef POSTGIS_RASTER_API_DEBUG
1522	ctx->info("rt_raster_dump_as_wktpolygons: Adding GDAL MEM raster band");
1523	#endif
1524
1525
1526	/**
1527	* Now, add the raster band
1528	*/
1529	pt = rt_band_get_pixtype(ctx, band);
1530
1531	switch(pt)
1532	{
1533	case PT_1BB: case PT_2BUI: case PT_4BUI: case PT_8BSI: case PT_8BUI:
1534	nPixelType = GDT_Byte;
1535	break;
1536
1537	case PT_16BSI: case PT_16BUI:
1538	if (pt == PT_16BSI)
1539	nPixelType = GDT_Int16;
1540	else
1541	nPixelType = GDT_UInt16;
1542	break;
1543
1544	case PT_32BSI: case PT_32BUI: case PT_32BF:
1545	if (pt == PT_32BSI)
1546	nPixelType = GDT_Int32;
1547	else if (pt == PT_32BUI)
1548	nPixelType = GDT_UInt32;
1549	else
1550	nPixelType = GDT_Float32;
1551	break;
1552
1553	case PT_64BF:
1554	nPixelType = GDT_Float64;
1555	break;
1556
1557	default:
1558	ctx->warn("Unknown pixel type for band\n");
1559	nPixelType = GDT_Unknown;
1560	break;
1561	}
1562
1563	void * pVoid = rt_band_get_data(ctx, band);
1564
1565	#ifdef POSTGIS_RASTER_API_DEBUG
1566	ctx->info("rt_raster_dump_as_wktpolygons: Band data is at pos %p", pVoid);
1567	#endif
1568
1569	/**
1570	* Be careful!! If this pointer is defined as szDataPointer[20]
1571	* the sprintf crash! Probable because the memory is taken from
1572	* an invalid memory context for PostgreSQL.
1573	* And be careful with size too: 10 characters may be insufficient
1574	* to store 64bits memory addresses
1575	*/
1576	pszDataPointer = (char )ctx->alloc(20 sizeof(char));
1577	sprintf(pszDataPointer, "%p", pVoid);
1578
1579	#ifdef POSTGIS_RASTER_API_DEBUG
1580	ctx->info("rt_raster_dump_as_wktpolygons: szDatapointer is %p",
1581	pszDataPointer);
1582	#endif
1583
1584	if (strnicmp(pszDataPointer,"0x", 2) == 0)
1585	sprintf(szGdalOption, "DATAPOINTER=%s", pszDataPointer);
1586	else
1587	sprintf(szGdalOption, "DATAPOINTER=0x%s", pszDataPointer);
1588
1589	#ifdef POSTGIS_RASTER_API_DEBUG
1590	ctx->info("rt_raster_dump_as_wktpolygons: Storing info for GDAL MEM raster \
1591	band");
1592	#endif
1593
1594	apszOptions[0] = szGdalOption;
1595	apszOptions[1] = NULL; // pixel offset, not needed
1596	apszOptions[2] = NULL; // line offset, not needed
1597	apszOptions[3] = NULL;
1598
1599	/**
1600	* This memory must be deallocated because we own it. The GDALRasterBand
1601	* destructor will not deallocate it
1602	**/
1603	ctx->dealloc(pszDataPointer);
1604
1605	if (GDALAddBand(memdataset, nPixelType, apszOptions) == CE_Failure)
1606	{
1607	ctx->err("Couldn't transform WKT Raster band in GDALRasterBand format to polygonize it");
1608	GDALClose(memdataset);
1609	GDALDeregisterDriver(gdal_drv);
1610	GDALDestroyDriver(gdal_drv);
1611	OGRReleaseDataSource(memdatasource);
1612
1613	return 0;
1614	}
1615
1616	/* Checking */
1617	if (GDALGetRasterCount(memdataset) != 1)
1618	{
1619	ctx->err("Error creating GDAL MEM raster bands");
1620	GDALClose(memdataset);
1621	GDALDeregisterDriver(gdal_drv);
1622	GDALDestroyDriver(gdal_drv);
1623	OGRReleaseDataSource(memdatasource);
1624	return 0;
1625	}
1626
1627
1628	#ifdef POSTGIS_RASTER_API_DEBUG
1629	ctx->info("rt_raster_dump_as_wktpolygons: polygonizying GDAL MEM raster band");
1630	#endif
1631
1632	/*****************************
1633	* Polygonize the raster band
1634	*****************************/
1635
1636	/**
1637	* From GDALPolygonize function header: "Polygon features will be
1638	* created on the output layer, with polygon geometries representing
1639	* the polygons". So,the WKB geometry type should be "wkbPolygon"
1640	**/
1641	hLayer = OGR_DS_CreateLayer(memdatasource, "Polygonized layer", NULL,
1642	wkbPolygon, NULL);
1643
1644	if (NULL == hLayer)
1645	{
1646	ctx->err("Couldn't create layer to store polygons");
1647	GDALClose(memdataset);
1648	GDALDeregisterDriver(gdal_drv);
1649	GDALDestroyDriver(gdal_drv);
1650	OGRReleaseDataSource(memdatasource);
1651	return 0;
1652	}
1653
1654	/**
1655	* Create a new field in the layer, to store the px value
1656	*/
1657
1658	/* First, create a field definition to create the field */
1659	hFldDfn = OGR_Fld_Create("Pixel value", OFTInteger);
1660
1661	/* Second, create the field */
1662	if (OGR_L_CreateField(hLayer, hFldDfn, TRUE) !=
1663	OGRERR_NONE)
1664	{
1665	ctx->warn("Couldn't create a field in OGR Layer. The polygons generated won't be able to store the pixel value");
1666	iPixVal = -1;
1667	}
1668
1669	else
1670	{
1671	/* Index to the new field created in the layer */
1672	iPixVal = 0;
1673	}
1674
1675
1676	/* Get GDAL raster band */
1677	gdal_band = GDALGetRasterBand(memdataset, 1);
1678	if (NULL == gdal_band)
1679	{
1680	ctx->err("Couldn't get GDAL band to polygonize");
1681	GDALClose(memdataset);
1682	GDALDeregisterDriver(gdal_drv);
1683	GDALDestroyDriver(gdal_drv);
1684
1685	OGR_Fld_Destroy(hFldDfn);
1686	OGR_DS_DeleteLayer(memdatasource, 0);
1687	OGRReleaseDataSource(memdatasource);
1688
1689	return 0;
1690	}
1691
1692	iBandHasNodataValue = rt_band_get_hasnodata_flag(ctx, band);
1693	if (iBandHasNodataValue)
1694	{
1695	/* Add nodata value for band */
1696	dBandNoData = rt_band_get_nodata(ctx, band);
1697	if (GDALSetRasterNoDataValue(gdal_band, dBandNoData) != CE_None)
1698	ctx->warn("Couldn't set nodata value for band.");
1699	}
1700
1701	/**
1702	* We don't need a raster mask band. Each band has a nodata value.
1703	**/
1704	GDALPolygonize(gdal_band, NULL, hLayer, iPixVal, NULL, NULL, NULL);
1705
1706	/*********************************************************************
1707	* Transform OGR layers in WKT polygons
1708	* XXX jorgearevalo: GDALPolygonize does not set the coordinate system
1709	* on the output layer. Application code should do this when the layer
1710	* is created, presumably matching the raster coordinate system.
1711	* XXX jorgearevalo: modify GDALPolygonize to directly emit polygons
1712	* in WKT format?
1713	*********************************************************************/
1714	nFeatureCount = OGR_L_GetFeatureCount(hLayer, TRUE);
1715
1716
1717	/*********************************************************************
1718	* Now, we need to:
1719	* 1.- Count the number of polygons with a field's value distinct
1720	* from nodata value for this band.
1721	* 2.- Allocate memory for this number of rt_wktgeomval structures.
1722	* 3.- Fill these structures with the needed values of the selected
1723	* polygons.
1724	*
1725	* We can do it in, at least, 2 ways:
1726	* a) 2 loops. The first one to count the number of polygons and the
1727	* second one to fill the structures using only valid polygons
1728	* b) Allocate memory for one structure. Then, in a loop, for each
1729	* valid polygon, reallocate memory for one additional structure
1730	*
1731	* I think the b) option is less efficient, because too much sys calls
1732	* and the memory fragmentation (as result of many realloactions).
1733	* I choose a), for this reason
1734	*********************************************************************/
1735
1736
1737
1738	#ifdef POSTGIS_RASTER_API_DEBUG
1739	ctx->info("rt_raster_dump_as_wktpolygons: counting valid polygons");
1740	#endif
1741
1742	/*********************************************************************
1743	* Count the "valid" polygons. This is, the polygons with the "Pixel
1744	* Value" field value distinct from raster nodata value
1745	*********************************************************************/
1746	nValidPols = nFeatureCount;
1747	if (iBandHasNodataValue)
1748	{
1749	dBandNoData = GDALGetRasterNoDataValue(gdal_band, NULL);
1750	for(j = 0; j < nFeatureCount; j++)
1751	{
1752	hFeature = OGR_L_GetFeature(hLayer, j);
1753
1754
1755	/**
1756	* The field was stored as int, but we can use this function
1757	* because uses "atof" to transform the string representation of
1758	* the number into a double. We shouldn't have problems
1759	**/
1760	dValueToCompare = OGR_F_GetFieldAsDouble(hFeature, iPixVal);
1761
1762
1763	/**
1764	* Compare value obtained and nodata from band. If aren't equal,
1765	* the associated polygon is discarded
1766	**/
1767	dEpsilon = fabs(dValueToCompare - dBandNoData);
1768	if (dEpsilon <= FLT_EPSILON)
1769	nValidPols--;
1770
1771	OGR_F_Destroy(hFeature);
1772	}
1773	}
1774
1775
1776
1777	/* Only allocate for valid pols!!! */
1778
1779	pols = (rt_wktgeomval)ctx->alloc(nValidPols *
1780	sizeof(struct rt_wktgeomval_t));
1781
1782	if (NULL == pols)
1783	{
1784	ctx->err("Couldn't allocate memory for wktgeomval structure");
1785	GDALClose(memdataset);
1786	GDALDeregisterDriver(gdal_drv);
1787	GDALDestroyDriver(gdal_drv);
1788
1789	OGR_Fld_Destroy(hFldDfn);
1790	OGR_DS_DeleteLayer(memdatasource, 0);
1791	OGRReleaseDataSource(memdatasource);
1792
1793	return 0;
1794	}
1795
1796
1797	#ifdef POSTGIS_RASTER_API_DEBUG
1798	ctx->info("rt_raster_dump_as_wktpolygons: storing polygons (%d)",
1799	nFeatureCount);
1800	#endif
1801
1802	nCont = 0;
1803	if (pnElements)
1804	*pnElements = 0;
1805
1806	for(j = 0; j < nFeatureCount; j++)
1807	{
1808
1809	hFeature = OGR_L_GetFeature(hLayer, j);
1810	dValueToCompare = OGR_F_GetFieldAsDouble(hFeature, iPixVal);
1811	dEpsilon = fabs(dValueToCompare - dBandNoData);
1812
1813	if (dEpsilon > FLT_EPSILON \|\| !iBandHasNodataValue)
1814	{
1815	hGeom = OGR_F_GetGeometryRef(hFeature);
1816	OGR_G_ExportToWkt(hGeom, &pszSrcText);
1817
1818	pols[nCont].val = dValueToCompare;
1819	pols[nCont].srid = rt_raster_get_srid(ctx, raster);
1820	pols[nCont].wktGeom = (char *)ctx->alloc((1 + strlen(pszSrcText))
1821	* sizeof(char));
1822	strcpy(pols[nCont].wktGeom, pszSrcText);
1823
1824	#ifdef POSTGIS_RASTER_API_DEBUG
1825	ctx->info("Feature %d, WKT Polygon: %s", j, pols[nCont].wktGeom);
1826	ctx->info("Feature %d, value: %d", j, (int)(pols[nCont].val));
1827	ctx->info("Feature %d, srid: %d", j, pols[nCont].srid);
1828	#endif
1829	nCont++;
1830
1831	/**
1832	* We can't use deallocator from rt_context, because it comes from
1833	* postgresql backend, that uses pfree. This function needs a
1834	* postgresql memory context to work with, and the memory created
1835	* for pszSrcText is created outside this context.
1836	**/
1837	//ctx->dealloc(pszSrcText);
1838	free(pszSrcText);
1839	pszSrcText = NULL;
1840	}
1841
1842	OGR_F_Destroy(hFeature);
1843	}
1844
1845	if (pnElements)
1846	*pnElements = nCont;
1847
1848	#ifdef POSTGIS_RASTER_API_DEBUG
1849	ctx->info("rt_raster_dump_as_wktpolygons: destroying GDAL MEM raster");
1850	#endif
1851
1852	GDALClose(memdataset);
1853	GDALDeregisterDriver(gdal_drv);
1854	GDALDestroyDriver(gdal_drv);
1855
1856
1857	#ifdef POSTGIS_RASTER_API_DEBUG
1858	ctx->info("rt_raster_dump_as_wktpolygons: destroying OGR MEM vector");
1859	#endif
1860	OGR_Fld_Destroy(hFldDfn);
1861	OGR_DS_DeleteLayer(memdatasource, 0);
1862	OGRReleaseDataSource(memdatasource);
1863
1864	return pols;
1865	}
1866
1867
1868	LWPOLY*
1869	rt_raster_get_convex_hull(rt_context ctx, rt_raster raster)
1870	{
1871	POINTARRAY **rings = NULL;
1872	POINTARRAY *pts = NULL;
1873	LWPOLY* ret = NULL;
1874	POINT4D p4d;
1875
1876	assert(NULL != ctx);
1877	assert(NULL != raster);
1878
1879	#ifdef POSTGIS_RASTER_API_DEBUG
1880	ctx->info("rt_raster_get_convex_hull: raster is %dx%d",
1881	raster->width, raster->height);
1882	#endif
1883	if ( (!raster->width) \|\| (!raster->height) )
1884	{
1885	return 0;
1886	}
1887
1888	rings = (POINTARRAY *)ctx->alloc(sizeof(POINTARRAY));
1889	if ( ! rings )
1890	{
1891	ctx->err("Out of memory [%s:%d]", __FILE__, __LINE__);
1892	return 0;
1893	}
1894	rings[0] = ptarray_construct(0, 0, 5);
1895	/* TODO: handle error on ptarray construction */
1896	/* XXX jorgearevalo: the error conditions aren't managed in ptarray_construct */
1897	if ( ! rings[0] )
1898	{
1899	ctx->err("Out of memory [%s:%d]", __FILE__, __LINE__);
1900	return 0;
1901	}
1902	pts = rings[0];
1903
1904	/* Upper-left corner (first and last points) */
1905	rt_raster_cell_to_geopoint(ctx, raster,
1906	0, 0,
1907	&p4d.x, &p4d.y);
1908	setPoint4d(pts, 0, &p4d);
1909	setPoint4d(pts, 4, &p4d); /* needed for closing it? */
1910
1911	/* Upper-right corner (we go clockwise) */
1912	rt_raster_cell_to_geopoint(ctx, raster,
1913	raster->width, 0,
1914	&p4d.x, &p4d.y);
1915	setPoint4d(pts, 1, &p4d);
1916
1917	/* Lower-right corner */
1918	rt_raster_cell_to_geopoint(ctx, raster,
1919	raster->width, raster->height,
1920	&p4d.x, &p4d.y);
1921	setPoint4d(pts, 2, &p4d);
1922
1923	/* Lower-left corner */
1924	rt_raster_cell_to_geopoint(ctx, raster,
1925	0, raster->height,
1926	&p4d.x, &p4d.y);
1927	setPoint4d(pts, 3, &p4d);
1928
1929	ret = lwpoly_construct(raster->srid, 0, 1, rings);
1930
1931	return ret;
1932	}
1933
1934	/--------- WKB I/O ---------------------------------------------------/
1935
1936	static uint8_t
1937	isMachineLittleEndian(void)
1938	{
1939	static int endian_check_int = 1; /* dont modify this!!! */
1940	/* 0=big endian\|xdr -- 1=little endian\|ndr */
1941	return ((uint8_t ) &endian_check_int);
1942	}
1943
1944	static uint8_t
1945	read_uint8(const uint8_t** from)
1946	{
1947	assert(NULL != from);
1948
1949	return (from)++;
1950	}
1951
1952	/* unused up to now
1953	static void
1954	write_uint8(uint8_t** from, uint8_t v)
1955	{
1956	assert(NULL != from);
1957
1958	(from)++ = v;
1959	}
1960	*/
1961
1962	static int8_t
1963	read_int8(const uint8_t** from)
1964	{
1965	assert(NULL != from);
1966
1967	return (int8_t)read_uint8(from);
1968	}
1969
1970	/* unused up to now
1971	static void
1972	write_int8(uint8_t** from, int8_t v)
1973	{
1974	assert(NULL != from);
1975
1976	(from)++ = v;
1977	}
1978	*/
1979
1980	static uint16_t
1981	read_uint16(const uint8_t** from, uint8_t littleEndian)
1982	{
1983	uint16_t ret = 0;
1984
1985	assert(NULL != from);
1986
1987	if ( littleEndian )
1988	{
1989	ret = (*from)[0] \|
1990	(*from)[1] << 8;
1991	}
1992	else
1993	{
1994	/* big endian */
1995	ret = (*from)[0] << 8 \|
1996	(*from)[1];
1997	}
1998	*from += 2;
1999	return ret;
2000	}
2001
2002	static void
2003	write_uint16(uint8_t** to, uint8_t littleEndian, uint16_t v)
2004	{
2005	assert(NULL != to);
2006
2007	if ( littleEndian )
2008	{
2009	(*to)[0] = v & 0x00FF;
2010	(*to)[1] = v >> 8;
2011	}
2012	else
2013	{
2014	(*to)[1] = v & 0x00FF;
2015	(*to)[0] = v >> 8;
2016	}
2017	*to += 2;
2018	}
2019
2020	static int16_t
2021	read_int16(const uint8_t** from, uint8_t littleEndian)
2022	{
2023	assert(NULL != from);
2024
2025	return read_uint16(from, littleEndian);
2026	}
2027
2028	/* unused up to now
2029	static void
2030	write_int16(uint8_t** to, uint8_t littleEndian, int16_t v)
2031	{
2032	assert(NULL != to);
2033
2034	if ( littleEndian )
2035	{
2036	(*to)[0] = v & 0x00FF;
2037	(*to)[1] = v >> 8;
2038	}
2039	else
2040	{
2041	(*to)[1] = v & 0x00FF;
2042	(*to)[0] = v >> 8;
2043	}
2044	*to += 2;
2045	}
2046	*/
2047
2048	static uint32_t
2049	read_uint32(const uint8_t** from, uint8_t littleEndian)
2050	{
2051	uint32_t ret = 0;
2052
2053	assert(NULL != from);
2054
2055	if ( littleEndian )
2056	{
2057	ret = (uint32_t) ((*from)[0] & 0xff) \|
2058	(uint32_t) ((*from)[1] & 0xff) << 8 \|
2059	(uint32_t) ((*from)[2] & 0xff) << 16 \|
2060	(uint32_t) ((*from)[3] & 0xff) << 24;
2061	}
2062	else
2063	{
2064	/* big endian */
2065	ret = (uint32_t) ((*from)[3] & 0xff) \|
2066	(uint32_t) ((*from)[2] & 0xff) << 8 \|
2067	(uint32_t) ((*from)[1] & 0xff) << 16 \|
2068	(uint32_t) ((*from)[0] & 0xff) << 24;
2069	}
2070
2071	*from += 4;
2072	return ret;
2073	}
2074
2075	/* unused up to now
2076	static void
2077	write_uint32(uint8_t** to, uint8_t littleEndian, uint32_t v)
2078	{
2079	assert(NULL != to);
2080
2081	if ( littleEndian )
2082	{
2083	(*to)[0] = v & 0x000000FF;
2084	(*to)[1] = ( v & 0x0000FF00 ) >> 8;
2085	(*to)[2] = ( v & 0x00FF0000 ) >> 16;
2086	(*to)[3] = ( v & 0xFF000000 ) >> 24;
2087	}
2088	else
2089	{
2090	(*to)[3] = v & 0x000000FF;
2091	(*to)[2] = ( v & 0x0000FF00 ) >> 8;
2092	(*to)[1] = ( v & 0x00FF0000 ) >> 16;
2093	(*to)[0] = ( v & 0xFF000000 ) >> 24;
2094	}
2095	*to += 4;
2096	}
2097	*/
2098
2099	static int32_t
2100	read_int32(const uint8_t** from, uint8_t littleEndian)
2101	{
2102	assert(NULL != from);
2103
2104	return read_uint32(from, littleEndian);
2105	}
2106
2107	/* unused up to now
2108	static void
2109	write_int32(uint8_t** to, uint8_t littleEndian, int32_t v)
2110	{
2111	assert(NULL != to);
2112
2113	if ( littleEndian )
2114	{
2115	(*to)[0] = v & 0x000000FF;
2116	(*to)[1] = ( v & 0x0000FF00 ) >> 8;
2117	(*to)[2] = ( v & 0x00FF0000 ) >> 16;
2118	(*to)[3] = ( v & 0xFF000000 ) >> 24;
2119	}
2120	else
2121	{
2122	(*to)[3] = v & 0x000000FF;
2123	(*to)[2] = ( v & 0x0000FF00 ) >> 8;
2124	(*to)[1] = ( v & 0x00FF0000 ) >> 16;
2125	(*to)[0] = ( v & 0xFF000000 ) >> 24;
2126	}
2127	*to += 4;
2128	}
2129	*/
2130
2131	static float
2132	read_float32(const uint8_t** from, uint8_t littleEndian)
2133	{
2134	union {
2135	float f;
2136	uint32_t i;
2137	} ret;
2138
2139	ret.i = read_uint32(from, littleEndian);
2140
2141	return ret.f;
2142	}
2143
2144	/* unused up to now
2145	static void
2146	write_float32(uint8_t** from, uint8_t littleEndian, float f)
2147	{
2148	union {
2149	float f;
2150	uint32_t i;
2151	} u;
2152
2153	u.f = f;
2154	write_uint32(from, littleEndian, u.i);
2155	}
2156	*/
2157
2158	static double
2159	read_float64(const uint8_t** from, uint8_t littleEndian)
2160	{
2161	union {
2162	double d;
2163	uint64_t i;
2164	} ret;
2165
2166	assert(NULL != from);
2167
2168	if ( littleEndian )
2169	{
2170	ret.i = (uint64_t) ((*from)[0] & 0xff) \|
2171	(uint64_t) ((*from)[1] & 0xff) << 8 \|
2172	(uint64_t) ((*from)[2] & 0xff) << 16 \|
2173	(uint64_t) ((*from)[3] & 0xff) << 24 \|
2174	(uint64_t) ((*from)[4] & 0xff) << 32 \|
2175	(uint64_t) ((*from)[5] & 0xff) << 40 \|
2176	(uint64_t) ((*from)[6] & 0xff) << 48 \|
2177	(uint64_t) ((*from)[7] & 0xff) << 56;
2178	}
2179	else
2180	{
2181	/* big endian */
2182	ret.i = (uint64_t) ((*from)[7] & 0xff) \|
2183	(uint64_t) ((*from)[6] & 0xff) << 8 \|
2184	(uint64_t) ((*from)[5] & 0xff) << 16 \|
2185	(uint64_t) ((*from)[4] & 0xff) << 24 \|
2186	(uint64_t) ((*from)[3] & 0xff) << 32 \|
2187	(uint64_t) ((*from)[2] & 0xff) << 40 \|
2188	(uint64_t) ((*from)[1] & 0xff) << 48 \|
2189	(uint64_t) ((*from)[0] & 0xff) << 56;
2190	}
2191
2192	*from += 8;
2193	return ret.d;
2194	}
2195
2196	/* unused up to now
2197	static void
2198	write_float64(uint8_t** to, uint8_t littleEndian, double v)
2199	{
2200	union {
2201	double d;
2202	uint64_t i;
2203	} u;
2204
2205	assert(NULL != to);
2206
2207	u.d = v;
2208
2209	if ( littleEndian )
2210	{
2211	(*to)[0] = u.i & 0x00000000000000FFULL;
2212	(*to)[1] = ( u.i & 0x000000000000FF00ULL ) >> 8;
2213	(*to)[2] = ( u.i & 0x0000000000FF0000ULL ) >> 16;
2214	(*to)[3] = ( u.i & 0x00000000FF000000ULL ) >> 24;
2215	(*to)[4] = ( u.i & 0x000000FF00000000ULL ) >> 32;
2216	(*to)[5] = ( u.i & 0x0000FF0000000000ULL ) >> 40;
2217	(*to)[6] = ( u.i & 0x00FF000000000000ULL ) >> 48;
2218	(*to)[7] = ( u.i & 0xFF00000000000000ULL ) >> 56;
2219	}
2220	else
2221	{
2222	(*to)[7] = u.i & 0x00000000000000FFULL;
2223	(*to)[6] = ( u.i & 0x000000000000FF00ULL ) >> 8;
2224	(*to)[5] = ( u.i & 0x0000000000FF0000ULL ) >> 16;
2225	(*to)[4] = ( u.i & 0x00000000FF000000ULL ) >> 24;
2226	(*to)[3] = ( u.i & 0x000000FF00000000ULL ) >> 32;
2227	(*to)[2] = ( u.i & 0x0000FF0000000000ULL ) >> 40;
2228	(*to)[1] = ( u.i & 0x00FF000000000000ULL ) >> 48;
2229	(*to)[0] = ( u.i & 0xFF00000000000000ULL ) >> 56;
2230	}
2231	*to += 8;
2232	}
2233	*/
2234
2235	#define BANDTYPE_FLAGS_MASK 0xF0
2236	#define BANDTYPE_PIXTYPE_MASK 0x0F
2237	#define BANDTYPE_FLAG_OFFDB (1<<7)
2238	#define BANDTYPE_FLAG_HASNODATA (1<<6)
2239	#define BANDTYPE_FLAG_RESERVED2 (1<<5)
2240	#define BANDTYPE_FLAG_RESERVED3 (1<<4)
2241
2242	#define BANDTYPE_PIXTYPE(x) ((x)&BANDTYPE_PIXTYPE_MASK)
2243	#define BANDTYPE_IS_OFFDB(x) ((x)&BANDTYPE_FLAG_OFFDB)
2244	#define BANDTYPE_HAS_NODATA(x) ((x)&BANDTYPE_FLAG_HASNODATA)
2245
2246	/* Read band from WKB as at start of band */
2247	static rt_band
2248	rt_band_from_wkb(rt_context ctx, uint16_t width, uint16_t height,
2249	const uint8_t** ptr, const uint8_t* end,
2250	uint8_t littleEndian)
2251	{
2252	rt_band band = NULL;
2253	int pixbytes = 0;
2254	uint8_t type = 0;
2255	unsigned long sz = 0;
2256	uint32_t v = 0;
2257
2258	assert(NULL != ctx);
2259	assert(NULL != ptr);
2260	assert(NULL != end);
2261
2262	band = ctx->alloc(sizeof(struct rt_band_t));
2263	if ( ! band )
2264	{
2265	ctx->err("Out of memory allocating rt_band during WKB parsing");
2266	return 0;
2267	}
2268
2269	if ( end - *ptr < 1 ) {
2270	ctx->err("Premature end of WKB on band reading (%s:%d)",
2271	__FILE__, __LINE__);
2272	return 0;
2273	}
2274	type = read_uint8(ptr);
2275
2276	if ( (type & BANDTYPE_PIXTYPE_MASK) >= PT_END )
2277	{
2278	ctx->err("Invalid pixtype %d", type & BANDTYPE_PIXTYPE_MASK);
2279	ctx->dealloc(band);
2280	return 0;
2281	}
2282	assert(NULL != band);
2283
2284	band->pixtype = type & BANDTYPE_PIXTYPE_MASK;
2285	band->offline = BANDTYPE_IS_OFFDB(type) ? 1 : 0;
2286	band->hasnodata = BANDTYPE_HAS_NODATA(type) ? 1 : 0;
2287	band->width = width;
2288	band->height = height;
2289
2290	#ifdef POSTGIS_RASTER_API_DEBUG
2291	ctx->info(" Band pixtype:%s, offline:%d, hasnodata:%d",
2292	rt_pixtype_name(ctx, band->pixtype),
2293	band->offline,
2294	band->hasnodata);
2295	#endif
2296
2297	/* Check there's enough bytes to read nodata value */
2298
2299	pixbytes = rt_pixtype_size(ctx, band->pixtype);
2300	if (((*ptr)+pixbytes) >= end) {
2301	ctx->err("Premature end of WKB on band novalue reading");
2302	ctx->dealloc(band);
2303	return 0;
2304	}
2305
2306	/* Read nodata value */
2307	switch (band->pixtype)
2308	{
2309	case PT_1BB:
2310	{
2311	band->nodataval = ((int)read_uint8(ptr)) & 0x01;
2312	break;
2313	}
2314	case PT_2BUI:
2315	{
2316	band->nodataval = ((int)read_uint8(ptr)) & 0x03;
2317	break;
2318	}
2319	case PT_4BUI:
2320	{
2321	band->nodataval = ((int)read_uint8(ptr)) & 0x0F;
2322	break;
2323	}
2324	case PT_8BSI:
2325	{
2326	band->nodataval = read_int8(ptr);
2327	break;
2328	}
2329	case PT_8BUI:
2330	{
2331	band->nodataval = read_uint8(ptr);
2332	break;
2333	}
2334	case PT_16BSI:
2335	{
2336	band->nodataval = read_int16(ptr, littleEndian);
2337	break;
2338	}
2339	case PT_16BUI:
2340	{
2341	band->nodataval = read_uint16(ptr, littleEndian);
2342	break;
2343	}
2344	case PT_32BSI:
2345	{
2346	band->nodataval = read_int32(ptr, littleEndian);
2347	break;
2348	}
2349	case PT_32BUI:
2350	{
2351	band->nodataval = read_uint32(ptr, littleEndian);
2352	break;
2353	}
2354	case PT_32BF:
2355	{
2356	band->nodataval = read_float32(ptr, littleEndian);
2357	break;
2358	}
2359	case PT_64BF:
2360	{
2361	band->nodataval = read_float64(ptr, littleEndian);
2362	break;
2363	}
2364	default:
2365	{
2366	ctx->err("Unknown pixeltype %d", band->pixtype);
2367	ctx->dealloc(band);
2368	return 0;
2369	}
2370	}
2371
2372	#ifdef POSTGIS_RASTER_API_DEBUG
2373	ctx->info(" Nodata value: %g, pixbytes: %d, ptr @ %p, end @ %p",
2374	band->nodataval, pixbytes, *ptr, end);
2375	#endif
2376
2377	if ( band->offline )
2378	{
2379	if (((*ptr)+1) >= end) {
2380	ctx->err("Premature end of WKB on offline "
2381	"band data bandNum reading (%s:%d)",
2382	__FILE__, __LINE__);
2383	ctx->dealloc(band);
2384	return 0;
2385	}
2386	band->data.offline.bandNum = read_int8(ptr);
2387
2388	{
2389	/* check we have a NULL-termination */
2390	sz = 0;
2391	while ( (ptr)[sz] && &((ptr)[sz]) < end) ++sz;
2392	if ( &((*ptr)[sz]) >= end )
2393	{
2394	ctx->err("Premature end of WKB on band offline path reading");
2395	ctx->dealloc(band);
2396	return 0;
2397	}
2398
2399	band->ownsData = 1;
2400	band->data.offline.path = ctx->alloc(sz + 1);
2401
2402	memcpy(band->data.offline.path, *ptr, sz);
2403	band->data.offline.path[sz] = '\0';
2404
2405	#ifdef POSTGIS_RASTER_API_DEBUG
2406	ctx->info("OFFDB band path is %s (size is %d)",
2407	band->data.offline.path, sz);
2408	#endif
2409	*ptr += sz + 1;
2410	}
2411	return band;
2412	}
2413
2414	/* This is an on-disk band */
2415	sz = width * height * pixbytes;
2416	if (((*ptr)+sz) > end)
2417	{
2418	ctx->err("Premature end of WKB on band data reading (%s:%d)",
2419	__FILE__, __LINE__);
2420	ctx->dealloc(band);
2421	return 0;
2422	}
2423
2424	band->data.mem = ctx->alloc(sz);
2425	if ( ! band->data.mem )
2426	{
2427	ctx->err("Out of memory during band creation in WKB parser");
2428	ctx->dealloc(band);
2429	return 0;
2430	}
2431
2432	memcpy(band->data.mem, *ptr, sz);
2433	*ptr += sz;
2434
2435	/* Should now flip values if > 8bit and
2436	* littleEndian != isMachineLittleEndian */
2437	if ( pixbytes > 1 )
2438	{
2439	if ( isMachineLittleEndian() != littleEndian )
2440	{{
2441	void (flipper)(uint8_t) = 0;
2442	uint8_t *flipme = NULL;
2443
2444	if ( pixbytes == 2 ) flipper = flip_endian_16;
2445	else if ( pixbytes == 4 ) flipper = flip_endian_32;
2446	else if ( pixbytes == 8 ) flipper = flip_endian_64;
2447	else {
2448	ctx->err("Unexpected pix bytes %d", pixbytes);
2449	ctx->dealloc(band);
2450	ctx->dealloc(band->data.mem);
2451	return 0;
2452	}
2453
2454	flipme = band->data.mem;
2455	sz = width * height;
2456	for (v=0; v<sz; ++v)
2457	{
2458	flipper(flipme);
2459	flipme += pixbytes;
2460	}
2461	}}
2462	}
2463
2464	/* And should check for invalid values for < 8bit types */
2465	else if ( band->pixtype == PT_1BB \|\|
2466	band->pixtype == PT_2BUI \|\|
2467	band->pixtype == PT_4BUI )
2468	{{
2469	uint8_t maxVal = band->pixtype == PT_1BB ? 1 :
2470	band->pixtype == PT_2BUI ? 3 :
2471	15 ;
2472	uint8_t val;
2473
2474	sz = width*height;
2475	for (v=0; v<sz; ++v)
2476	{
2477	val = ((uint8_t*)band->data.mem)[v];
2478	if ( val > maxVal )
2479	{
2480	ctx->err("Invalid value %d for pixel of type %s",
2481	val, rt_pixtype_name(ctx, band->pixtype));
2482	ctx->dealloc(band->data.mem);
2483	ctx->dealloc(band);
2484	return 0;
2485	}
2486	}
2487	}}
2488
2489	return band;
2490	}
2491
2492	/* -4 for size, +1 for endian */
2493	#define RT_WKB_HDR_SZ (sizeof(struct rt_raster_serialized_t)-4+1)
2494
2495	rt_raster
2496	rt_raster_from_wkb(rt_context ctx, const uint8_t* wkb, uint32_t wkbsize)
2497	{
2498	const uint8_t *ptr = wkb;
2499	const uint8_t *wkbend = NULL;
2500	rt_raster rast = NULL;
2501	uint8_t endian = 0;
2502	uint16_t version = 0;
2503	uint16_t i = 0;
2504
2505	assert(NULL != ctx);
2506	assert(NULL != ptr);
2507
2508	/* Check that wkbsize is >= sizeof(rt_raster_serialized) */
2509	if ( wkbsize < RT_WKB_HDR_SZ )
2510	{
2511	ctx->err("rt_raster_from_wkb: wkb size < min size (%d)",
2512	RT_WKB_HDR_SZ);
2513	return 0;
2514	}
2515	wkbend = wkb + wkbsize;
2516
2517	#ifdef POSTGIS_RASTER_API_DEBUG
2518	ctx->info("Parsing header from wkb position %d (expected 0)",
2519	d_binptr_to_pos(ptr, wkbend, wkbsize));
2520	#endif
2521
2522	CHECK_BINPTR_POSITION(ptr, wkbend, wkbsize, 0);
2523
2524	/* Read endianness */
2525	endian = *ptr;
2526	ptr += 1;
2527
2528	/* Read version of protocol */
2529	version = read_uint16(&ptr, endian);
2530	if ( version != 0 )
2531	{
2532	ctx->err("rt_raster_from_wkb: WKB version %d unsupported", version);
2533	return 0;
2534	}
2535
2536	/* Read other components of raster header */
2537	rast = (rt_raster)ctx->alloc(sizeof(struct rt_raster_t));
2538	if ( ! rast )
2539	{
2540	ctx->err("Out of memory allocating raster for wkb input");
2541	return 0;
2542	}
2543	rast->numBands = read_uint16(&ptr, endian);
2544	rast->scaleX = read_float64(&ptr, endian);
2545	rast->scaleY = read_float64(&ptr, endian);
2546	rast->ipX = read_float64(&ptr, endian);
2547	rast->ipY = read_float64(&ptr, endian);
2548	rast->skewX = read_float64(&ptr, endian);
2549	rast->skewY = read_float64(&ptr, endian);
2550	rast->srid = read_int32(&ptr, endian);
2551	rast->width = read_uint16(&ptr, endian);
2552	rast->height = read_uint16(&ptr, endian);
2553
2554	/* Consistency checking, should have been checked before */
2555	assert(ptr <= wkbend);
2556
2557	#ifdef POSTGIS_RASTER_API_DEBUG
2558	ctx->info("rt_raster_from_wkb: Raster numBands: %d",
2559	rast->numBands);
2560	ctx->info("rt_raster_from_wkb: Raster scale: %gx%g",
2561	rast->scaleX, rast->scaleY);
2562	ctx->info("rt_raster_from_wkb: Raster ip: %gx%g",
2563	rast->ipX, rast->ipY);
2564	ctx->info("rt_raster_from_wkb: Raster skew: %gx%g",
2565	rast->skewX, rast->skewY );
2566	ctx->info("rt_raster_from_wkb: Raster srid: %d",
2567	rast->srid);
2568	ctx->info("rt_raster_from_wkb: Raster dims: %dx%d",
2569	rast->width, rast->height);
2570	ctx->info("Parsing raster header finished at wkb position %d (expected 61)",
2571	d_binptr_to_pos(ptr, wkbend, wkbsize));
2572	#endif
2573
2574	CHECK_BINPTR_POSITION(ptr, wkbend, wkbsize, 61);
2575
2576	/* Read all bands of raster */
2577
2578	if ( ! rast->numBands )
2579	{
2580	/* Here ptr should have been left to right after last used byte */
2581	if (ptr < wkbend)
2582	{
2583	ctx->info("%d bytes of WKB remained unparsed", wkbend-ptr);
2584	}
2585	else if ( ptr > wkbend )
2586	{
2587	/* Easier to get a segfault before I guess */
2588	ctx->warn("We parsed %d bytes more then available!", ptr-wkbend);
2589	}
2590	rast->bands = 0;
2591	return rast;
2592	}
2593
2594	/* Now read the bands */
2595	rast->bands = (rt_band)ctx->alloc(sizeof(rt_band) rast->numBands);
2596	if ( ! rast->bands )
2597	{
2598	ctx->err("Out of memory allocating bands for WKB raster decoding");
2599	ctx->dealloc(rast);
2600	return 0;
2601	}
2602
2603	/* ptr should now point to start of first band */
2604	assert ( ptr <= wkbend ); /* we should have checked this before */
2605
2606	for (i = 0; i < rast->numBands; ++i)
2607	{
2608	#ifdef POSTGIS_RASTER_API_DEBUG
2609	ctx->info("Parsing band %d from wkb position %d", i,
2610	d_binptr_to_pos(ptr, wkbend, wkbsize));
2611	#endif
2612	rt_band band = rt_band_from_wkb(ctx, rast->width, rast->height,
2613	&ptr, wkbend, endian);
2614	if ( ! band )
2615	{
2616	ctx->err("Error reading WKB form of band %d", i);
2617	ctx->dealloc(rast);
2618	/* TODO: dealloc any previously allocated band too ! */
2619	return 0;
2620	}
2621	rast->bands[i] = band;
2622	}
2623
2624	/* Here ptr should have been left to right after last used byte */
2625	if (ptr < wkbend)
2626	{
2627	ctx->info("%d bytes of WKB remained unparsed", wkbend-ptr);
2628	}
2629	else if ( ptr > wkbend )
2630	{
2631	/* Easier to get a segfault before I guess */
2632	ctx->warn("We parsed %d bytes more then available!",
2633	ptr-wkbend);
2634	}
2635
2636	assert(NULL != rast);
2637	return rast;
2638	}
2639
2640	rt_raster
2641	rt_raster_from_hexwkb(rt_context ctx, const char* hexwkb,
2642	uint32_t hexwkbsize)
2643	{
2644	uint8_t* wkb = NULL;
2645	uint32_t wkbsize = 0;
2646	uint32_t i = 0;
2647
2648	assert(NULL != ctx);
2649	assert(NULL != hexwkb);
2650
2651	#ifdef POSTGIS_RASTER_API_DEBUG
2652	ctx->info("rt_raster_from_hexwkb: input wkb: %s", hexwkb);
2653	#endif
2654
2655	if ( hexwkbsize % 2 )
2656	{
2657	ctx->err("Raster HEXWKB input must have an even number of characters");
2658	return 0;
2659	}
2660	wkbsize = hexwkbsize / 2;
2661
2662	wkb = ctx->alloc(wkbsize);
2663	if ( ! wkb )
2664	{
2665	ctx->err("Out of memory allocating memory for decoding HEXWKB");
2666	return 0;
2667	}
2668
2669	for (i=0; i<wkbsize; ++i) /* parse full hex */
2670	{
2671	wkb[i] = parse_hex((char) & (hexwkb[i 2]));
2672	}
2673
2674	rt_raster ret = rt_raster_from_wkb(ctx, wkb, wkbsize);
2675
2676	ctx->dealloc(wkb); /* as long as rt_raster_from_wkb copies memory */
2677
2678	return ret;
2679	}
2680
2681	static uint32_t
2682	rt_raster_wkb_size(rt_context ctx, rt_raster raster)
2683	{
2684	uint32_t size = RT_WKB_HDR_SZ;
2685	uint16_t i = 0;
2686
2687	assert(NULL != ctx);
2688	assert(NULL != raster);
2689
2690	#ifdef POSTGIS_RASTER_API_DEBUG
2691	ctx->info("rt_raster_wkb_size: computing size for %d bands",
2692	raster->numBands);
2693	#endif
2694
2695	for (i = 0; i < raster->numBands; ++i)
2696	{
2697	rt_band band = raster->bands[i];
2698	rt_pixtype pixtype = band->pixtype;
2699	int pixbytes = rt_pixtype_size(ctx, pixtype);
2700
2701	#ifdef POSTGIS_RASTER_API_DEBUG
2702	ctx->info("rt_raster_wkb_size: adding size of band %d", i);
2703	#endif
2704
2705	if ( pixbytes < 1 ) {
2706	ctx->err("Corrupted band: unkonwn pixtype");
2707	return 0;
2708	}
2709
2710	/* Add space for band type */
2711	size += 1;
2712
2713	/* Add space for nodata value */
2714	size += pixbytes;
2715
2716	if ( band->offline )
2717	{
2718	/* Add space for band number */
2719	size += 1;
2720
2721	/* Add space for null-terminated path */
2722	size += strlen(band->data.offline.path)+1;
2723	}
2724	else
2725	{
2726	/* Add space for actual data */
2727	size += pixbytesraster->widthraster->height;
2728	}
2729
2730	}
2731
2732	return size;
2733	}
2734
2735	uint8_t *
2736	rt_raster_to_wkb(rt_context ctx, rt_raster raster, uint32_t *wkbsize)
2737	{
2738	#ifdef POSTGIS_RASTER_API_DEBUG
2739	const uint8_t *wkbend = NULL;
2740	#endif
2741	uint8_t *wkb = NULL;
2742	uint8_t *ptr = NULL;
2743	uint16_t i = 0;
2744	uint8_t littleEndian = isMachineLittleEndian();
2745
2746	assert(NULL != ctx);
2747	assert(NULL != raster);
2748	assert(NULL != wkbsize);
2749
2750	#ifdef POSTGIS_RASTER_API_DEBUG
2751	ctx->info("rt_raster_to_wkb: about to call rt_raster_wkb_size");
2752	#endif
2753
2754	*wkbsize = rt_raster_wkb_size(ctx, raster);
2755
2756	#ifdef POSTGIS_RASTER_API_DEBUG
2757	ctx->info("rt_raster_to_wkb: found size: %d", *wkbsize);
2758	#endif
2759
2760	wkb = (uint8_t)ctx->alloc(wkbsize);
2761	if ( ! wkb )
2762	{
2763	ctx->err("Out of memory allocating WKB for raster");
2764	return 0;
2765	}
2766
2767	ptr = wkb;
2768
2769	#ifdef POSTGIS_RASTER_API_DEBUG
2770	wkbend = ptr + (*wkbsize);
2771	ctx->info("Writing raster header to wkb on position %d (expected 0)",
2772	d_binptr_to_pos(ptr, wkbend, *wkbsize));
2773	#endif
2774
2775	/* Write endianness */
2776	*ptr = littleEndian;
2777	ptr += 1;
2778
2779	/* Write version(size - (end - ptr)) */
2780	write_uint16(&ptr, littleEndian, 0);
2781
2782	/* Copy header (from numBands up) */
2783	memcpy(ptr, &(raster->numBands), sizeof(struct rt_raster_serialized_t) - 6);
2784	ptr += sizeof(struct rt_raster_serialized_t) - 6;
2785
2786	#ifdef POSTGIS_RASTER_API_DEBUG
2787	ctx->info("Writing bands header to wkb position %d (expected 61)",
2788	d_binptr_to_pos(ptr, wkbend, *wkbsize));
2789	#endif
2790
2791	/* Serialize bands now */
2792	for (i = 0; i < raster->numBands; ++i)
2793	{
2794	rt_band band = raster->bands[i];
2795	rt_pixtype pixtype = band->pixtype;
2796	int pixbytes = rt_pixtype_size(ctx, pixtype);
2797
2798	#ifdef POSTGIS_RASTER_API_DEBUG
2799	ctx->info("Writing WKB for band %d", i);
2800	ctx->info("Writing band pixel type to wkb position %d",
2801	d_binptr_to_pos(ptr, wkbend, *wkbsize));
2802	#endif
2803
2804	if ( pixbytes < 1 ) {
2805	ctx->err("Corrupted band: unkonwn pixtype");
2806	return 0;
2807	}
2808
2809	/* Add band type */
2810	*ptr = band->pixtype;
2811	if ( band->offline ) *ptr \|= BANDTYPE_FLAG_OFFDB;
2812	if ( band->hasnodata ) *ptr \|= BANDTYPE_FLAG_HASNODATA;
2813	ptr += 1;
2814
2815	#if 0 // no padding required for WKB
2816	/* Add padding (if needed) */
2817	if ( pixbytes > 1 ) {
2818	memset(ptr, '\0', pixbytes-1);
2819	ptr += pixbytes-1;
2820	}
2821	/* Consistency checking (ptr is pixbytes-aligned) */
2822	assert(! (((uint64_t)ptr)%pixbytes) );
2823	#endif
2824
2825	#ifdef POSTGIS_RASTER_API_DEBUG
2826	ctx->info("Writing band nodata to wkb position %d",
2827	d_binptr_to_pos(ptr, wkbend, *wkbsize));
2828	#endif
2829
2830	/* Add nodata value */
2831	switch (pixtype)
2832	{
2833	case PT_1BB:
2834	case PT_2BUI:
2835	case PT_4BUI:
2836	case PT_8BUI:
2837	{
2838	uint8_t v = band->nodataval;
2839	*ptr = v; ptr+=1;
2840	break;
2841	}
2842	case PT_8BSI:
2843	{
2844	int8_t v = band->nodataval;
2845	*ptr = v; ptr+=1;
2846	break;
2847	}
2848	case PT_16BSI:
2849	case PT_16BUI:
2850	{
2851	uint16_t v = band->nodataval;
2852	memcpy(ptr, &v, 2); ptr+=2;
2853	break;
2854	}
2855	case PT_32BSI:
2856	case PT_32BUI:
2857	{
2858	uint32_t v = band->nodataval;
2859	memcpy(ptr, &v, 4); ptr+=4;
2860	break;
2861	}
2862	case PT_32BF:
2863	{
2864	float v = band->nodataval;
2865	memcpy(ptr, &v, 4); ptr+=4;
2866	break;
2867	}
2868	case PT_64BF:
2869	{
2870	memcpy(ptr, &band->nodataval, 8); ptr+=8;
2871	break;
2872	}
2873	default:
2874	ctx->err("Fatal error caused by unknown pixel type. Aborting.");
2875	abort(); /* shoudn't happen */
2876	return 0;
2877	}
2878
2879	#if 0 // no padding for WKB
2880	/* Consistency checking (ptr is pixbytes-aligned) */
2881	assert(! ((uint64_t)ptr%pixbytes) );
2882	#endif
2883
2884	if ( band->offline )
2885	{
2886	/* Write band number */
2887	*ptr = band->data.offline.bandNum; ptr += 1;
2888
2889	/* Write path */
2890	strcpy((char*)ptr, band->data.offline.path);
2891	ptr += strlen(band->data.offline.path)+1;
2892	}
2893	else
2894	{
2895	/* Write data */
2896	{
2897	uint32_t datasize = raster->width * raster->height * pixbytes;
2898	memcpy(ptr, band->data.mem, datasize);
2899	ptr += datasize;
2900	}
2901	}
2902
2903	#if 0 // no padding for WKB
2904	/* Pad up to 8-bytes boundary */
2905	while ( (uint64_t)ptr%8 ) {
2906	*ptr = 0; ++ptr;
2907	}
2908
2909	/* Consistency checking (ptr is pixbytes-aligned) */
2910	assert(! ((uint64_t)ptr%pixbytes) );
2911	#endif
2912	}
2913
2914	return wkb;
2915	}
2916
2917	char *
2918	rt_raster_to_hexwkb(rt_context ctx, rt_raster raster, uint32_t *hexwkbsize)
2919	{
2920	uint8_t *wkb = NULL;
2921	char* hexwkb = NULL;
2922	uint32_t i = 0;
2923	uint32_t wkbsize = 0;
2924
2925	assert(NULL != ctx);
2926	assert(NULL != raster);
2927	assert(NULL != hexwkbsize);
2928
2929	#ifdef POSTGIS_RASTER_API_DEBUG
2930	ctx->info("rt_raster_to_hexwkb: calling rt_raster_to_wkb");
2931	#endif
2932
2933	wkb = rt_raster_to_wkb(ctx, raster, &wkbsize);
2934
2935	#ifdef POSTGIS_RASTER_API_DEBUG
2936	ctx->info("rt_raster_to_hexwkb: rt_raster_to_wkb returned");
2937	#endif
2938
2939	hexwkbsize = wkbsize 2; /* hex is 2 times bytes */
2940	hexwkb = (char)ctx->alloc( (hexwkbsize)+1);
2941	if ( ! hexwkb )
2942	{
2943	ctx->dealloc(wkb);
2944	ctx->err("Out of memory hexifying raster WKB");
2945	return 0;
2946	}
2947	hexwkb[hexwkbsize] = '\0'; / Null-terminate */
2948
2949	for (i = 0; i < wkbsize; ++i)
2950	{
2951	deparse_hex(wkb[i], &(hexwkb[2 * i]));
2952	}
2953
2954	ctx->dealloc(wkb); /* we don't need this anymore */
2955
2956	#ifdef POSTGIS_RASTER_API_DEBUG
2957	ctx->info("rt_raster_to_hexwkb: output wkb: %s", hexwkb);
2958	#endif
2959
2960	return hexwkb;
2961	}
2962
2963
2964	/--------- Serializer/Deserializer --------------------------------------/
2965
2966	static uint32_t
2967	rt_raster_serialized_size(rt_context ctx, rt_raster raster)
2968	{
2969	uint32_t size = sizeof(struct rt_raster_serialized_t);
2970	uint16_t i = 0;
2971
2972	assert(NULL != ctx);
2973	assert(NULL != raster);
2974
2975	#ifdef POSTGIS_RASTER_API_DEBUG
2976	ctx->info("Serialized size with just header:%d - now adding size of %d bands",
2977	size, raster->numBands);
2978	#endif
2979
2980	for (i = 0; i < raster->numBands; ++i)
2981	{
2982	rt_band band = raster->bands[i];
2983	rt_pixtype pixtype = band->pixtype;
2984	int pixbytes = rt_pixtype_size(ctx, pixtype);
2985
2986	if ( pixbytes < 1 ) {
2987	ctx->err("Corrupted band: unknown pixtype");
2988	return 0;
2989	}
2990
2991	/* Add space for band type, hasnodata flag and data padding */
2992	size += pixbytes;
2993
2994	/* Add space for nodata value */
2995	size += pixbytes;
2996
2997	if ( band->offline )
2998	{
2999	/* Add space for band number */
3000	size += 1;
3001
3002	/* Add space for null-terminated path */
3003	size += strlen(band->data.offline.path)+1;
3004	}
3005	else
3006	{
3007	/* Add space for raster band data */
3008	size += pixbytes * raster->width * raster->height;
3009	}
3010
3011	#ifdef POSTGIS_RASTER_API_DEBUG
3012	ctx->info("Size before alignment is %d", size);
3013	#endif
3014
3015	/* Align size to 8-bytes boundary (trailing padding) */
3016	size += 8 - (size % 8);
3017
3018	#ifdef POSTGIS_RASTER_API_DEBUG
3019	ctx->info("Size after alignment is %d", size);
3020	#endif
3021	}
3022
3023	return size;
3024	}
3025
3026	void*
3027	rt_raster_serialize(rt_context ctx, rt_raster raster)
3028	{
3029	uint32_t size = rt_raster_serialized_size(ctx, raster);
3030	uint8_t* ret = NULL;
3031	uint8_t* ptr = NULL;
3032	uint16_t i = 0;
3033
3034	assert(NULL != ctx);
3035	assert(NULL != ctx->alloc);
3036	assert(NULL != raster);
3037
3038	ret = (uint8_t*)ctx->alloc(size);
3039	if ( ! ret )
3040	{
3041	ctx->err("Out of memory allocating %d bytes for serializing a raster",
3042	size);
3043	return 0;
3044	}
3045	memset(ret, '-', size);
3046
3047	ptr = ret;
3048
3049	#ifdef POSTGIS_RASTER_API_DEBUG
3050	ctx->info("sizeof(struct rt_raster_serialized_t):%u",
3051	sizeof(struct rt_raster_serialized_t));
3052	ctx->info("sizeof(struct rt_raster_t):%u",
3053	sizeof(struct rt_raster_t));
3054	ctx->info("serialized size:%lu", (long unsigned)size);
3055	#endif
3056
3057	/* Set size */
3058	/* NOTE: Value of rt_raster.size may be updated in
3059	* returned object, for instance, by rt_pg layer to
3060	* store value calculated by SET_VARSIZE.
3061	*/
3062	raster->size = size;
3063
3064	/* Set version */
3065	raster->version = 0;
3066
3067	/* Copy header */
3068	memcpy(ptr, raster, sizeof(struct rt_raster_serialized_t));
3069
3070	#ifdef POSTGIS_RASTER_API_DEBUG
3071	ctx->info("Start hex dump of raster being serialized using 0x2D to mark non-written bytes");
3072	uint8_t* dbg_ptr = ptr;
3073	d_print_binary_hex(ctx, "HEADER", dbg_ptr, size);
3074	#endif
3075
3076	ptr += sizeof(struct rt_raster_serialized_t);
3077
3078	/* Serialize bands now */
3079	for (i = 0; i < raster->numBands; ++i)
3080	{
3081	rt_band band = raster->bands[i];
3082	assert(NULL != band);
3083
3084	rt_pixtype pixtype = band->pixtype;
3085	int pixbytes = rt_pixtype_size(ctx, pixtype);
3086	if ( pixbytes < 1 )
3087	{
3088	ctx->err("Corrupted band: unkonwn pixtype");
3089	return 0;
3090	}
3091
3092	/* Add band type */
3093	*ptr = band->pixtype;
3094	if ( band->offline ) { *ptr \|= BANDTYPE_FLAG_OFFDB; }
3095	if ( band->hasnodata ) { *ptr \|= BANDTYPE_FLAG_HASNODATA; }
3096
3097	#ifdef POSTGIS_RASTER_API_DEBUG
3098	d_print_binary_hex(ctx, "PIXTYPE", dbg_ptr, size);
3099	#endif
3100
3101	ptr += 1;
3102
3103	/* Add padding (if needed) */
3104	if ( pixbytes > 1 ) {
3105	memset(ptr, '\0', pixbytes - 1);
3106	ptr += pixbytes - 1;
3107	}
3108
3109	#ifdef POSTGIS_RASTER_API_DEBUG
3110	d_print_binary_hex(ctx, "PADDING", dbg_ptr, size);
3111	#endif
3112
3113	/* Consistency checking (ptr is pixbytes-aligned) */
3114	assert(! (((uintptr_t)ptr) % pixbytes) );
3115
3116	/* Add nodata value */
3117	switch (pixtype)
3118	{
3119	case PT_1BB:
3120	case PT_2BUI:
3121	case PT_4BUI:
3122	case PT_8BUI:
3123	{
3124	uint8_t v = band->nodataval;
3125	*ptr = v; ptr += 1;
3126	break;
3127	}
3128	case PT_8BSI:
3129	{
3130	int8_t v = band->nodataval;
3131	*ptr = v; ptr += 1;
3132	break;
3133	}
3134	case PT_16BSI:
3135	case PT_16BUI:
3136	{
3137	uint16_t v = band->nodataval;
3138	memcpy(ptr, &v, 2); ptr += 2;
3139	break;
3140	}
3141	case PT_32BSI:
3142	case PT_32BUI:
3143	{
3144	uint32_t v = band->nodataval;
3145	memcpy(ptr, &v, 4); ptr += 4;
3146	break;
3147	}
3148	case PT_32BF:
3149	{
3150	float v = band->nodataval;
3151	memcpy(ptr, &v, 4); ptr += 4;
3152	break;
3153	}
3154	case PT_64BF:
3155	{
3156	memcpy(ptr, &band->nodataval, 8); ptr += 8;
3157	break;
3158	}
3159	default:
3160	ctx->err("Fatal error caused by unknown pixel type. Aborting.");
3161	abort(); /* shoudn't happen */
3162	return 0;
3163	}
3164
3165	/* Consistency checking (ptr is pixbytes-aligned) */
3166	assert(! ((uintptr_t)ptr % pixbytes) );
3167
3168	#ifdef POSTGIS_RASTER_API_DEBUG
3169	d_print_binary_hex(ctx, "NODATA", dbg_ptr, size);
3170	#endif
3171
3172	if ( band->offline )
3173	{
3174	/* Write band number */
3175	*ptr = band->data.offline.bandNum;
3176	ptr += 1;
3177
3178	/* Write path */
3179	strcpy((char*)ptr, band->data.offline.path);
3180	ptr += strlen(band->data.offline.path) + 1;
3181	}
3182	else
3183	{
3184	/* Write data */
3185	uint32_t datasize = raster->width * raster->height * pixbytes;
3186	memcpy(ptr, band->data.mem, datasize);
3187	ptr += datasize;
3188	}
3189
3190	#ifdef POSTGIS_RASTER_API_DEBUG
3191	d_print_binary_hex(ctx, "BAND", dbg_ptr, size);
3192	#endif
3193
3194	/* Pad up to 8-bytes boundary */
3195	while ((uintptr_t)ptr % 8)
3196	{
3197	*ptr = 0;
3198	++ptr;
3199
3200	#ifdef POSTGIS_RASTER_API_DEBUG
3201	ctx->info("PAD at %d", (uint64_t)ptr % 8);
3202	#endif
3203	}
3204
3205	/* Consistency checking (ptr is pixbytes-aligned) */
3206	assert(! ((uintptr_t)ptr % pixbytes) );
3207
3208	} /* for-loop over bands */
3209
3210	#ifdef POSTGIS_RASTER_API_DEBUG
3211	d_print_binary_hex(ctx, "SERIALIZED RASTER", dbg_ptr, size);
3212	#endif
3213
3214	return ret;
3215	}
3216
3217	rt_raster
3218	rt_raster_deserialize(rt_context ctx, void* serialized)
3219	{
3220	rt_raster rast = NULL;
3221	const uint8_t *ptr = NULL;
3222	const uint8_t *beg = NULL;
3223	uint16_t i = 0;
3224	uint8_t littleEndian = isMachineLittleEndian();
3225
3226	assert(NULL != ctx);
3227	assert(NULL != serialized);
3228
3229	/* NOTE: Value of rt_raster.size may be different
3230	* than actual size of raster data being read.
3231	* See note on SET_VARSIZE in rt_raster_serialize function above.
3232	*/
3233
3234	/* Allocate memory for deserialized raster header */
3235	rast = (rt_raster)ctx->alloc(sizeof(struct rt_raster_t));
3236	if ( ! rast )
3237	{
3238	ctx->err("Out of memory allocating raster for deserialization");
3239	return 0;
3240	}
3241
3242	/* Deserialize raster header */
3243	memcpy(rast, serialized, sizeof(struct rt_raster_serialized_t));
3244
3245	if ( ! rast->numBands )
3246	{
3247	rast->bands = 0;
3248	return rast;
3249	}
3250	beg = (const uint8_t*)serialized;
3251
3252	/* Allocate registry of raster bands */
3253	rast->bands = ctx->alloc(rast->numBands * sizeof(rt_band));
3254
3255	#ifdef POSTGIS_RASTER_API_DEBUG
3256	ctx->info("rt_raster_deserialize: %d bands", rast->numBands);
3257	#endif
3258
3259	/* Move to the beginning of first band */
3260	ptr = beg;
3261	ptr += sizeof(struct rt_raster_serialized_t);
3262
3263	/* Deserialize bands now */
3264	for (i = 0; i < rast->numBands; ++i)
3265	{
3266	rt_band band = NULL;
3267	uint8_t type = 0;
3268	int pixbytes = 0;
3269
3270	band = ctx->alloc(sizeof(struct rt_band_t));
3271	if ( ! band )
3272	{
3273	ctx->err("Out of memory allocating rt_band during deserialization");
3274	return 0;
3275	}
3276
3277	rast->bands[i] = band;
3278
3279	type = *ptr;
3280	ptr++;
3281	band->pixtype = type & BANDTYPE_PIXTYPE_MASK;
3282
3283	#ifdef POSTGIS_RASTER_API_DEBUG
3284	ctx->info("rt_raster_deserialize: band %d with pixel type %s", i, rt_pixtype_name(ctx, band->pixtype));
3285	#endif
3286
3287	band->offline = BANDTYPE_IS_OFFDB(type) ? 1 : 0;
3288	band->hasnodata = BANDTYPE_HAS_NODATA(type) ? 1: 0;
3289	band->width = rast->width;
3290	band->height = rast->height;
3291	band->ownsData = 0;
3292
3293	/* Advance by data padding */
3294	pixbytes = rt_pixtype_size(ctx, band->pixtype);
3295	ptr += pixbytes - 1;
3296
3297	/* Read nodata value */
3298	switch (band->pixtype)
3299	{
3300	case PT_1BB:
3301	{
3302	band->nodataval = ((int)read_uint8(&ptr)) & 0x01;
3303	break;
3304	}
3305	case PT_2BUI:
3306	{
3307	band->nodataval = ((int)read_uint8(&ptr)) & 0x03;
3308	break;
3309	}
3310	case PT_4BUI:
3311	{
3312	band->nodataval = ((int)read_uint8(&ptr)) & 0x0F;
3313	break;
3314	}
3315	case PT_8BSI:
3316	{
3317	band->nodataval = read_int8(&ptr);
3318	break;
3319	}
3320	case PT_8BUI:
3321	{
3322	band->nodataval = read_uint8(&ptr);
3323	break;
3324	}
3325	case PT_16BSI:
3326	{
3327	band->nodataval = read_int16(&ptr, littleEndian);
3328	break;
3329	}
3330	case PT_16BUI:
3331	{
3332	band->nodataval = read_uint16(&ptr, littleEndian);
3333	break;
3334	}
3335	case PT_32BSI:
3336	{
3337	band->nodataval = read_int32(&ptr, littleEndian);
3338	break;
3339	}
3340	case PT_32BUI:
3341	{
3342	band->nodataval = read_uint32(&ptr, littleEndian);
3343	break;
3344	}
3345	case PT_32BF:
3346	{
3347	band->nodataval = read_float32(&ptr, littleEndian);
3348	break;
3349	}
3350	case PT_64BF:
3351	{
3352	band->nodataval = read_float64(&ptr, littleEndian);
3353	break;
3354	}
3355	default:
3356	{
3357	ctx->err("Unknown pixeltype %d", band->pixtype);
3358	ctx->dealloc(band);
3359	ctx->dealloc(rast);
3360	return 0;
3361	}
3362	}
3363
3364	#ifdef POSTGIS_RASTER_API_DEBUG
3365	ctx->info("rt_raster_deserialize: has NODATA flag %d", band->hasnodata);
3366	ctx->info("rt_raster_deserialize: NODATA value %g", band->nodataval);
3367	#endif
3368	/* Consistency checking (ptr is pixbytes-aligned) */
3369	assert(! (((uintptr_t)ptr) % pixbytes) );
3370
3371	if ( band->offline )
3372	{
3373	/* Read band number */
3374	band->data.offline.bandNum = *ptr;
3375	ptr += 1;
3376
3377	/* Register path */
3378	band->data.offline.path = (char*)ptr;
3379	ptr += strlen(band->data.offline.path) + 1;
3380	}
3381	else
3382	{
3383	/* Register data */
3384	const uint32_t datasize = rast->width * rast->height * pixbytes;
3385	band->data.mem = (uint8_t*)ptr;
3386	ptr += datasize;
3387	}
3388
3389	/* Skip bytes of padding up to 8-bytes boundary */
3390	#ifdef POSTGIS_RASTER_API_DEBUG
3391	const uint8_t *padbeg = ptr;
3392	#endif
3393	while ( 0 != ((ptr - beg) % 8) )
3394
3395	{
3396	++ptr;
3397	}
3398
3399	#ifdef POSTGIS_RASTER_API_DEBUG
3400	ctx->info("rt_raster_deserialize: skip %d bytes of 8-bytes boundary padding", ptr - padbeg);
3401	#endif
3402	/* Consistency checking (ptr is pixbytes-aligned) */
3403	assert(! ((uintptr_t)ptr % pixbytes) );
3404	}
3405
3406	return rast;
3407	}

Note: See TracBrowser for help on using the browser.

Context Navigation

root/spike/wktraster/rt_core/rt_api.c

Download in other formats: