1 : /******************************************************************************
2 : * Project: MG4 Lidar GDAL Plugin
3 : * Purpose: Provide an orthographic view of MG4-encoded Lidar dataset
4 : * for use with LizardTech Lidar SDK version 1.1.0
5 : * Author: Michael Rosen, mrosen@lizardtech.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2010, LizardTech
9 : * All rights reserved.
10 :
11 : * Redistribution and use in source and binary forms, with or without
12 : * modification, are permitted provided that the following conditions are met:
13 : *
14 : * Redistributions of source code must retain the above copyright notice,
15 : * this list of conditions and the following disclaimer.
16 : *
17 : * Redistributions in binary form must reproduce the above copyright notice,
18 : * this list of conditions and the following disclaimer in the documentation
19 : * and/or other materials provided with the distribution.
20 : *
21 : * Neither the name of the LizardTech nor the names of its contributors may
22 : * be used to endorse or promote products derived from this software without
23 : * specific prior written permission.
24 : *
25 : * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
26 : * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 : * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
28 : * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
29 : * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
30 : * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
31 : * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
32 : * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
33 : * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
34 : * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
35 : * POSSIBILITY OF SUCH DAMAGE.
36 : ****************************************************************************/
37 : #include "lidar/MG4PointReader.h"
38 : #include "lidar/FileIO.h"
39 : #include "lidar/Error.h"
40 : #include "lidar/Version.h"
41 : #include <float.h>
42 : LT_USE_LIDAR_NAMESPACE
43 :
44 : #include "gdal_pam.h"
45 : // #include "gdal_alg.h" // 1.6 and later have gridding algorithms
46 :
47 : CPL_C_START
48 : //void __declspec(dllexport) GDALRegister_MG4Lidar(void);
49 : void CPL_DLL GDALRegister_MG4Lidar(void);
50 : CPL_C_END
51 :
52 : /************************************************************************/
53 : /* ==================================================================== */
54 : /* MG4LidarDataset */
55 : /* ==================================================================== */
56 : /************************************************************************/
57 :
58 : // Resolution Ratio between adjacent levels.
59 : #define RESOLUTION_RATIO 2.0
60 :
61 : class MG4LidarRasterBand;
62 : class CropableMG4PointReader : public MG4PointReader
63 : {
64 : CONCRETE_OBJECT(CropableMG4PointReader);
65 1 : void init (IO *io, Bounds *bounds)
66 : {
67 1 : MG4PointReader::init(io);
68 1 : if (bounds != NULL)
69 1 : setBounds(*bounds);
70 1 : }
71 : };
72 1 : CropableMG4PointReader::CropableMG4PointReader() : MG4PointReader() {};
73 1 : CropableMG4PointReader::~CropableMG4PointReader() {};
74 :
75 1 : IMPLEMENT_OBJECT_CREATE(CropableMG4PointReader);
76 :
77 : static double MaxRasterSize = 2048.0;
78 : static double MaxBlockSideSize = 1024.0;
79 : class MG4LidarDataset : public GDALPamDataset
80 : {
81 : friend class MG4LidarRasterBand;
82 : public:
83 : MG4LidarDataset();
84 : ~MG4LidarDataset();
85 : static GDALDataset *Open( GDALOpenInfo * );
86 : CPLErr GetGeoTransform( double * padfTransform );
87 : const char *GetProjectionRef();
88 : protected:
89 : MG4PointReader *reader;
90 : FileIO *fileIO;
91 : CPLErr OpenZoomLevel( int Zoom );
92 : PointInfo requiredChannels;
93 : int nOverviewCount;
94 : MG4LidarDataset **papoOverviewDS;
95 : CPLXMLNode *poXMLPCView;
96 : bool ownsXML;
97 : int nBlockXSize, nBlockYSize;
98 : int iLevel;
99 : };
100 :
101 : /* ======================================== */
102 : /* MG4LidarRasterBand */
103 : /* ======================================== */
104 :
105 : class MG4LidarRasterBand : public GDALPamRasterBand
106 : {
107 : friend class MG4LidarDataset;
108 :
109 : public:
110 :
111 : MG4LidarRasterBand( MG4LidarDataset *, int, CPLXMLNode *, const char * );
112 : ~MG4LidarRasterBand();
113 :
114 : virtual CPLErr GetStatistics( int bApproxOK, int bForce, double *pdfMin, double *pdfMax, double *pdfMean, double *padfStdDev );
115 : virtual int GetOverviewCount();
116 : virtual GDALRasterBand * GetOverview( int i );
117 : virtual CPLErr IReadBlock( int, int, void * );
118 : virtual double GetNoDataValue( int *pbSuccess = NULL );
119 :
120 : protected:
121 : double getMaxValue();
122 : double nodatavalue;
123 : virtual const bool ElementPassesFilter(const PointData &, size_t);
124 : template<typename DTYPE>
125 : CPLErr doReadBlock(int, int, void *);
126 : CPLXMLNode *poxmlBand;
127 : char **papszFilterClassCodes;
128 : char **papszFilterReturnNums;
129 : const char * Aggregation;
130 : CPLString ChannelName;
131 : };
132 :
133 :
134 : /************************************************************************/
135 : /* MG4LidarRasterBand() */
136 : /************************************************************************/
137 :
138 2 : MG4LidarRasterBand::MG4LidarRasterBand( MG4LidarDataset *pods, int nband, CPLXMLNode *xmlBand, const char * name )
139 : {
140 2 : this->poDS = pods;
141 2 : this->nBand = nband;
142 2 : this->poxmlBand = xmlBand;
143 2 : this->ChannelName = name;
144 2 : this->Aggregation = NULL;
145 2 : nBlockXSize = pods->nBlockXSize;
146 2 : nBlockYSize = pods->nBlockYSize;
147 :
148 2 : switch(pods->reader->getChannel(name)->getDataType())
149 : {
150 : #define DO_CASE(ltdt, gdt) case (ltdt): eDataType = gdt; break;
151 2 : DO_CASE(DATATYPE_FLOAT64, GDT_Float64);
152 0 : DO_CASE(DATATYPE_FLOAT32, GDT_Float32);
153 0 : DO_CASE(DATATYPE_SINT32, GDT_Int32);
154 0 : DO_CASE(DATATYPE_UINT32, GDT_UInt32);
155 0 : DO_CASE(DATATYPE_SINT16, GDT_Int16);
156 0 : DO_CASE(DATATYPE_UINT16, GDT_UInt16);
157 0 : DO_CASE(DATATYPE_SINT8, GDT_Byte);
158 0 : DO_CASE(DATATYPE_UINT8, GDT_Byte);
159 : default:
160 : CPLError(CE_Failure, CPLE_AssertionFailed,
161 0 : "Invalid datatype in MG4 file");
162 : break;
163 : #undef DO_CASE
164 : }
165 : // Coerce datatypes as required.
166 2 : const char * ForceDataType = CPLGetXMLValue(pods->poXMLPCView, "Datatype", NULL);
167 :
168 2 : if (ForceDataType != NULL)
169 : {
170 2 : GDALDataType dt = GDALGetDataTypeByName(ForceDataType);
171 2 : if (dt != GDT_Unknown)
172 2 : eDataType = dt;
173 : }
174 :
175 2 : CPLXMLNode *poxmlFilter = CPLGetXMLNode(poxmlBand, "ClassificationFilter");
176 2 : if( poxmlFilter == NULL )
177 0 : poxmlFilter = CPLGetXMLNode(pods->poXMLPCView, "ClassificationFilter");
178 2 : if (poxmlFilter == NULL || poxmlFilter->psChild == NULL ||
179 : poxmlFilter->psChild->pszValue == NULL)
180 0 : papszFilterClassCodes = NULL;
181 : else
182 2 : papszFilterClassCodes = CSLTokenizeString(poxmlFilter->psChild->pszValue);
183 :
184 2 : poxmlFilter = CPLGetXMLNode(poxmlBand, "ReturnNumberFilter");
185 2 : if( poxmlFilter == NULL )
186 2 : poxmlFilter = CPLGetXMLNode(pods->poXMLPCView, "ReturnNumberFilter");
187 4 : if (poxmlFilter == NULL || poxmlFilter->psChild == NULL ||
188 : poxmlFilter->psChild->pszValue == NULL)
189 2 : papszFilterReturnNums = NULL;
190 : else
191 0 : papszFilterReturnNums = CSLTokenizeString(poxmlFilter->psChild->pszValue);
192 :
193 :
194 2 : CPLXMLNode * poxmlAggregation = CPLGetXMLNode(poxmlBand, "AggregationMethod");
195 2 : if( poxmlAggregation == NULL )
196 2 : poxmlAggregation = CPLGetXMLNode(pods->poXMLPCView, "AggregationMethod");
197 4 : if (poxmlAggregation == NULL || poxmlAggregation->psChild == NULL ||
198 : poxmlAggregation->psChild->pszValue == NULL)
199 2 : Aggregation = "Mean";
200 : else
201 0 : Aggregation = poxmlAggregation->psChild->pszValue;
202 :
203 2 : nodatavalue = getMaxValue();
204 :
205 2 : CPLXMLNode * poxmlIntepolation = CPLGetXMLNode(poxmlBand, "InterpolationMethod");
206 2 : if( poxmlIntepolation == NULL )
207 2 : poxmlIntepolation = CPLGetXMLNode(pods->poXMLPCView, "InterpolationMethod");
208 2 : if (poxmlIntepolation != NULL )
209 : {
210 0 : CPLXMLNode * poxmlMethod= NULL;
211 0 : char ** papszParams = NULL;
212 0 : if (((poxmlMethod = CPLSearchXMLNode(poxmlIntepolation, "None")) != NULL) &&
213 : poxmlMethod->psChild != NULL && poxmlMethod->psChild->pszValue != NULL)
214 : {
215 0 : papszParams = CSLTokenizeString(poxmlMethod->psChild->pszValue);
216 0 : if (!EQUAL(papszParams[0], "MAX"))
217 0 : nodatavalue = atof(papszParams[0]);
218 : }
219 : // else if .... Add support for other interpolation methods here.
220 0 : CSLDestroy(papszParams);
221 : }
222 2 : const char * filter = NULL;
223 2 : if (papszFilterClassCodes != NULL && papszFilterReturnNums != NULL)
224 0 : filter ="Classification and Return";
225 2 : if (papszFilterClassCodes != NULL)
226 2 : filter = "Classification";
227 0 : else if (papszFilterReturnNums != NULL)
228 0 : filter = "Return";
229 2 : CPLString osDesc;
230 2 : if (filter)
231 2 : osDesc.Printf("%s of %s (filtered by %s)", Aggregation, ChannelName.c_str(), filter);
232 : else
233 0 : osDesc.Printf("%s of %s", Aggregation, ChannelName.c_str());
234 2 : SetDescription(osDesc);
235 2 : }
236 :
237 : /************************************************************************/
238 : /* ~MG4LidarRasterBand() */
239 : /************************************************************************/
240 2 : MG4LidarRasterBand::~MG4LidarRasterBand()
241 : {
242 2 : CSLDestroy(papszFilterClassCodes);
243 2 : CSLDestroy(papszFilterReturnNums);
244 2 : }
245 :
246 : /************************************************************************/
247 : /* GetOverviewCount() */
248 : /************************************************************************/
249 :
250 0 : int MG4LidarRasterBand::GetOverviewCount()
251 : {
252 0 : MG4LidarDataset *poGDS = (MG4LidarDataset *) poDS;
253 0 : return poGDS->nOverviewCount;
254 : }
255 :
256 : /************************************************************************/
257 : /* GetOverview() */
258 : /************************************************************************/
259 1 : GDALRasterBand *MG4LidarRasterBand::GetOverview( int i )
260 : {
261 1 : MG4LidarDataset *poGDS = (MG4LidarDataset *) poDS;
262 :
263 1 : if( i < 0 || i >= poGDS->nOverviewCount )
264 0 : return NULL;
265 : else
266 1 : return poGDS->papoOverviewDS[i]->GetRasterBand( nBand );
267 : }
268 : template<typename DTYPE>
269 407623 : const DTYPE GetChannelElement(const ChannelData &channel, size_t idx)
270 : {
271 407623 : DTYPE retval = static_cast<DTYPE>(0);
272 407623 : switch (channel.getDataType())
273 : {
274 : case (DATATYPE_FLOAT64):
275 112133 : retval = static_cast<DTYPE>(static_cast<const double*>(channel.getData())[idx]);
276 112133 : break;
277 : case (DATATYPE_FLOAT32):
278 0 : retval = static_cast<DTYPE>(static_cast<const float *>(channel.getData())[idx]);
279 0 : break;
280 : case (DATATYPE_SINT32):
281 0 : retval = static_cast<DTYPE>(static_cast<const long *>(channel.getData())[idx]);
282 0 : break;
283 : case (DATATYPE_UINT32):
284 0 : retval = static_cast<DTYPE>(static_cast<const unsigned long *>(channel.getData())[idx]);
285 0 : break;
286 : case (DATATYPE_SINT16):
287 0 : retval = static_cast<DTYPE>(static_cast<const short *>(channel.getData())[idx]);
288 0 : break;
289 : case (DATATYPE_UINT16):
290 0 : retval = static_cast<DTYPE>(static_cast<const unsigned short *>(channel.getData())[idx]);
291 0 : break;
292 : case (DATATYPE_SINT8):
293 0 : retval = static_cast<DTYPE>(static_cast<const char *>(channel.getData())[idx]);
294 0 : break;
295 : case (DATATYPE_UINT8):
296 295490 : retval = static_cast<DTYPE>(static_cast<const unsigned char *>(channel.getData())[idx]);
297 295490 : break;
298 : case (DATATYPE_SINT64):
299 0 : retval = static_cast<DTYPE>(static_cast<const GIntBig *>(channel.getData())[idx]);
300 0 : break;
301 : case (DATATYPE_UINT64):
302 0 : retval = static_cast<DTYPE>(static_cast<const GUIntBig *>(channel.getData())[idx]);
303 : break;
304 : default:
305 : break;
306 : }
307 407623 : return retval;
308 : }
309 :
310 :
311 295490 : const bool MG4LidarRasterBand::ElementPassesFilter(const PointData &pointdata, size_t i)
312 : {
313 295490 : bool bClassificationOK = true;
314 295490 : bool bReturnNumOK = true;
315 :
316 : // Check if classification code is ok: it was requested and it does match one of the requested codes
317 295490 : const int classcode = GetChannelElement<int>(*pointdata.getChannel(CHANNEL_NAME_ClassId), i);
318 : char bufCode[16];
319 295490 : sprintf(bufCode, "%d", classcode);
320 : bClassificationOK = (papszFilterClassCodes == NULL ? true :
321 295490 : (CSLFindString(papszFilterClassCodes,bufCode)!=-1));
322 :
323 295490 : if (bClassificationOK)
324 : {
325 : // Check if return num is ok: it was requested and it does match one of the requested return numbers
326 112133 : const long returnnum= static_cast<const unsigned char *>(pointdata.getChannel(CHANNEL_NAME_ReturnNum)->getData())[i];
327 112133 : sprintf(bufCode, "%d", (int)returnnum);
328 : bReturnNumOK = (papszFilterReturnNums == NULL ? true :
329 112133 : (CSLFindString(papszFilterReturnNums, bufCode)!=-1));
330 112133 : if (!bReturnNumOK && CSLFindString(papszFilterReturnNums, "Last")!=-1)
331 : { // Didn't find an explicit match (e.g. return number "1") so we handle a request for "Last" returns
332 0 : const long numreturns= GetChannelElement<long>(*pointdata.getChannel(CHANNEL_NAME_NumReturns), i);
333 0 : bReturnNumOK = (returnnum == numreturns);
334 : }
335 : }
336 :
337 295490 : return bReturnNumOK && bClassificationOK;
338 :
339 : }
340 :
341 :
342 : template<typename DTYPE>
343 2 : CPLErr MG4LidarRasterBand::doReadBlock(int nBlockXOff, int nBlockYOff, void * pImage)
344 : {
345 2 : MG4LidarDataset * poGDS = (MG4LidarDataset *)poDS;
346 2 : MG4PointReader *reader =poGDS->reader;
347 :
348 : struct Accumulator_t
349 : {
350 : DTYPE value;
351 : int count;
352 : } ;
353 2 : Accumulator_t * Accumulator = NULL;
354 2 : if (EQUAL(Aggregation, "Mean"))
355 : {
356 2 : Accumulator = new Accumulator_t[nBlockXSize*nBlockYSize];
357 2 : memset (Accumulator, 0, sizeof(Accumulator_t)*nBlockXSize*nBlockYSize);
358 : }
359 688 : for (int i = 0; i < nBlockXSize; i++)
360 : {
361 247117 : for (int j = 0; j < nBlockYSize; j++)
362 : {
363 246431 : static_cast<DTYPE*>(pImage)[i*nBlockYSize+j] = static_cast<DTYPE>(nodatavalue);
364 : }
365 : }
366 :
367 : double geoTrans[6];
368 2 : poGDS->GetGeoTransform(geoTrans);
369 2 : double xres = geoTrans[1];
370 2 : double yres = geoTrans[5];
371 :
372 : // Get the extent of the requested block.
373 2 : double xmin = geoTrans[0] + (nBlockXOff *nBlockXSize* xres);
374 2 : double xmax = xmin + nBlockXSize* xres;
375 2 : double ymax = reader->getBounds().y.max - (nBlockYOff * nBlockYSize* -yres);
376 2 : double ymin = ymax - nBlockYSize* -yres;
377 2 : Bounds bounds(xmin, xmax, ymin, ymax, -HUGE_VAL, +HUGE_VAL);
378 2 : PointData pointdata;
379 2 : pointdata.init(reader->getPointInfo(), 4096);
380 2 : double fraction = 1.0/pow(RESOLUTION_RATIO, poGDS->iLevel);
381 2 : CPLDebug( "MG4Lidar", "IReadBlock(x=%d y=%d, level=%d, fraction=%f)", nBlockXOff, nBlockYOff, poGDS->iLevel, fraction);
382 2 : Scoped<PointIterator> iter(reader->createIterator(bounds, fraction, reader->getPointInfo(), NULL));
383 :
384 2 : const double * x = pointdata.getX();
385 2 : const double * y = pointdata.getY();
386 : size_t nPoints;
387 78 : while ( (nPoints = iter->getNextPoints(pointdata)) != 0)
388 : {
389 295564 : for( size_t i = 0; i < nPoints; i++ )
390 : {
391 295490 : const ChannelData * channel = pointdata.getChannel(ChannelName);
392 295490 : if (papszFilterClassCodes || papszFilterReturnNums)
393 : {
394 295490 : if (!ElementPassesFilter(pointdata, i))
395 183357 : continue;
396 : }
397 112133 : double col = (x[i] - xmin) / xres;
398 112133 : double row = (ymax - y[i]) / -yres;
399 112133 : col = floor (col);
400 112133 : row = floor (row);
401 :
402 112133 : if (row < 0)
403 0 : row = 0;
404 112133 : else if (row >= nBlockYSize)
405 0 : row = nBlockYSize - 1;
406 112133 : if (col < 0)
407 0 : col = 0;
408 112133 : else if (col >= nBlockXSize )
409 1 : col = nBlockXSize - 1;
410 :
411 112133 : int iCol = (int) (col);
412 112133 : int iRow = (int) (row);
413 112133 : const int offset =iRow* nBlockXSize + iCol;
414 112133 : if (EQUAL(Aggregation, "Max"))
415 : {
416 0 : DTYPE value = GetChannelElement<DTYPE>(*channel, i);
417 0 : if (static_cast<DTYPE *>(pImage)[offset] == static_cast<DTYPE>(nodatavalue) ||
418 : static_cast<DTYPE *>(pImage)[offset] < value)
419 0 : static_cast<DTYPE *>(pImage)[offset] = value;
420 : }
421 112133 : else if (EQUAL(Aggregation, "Min"))
422 : {
423 0 : DTYPE value = GetChannelElement<DTYPE>(*channel, i);
424 0 : if (static_cast<DTYPE *>(pImage)[offset] == static_cast<DTYPE>(nodatavalue) ||
425 : static_cast<DTYPE *>(pImage)[offset] > value)
426 0 : static_cast<DTYPE *>(pImage)[offset] = value;
427 : }
428 112133 : else if (EQUAL(Aggregation, "Mean"))
429 : {
430 112133 : DTYPE value = GetChannelElement<DTYPE>(*channel, i);
431 112133 : Accumulator[offset].count++;
432 112133 : Accumulator[offset].value += value;
433 112133 : static_cast<DTYPE*>(pImage)[offset] = static_cast<DTYPE>(Accumulator[offset].value / Accumulator[offset].count);
434 : }
435 : }
436 : }
437 :
438 2 : delete[] Accumulator;
439 2 : return CE_None;
440 : }
441 :
442 2 : double MG4LidarRasterBand::getMaxValue()
443 : {
444 : double retval;
445 2 : switch(eDataType)
446 : {
447 : #define DO_CASE(gdt, largestvalue) case (gdt):\
448 : retval = static_cast<double>(largestvalue); \
449 : break;
450 :
451 0 : DO_CASE (GDT_Float64, DBL_MAX);
452 2 : DO_CASE (GDT_Float32, FLT_MAX);
453 0 : DO_CASE (GDT_Int32, INT_MAX);
454 0 : DO_CASE (GDT_UInt32, UINT_MAX);
455 0 : DO_CASE (GDT_Int16, SHRT_MAX);
456 0 : DO_CASE (GDT_UInt16, USHRT_MAX);
457 0 : DO_CASE (GDT_Byte, UCHAR_MAX);
458 : #undef DO_CASE
459 : default:
460 0 : retval = 0;
461 : break;
462 : }
463 2 : return retval;
464 : }
465 : /************************************************************************/
466 : /* IReadBlock() */
467 : /************************************************************************/
468 :
469 2 : CPLErr MG4LidarRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
470 : void * pImage )
471 : {
472 : CPLErr eErr;
473 2 : switch(eDataType)
474 : {
475 : #define DO_CASE(gdt, nativedt) case (gdt):\
476 : eErr = doReadBlock<nativedt>(nBlockXOff, nBlockYOff, pImage); \
477 : break;
478 :
479 0 : DO_CASE (GDT_Float64, double);
480 2 : DO_CASE (GDT_Float32, float);
481 0 : DO_CASE (GDT_Int32, long);
482 0 : DO_CASE (GDT_UInt32, unsigned long);
483 0 : DO_CASE (GDT_Int16, short);
484 0 : DO_CASE (GDT_UInt16, unsigned short);
485 0 : DO_CASE (GDT_Byte, char);
486 : #undef DO_CASE
487 : default:
488 0 : return CE_Failure;
489 : break;
490 :
491 : }
492 2 : return CE_None;
493 : }
494 :
495 : /************************************************************************/
496 : /* GetStatistics() */
497 : /************************************************************************/
498 :
499 0 : CPLErr MG4LidarRasterBand::GetStatistics( int bApproxOK, int bForce,
500 : double *pdfMin, double *pdfMax,
501 : double *pdfMean, double *pdfStdDev )
502 :
503 : {
504 0 : bApproxOK = TRUE;
505 0 : bForce = TRUE;
506 :
507 : return GDALPamRasterBand::GetStatistics( bApproxOK, bForce,
508 : pdfMin, pdfMax,
509 0 : pdfMean, pdfStdDev );
510 :
511 : }
512 : /************************************************************************/
513 : /* GetNoDataValue() */
514 : /************************************************************************/
515 :
516 0 : double MG4LidarRasterBand::GetNoDataValue( int *pbSuccess )
517 : {
518 0 : if (pbSuccess)
519 0 : *pbSuccess = TRUE;
520 0 : return nodatavalue;
521 : }
522 :
523 :
524 : /************************************************************************/
525 : /* MG4LidarDataset() */
526 : /************************************************************************/
527 2 : MG4LidarDataset::MG4LidarDataset()
528 : {
529 2 : reader = NULL;
530 2 : fileIO = NULL;
531 :
532 2 : poXMLPCView = NULL;
533 2 : ownsXML = false;
534 2 : nOverviewCount = 0;
535 2 : papoOverviewDS = NULL;
536 :
537 2 : }
538 : /************************************************************************/
539 : /* ~MG4LidarDataset() */
540 : /************************************************************************/
541 :
542 2 : MG4LidarDataset::~MG4LidarDataset()
543 :
544 : {
545 2 : FlushCache();
546 2 : if ( papoOverviewDS )
547 : {
548 2 : for( int i = 0; i < nOverviewCount; i++ )
549 1 : delete papoOverviewDS[i];
550 1 : CPLFree( papoOverviewDS );
551 : }
552 2 : if (ownsXML)
553 1 : CPLDestroyXMLNode(poXMLPCView);
554 :
555 2 : RELEASE(reader);
556 2 : RELEASE(fileIO);
557 2 : }
558 :
559 : /************************************************************************/
560 : /* GetGeoTransform() */
561 : /************************************************************************/
562 :
563 3 : CPLErr MG4LidarDataset::GetGeoTransform( double * padfTransform )
564 :
565 : {
566 3 : padfTransform[0] = reader->getBounds().x.min ;// Upper left X, Y
567 3 : padfTransform[3] = reader->getBounds().y.max; //
568 3 : padfTransform[1] = reader->getBounds().x.length()/GetRasterXSize(); //xRes
569 3 : padfTransform[2] = 0.0;
570 :
571 3 : padfTransform[4] = 0.0;
572 3 : padfTransform[5] = -1 * reader->getBounds().y.length()/GetRasterYSize(); //yRes
573 :
574 3 : return CE_None;
575 : }
576 :
577 : /************************************************************************/
578 : /* GetProjectionRef() */
579 : /************************************************************************/
580 :
581 1 : const char *MG4LidarDataset::GetProjectionRef()
582 :
583 : {
584 1 : const char * wkt = CPLGetXMLValue(poXMLPCView, "GeoReference", NULL);
585 1 : if (wkt == NULL)
586 1 : wkt = reader->getWKT();
587 1 : return(wkt);
588 : }
589 :
590 : /************************************************************************/
591 : /* OpenZoomLevel() */
592 : /************************************************************************/
593 :
594 2 : CPLErr MG4LidarDataset::OpenZoomLevel( int iZoom )
595 : {
596 : /* -------------------------------------------------------------------- */
597 : /* Get image geometry. */
598 : /* -------------------------------------------------------------------- */
599 2 : iLevel = iZoom;
600 :
601 : // geo dimensions
602 2 : const double gWidth = reader->getBounds().x.length() ;
603 2 : const double gHeight = reader->getBounds().y.length() ;
604 :
605 : // geo res
606 2 : double xRes = pow(RESOLUTION_RATIO, iZoom) * gWidth / MaxRasterSize ;
607 2 : double yRes = pow(RESOLUTION_RATIO, iZoom) * gHeight / MaxRasterSize ;
608 2 : xRes = yRes = MAX(xRes, yRes);
609 :
610 : // pixel dimensions
611 2 : nRasterXSize = static_cast<int>(gWidth / xRes + 0.5);
612 2 : nRasterYSize = static_cast<int>(gHeight / yRes + 0.5);
613 :
614 2 : nBlockXSize = static_cast<int>(MIN(MaxBlockSideSize , GetRasterXSize()));
615 2 : nBlockYSize = static_cast<int>(MIN(MaxBlockSideSize , GetRasterYSize()));
616 :
617 : CPLDebug( "MG4Lidar", "Opened zoom level %d with size %dx%d.\n",
618 2 : iZoom, nRasterXSize, nRasterYSize );
619 :
620 :
621 :
622 : /* -------------------------------------------------------------------- */
623 : /* Handle sample type and color space. */
624 : /* -------------------------------------------------------------------- */
625 : //eColorSpace = poImageReader->getColorSpace();
626 : /* -------------------------------------------------------------------- */
627 : /* Create band information objects. */
628 : /* -------------------------------------------------------------------- */
629 2 : size_t BandCount = 0;
630 2 : CPLXMLNode* xmlBand = poXMLPCView;
631 2 : bool bClass = false;
632 2 : bool bNumRets = false;
633 2 : bool bRetNum = false;
634 6 : while ((xmlBand = CPLSearchXMLNode(xmlBand, "Band")) != NULL)
635 : {
636 2 : CPLXMLNode * xmlChannel = CPLSearchXMLNode(xmlBand, "Channel");
637 2 : const char * name = "Z";
638 2 : if (xmlChannel && xmlChannel->psChild && xmlChannel->psChild->pszValue)
639 2 : name = xmlChannel->psChild->pszValue;
640 :
641 2 : BandCount++;
642 2 : MG4LidarRasterBand *band = new MG4LidarRasterBand(this, BandCount, xmlBand, name);
643 2 : SetBand(BandCount, band);
644 4 : if (band->papszFilterClassCodes) bClass = true;
645 2 : if (band->papszFilterReturnNums) bNumRets = true;
646 2 : if (bNumRets && CSLFindString(band->papszFilterReturnNums, "Last")) bRetNum = true;
647 2 : xmlBand = xmlBand->psNext;
648 : }
649 2 : nBands = BandCount;
650 2 : int nSDKChannels = BandCount + (bClass ? 1 : 0) + (bNumRets ? 1 : 0) + (bRetNum ? 1 : 0);
651 2 : if (BandCount == 0) // default if no bands specified.
652 : {
653 0 : MG4LidarRasterBand *band = new MG4LidarRasterBand(this, 1, NULL, CHANNEL_NAME_Z);
654 0 : SetBand(1, band);
655 0 : nBands = 1;
656 0 : nSDKChannels = 1;
657 : }
658 2 : requiredChannels.init(nSDKChannels);
659 2 : const ChannelInfo *ci = NULL;
660 4 : for (int i=0; i<nBands; i++)
661 : {
662 2 : ci = reader->getChannel(dynamic_cast<MG4LidarRasterBand*>(papoBands[i])->ChannelName);
663 2 : requiredChannels.getChannel(i).init(*ci);
664 : }
665 2 : int iSDKChannels = nBands;
666 2 : if (bClass)
667 : {
668 2 : ci = reader->getChannel(CHANNEL_NAME_ClassId);
669 2 : requiredChannels.getChannel(iSDKChannels++).init(*ci);
670 : }
671 2 : if (bRetNum)
672 : {
673 0 : ci = reader->getChannel(CHANNEL_NAME_ReturnNum);
674 0 : requiredChannels.getChannel(iSDKChannels++).init(*ci);
675 : }
676 2 : if (bNumRets)
677 : {
678 0 : ci = reader->getChannel(CHANNEL_NAME_NumReturns);
679 0 : requiredChannels.getChannel(iSDKChannels++).init(*ci);
680 : }
681 :
682 2 : return CE_None;
683 : }
684 :
685 : /************************************************************************/
686 : /* Open() */
687 : /************************************************************************/
688 :
689 12617 : GDALDataset *MG4LidarDataset::Open( GDALOpenInfo * poOpenInfo )
690 :
691 : {
692 : #ifdef notdef
693 : CPLPushErrorHandler( CPLLoggingErrorHandler );
694 : CPLSetConfigOption( "CPL_DEBUG", "ON" );
695 : CPLSetConfigOption( "CPL_LOG", "C:\\ArcGIS_GDAL\\jdem\\cpl.log" );
696 : #endif
697 :
698 12617 : if( poOpenInfo->fp == NULL || poOpenInfo->nHeaderBytes < 32 )
699 11682 : return NULL;
700 :
701 : CPLXMLNode *pxmlPCView;
702 :
703 : // do something sensible for .sid files without a .view
704 935 : if( EQUALN((const char *) poOpenInfo->pabyHeader, "msid", 4) )
705 : {
706 : int gen;
707 : bool raster;
708 0 : if( !Version::getMrSIDFileVersion(poOpenInfo->pabyHeader, gen, raster)
709 : || raster )
710 0 : return NULL;
711 :
712 0 : CPLString xmltmp( "<PointCloudView><InputFile>" );
713 0 : xmltmp.append( poOpenInfo->pszFilename );
714 0 : xmltmp.append( "</InputFile></PointCloudView>" );
715 0 : pxmlPCView = CPLParseXMLString( xmltmp );
716 0 : if (pxmlPCView == NULL)
717 0 : return NULL;
718 : }
719 : else
720 : {
721 : // support .view xml
722 935 : if( !EQUALN((const char *) poOpenInfo->pabyHeader, "<PointCloudView", 15 ) )
723 934 : return NULL;
724 :
725 1 : pxmlPCView = CPLParseXMLFile( poOpenInfo->pszFilename );
726 1 : if (pxmlPCView == NULL)
727 0 : return NULL;
728 : }
729 :
730 1 : CPLXMLNode *psInputFile = CPLGetXMLNode( pxmlPCView, "InputFile" );
731 1 : if( psInputFile == NULL )
732 : {
733 : CPLError( CE_Failure, CPLE_OpenFailed,
734 0 : "Failed to find <InputFile> in document." );
735 0 : CPLDestroyXMLNode(pxmlPCView);
736 0 : return NULL;
737 : }
738 1 : CPLString sidInputName(psInputFile->psChild->pszValue);
739 1 : if (CPLIsFilenameRelative(sidInputName))
740 : {
741 1 : CPLString dirname(CPLGetDirname(poOpenInfo->pszFilename));
742 1 : sidInputName = CPLString(CPLFormFilename(dirname, sidInputName, NULL));
743 : }
744 1 : GDALOpenInfo openinfo(sidInputName, GA_ReadOnly);
745 :
746 : /* -------------------------------------------------------------------- */
747 : /* Check that particular fields in the header are valid looking */
748 : /* dates. */
749 : /* -------------------------------------------------------------------- */
750 1 : if( openinfo.fp == NULL || openinfo.nHeaderBytes < 50 )
751 : {
752 0 : CPLDestroyXMLNode(pxmlPCView);
753 0 : return NULL;
754 : }
755 :
756 : /* check magic */
757 : // to do: SDK should provide an API for this.
758 1 : if( !EQUALN((const char *) openinfo.pabyHeader, "msid", 4)
759 : || (*(openinfo.pabyHeader+4) != 0x4 )) // Generation 4. ... is there more we can check?
760 : {
761 0 : CPLDestroyXMLNode(pxmlPCView);
762 0 : return NULL;
763 : }
764 :
765 : /* -------------------------------------------------------------------- */
766 : /* Create a corresponding GDALDataset. */
767 : /* -------------------------------------------------------------------- */
768 : MG4LidarDataset *poDS;
769 :
770 1 : poDS = new MG4LidarDataset();
771 1 : poDS->poXMLPCView = pxmlPCView;
772 1 : poDS->ownsXML = true;
773 2 : poDS->reader = CropableMG4PointReader::create();
774 1 : poDS->fileIO = FileIO::create();
775 :
776 1 : const char * pszClipExtent = CPLGetXMLValue(pxmlPCView, "ClipBox", NULL);
777 1 : MG4PointReader *r = MG4PointReader::create();
778 1 : FileIO* io = FileIO::create();
779 :
780 : #if (defined(WIN32) && _MSC_VER >= 1310) || __MSVCRT_VERSION__ >= 0x0601
781 : bool bIsUTF8 = CSLTestBoolean( CPLGetConfigOption( "GDAL_FILENAME_IS_UTF8", "YES" ) );
782 : wchar_t *pwszFilename = NULL;
783 : if (bIsUTF8)
784 : {
785 : pwszFilename = CPLRecodeToWChar(openinfo.pszFilename, CPL_ENC_UTF8, CPL_ENC_UCS2);
786 : if (!pwszFilename)
787 : {
788 : RELEASE(r);
789 : RELEASE(io);
790 : return NULL;
791 : }
792 : io->init(pwszFilename, "r");
793 : }
794 : else
795 : io->init(openinfo.pszFilename, "r");
796 : #else
797 1 : io->init(openinfo.pszFilename, "r");
798 : #endif
799 1 : r->init(io);
800 1 : Bounds bounds = r->getBounds();
801 1 : if (pszClipExtent)
802 : {
803 0 : char ** papszClipExtent = CSLTokenizeString(pszClipExtent);
804 0 : int cslcount = CSLCount(papszClipExtent);
805 0 : if (cslcount != 4 && cslcount != 6)
806 : {
807 : CPLError( CE_Failure, CPLE_OpenFailed,
808 0 : "Invalid ClipBox. Must contain 4 or 6 floats." );
809 0 : CSLDestroy(papszClipExtent);
810 0 : delete poDS;
811 0 : RELEASE(r);
812 0 : RELEASE(io);
813 : #if (defined(WIN32) && _MSC_VER >= 1310) || __MSVCRT_VERSION__ >= 0x0601
814 : if ( pwszFilename )
815 : CPLFree( pwszFilename );
816 : #endif
817 0 : return NULL;
818 : }
819 0 : if (!EQUAL(papszClipExtent[0], "NOFILTER"))
820 0 : bounds.x.min = atof(papszClipExtent[0]);
821 0 : if (!EQUAL(papszClipExtent[1], "NOFILTER"))
822 0 : bounds.x.max = atof(papszClipExtent[1]);
823 0 : if (!EQUAL(papszClipExtent[2], "NOFILTER"))
824 0 : bounds.y.min = atof(papszClipExtent[2]);
825 0 : if (!EQUAL(papszClipExtent[3], "NOFILTER"))
826 0 : bounds.y.max = atof(papszClipExtent[3]);
827 0 : if (cslcount == 6)
828 : {
829 0 : if (!EQUAL(papszClipExtent[4], "NOFILTER"))
830 0 : bounds.z.min = atof(papszClipExtent[4]);
831 0 : if (!EQUAL(papszClipExtent[5], "NOFILTER"))
832 0 : bounds.z.max = atof(papszClipExtent[5]);
833 : }
834 0 : CSLDestroy(papszClipExtent);
835 : }
836 :
837 : #if (defined(WIN32) && _MSC_VER >= 1310) || __MSVCRT_VERSION__ >= 0x0601
838 : if (bIsUTF8)
839 : {
840 : poDS->fileIO->init(pwszFilename, "r");
841 : CPLFree(pwszFilename);
842 : }
843 : else
844 : poDS->fileIO->init( openinfo.pszFilename, "r" );
845 : #else
846 1 : poDS->fileIO->init( openinfo.pszFilename, "r" );
847 : #endif
848 1 : dynamic_cast<CropableMG4PointReader *>(poDS->reader)->init(poDS->fileIO, &bounds);
849 1 : poDS->SetDescription(poOpenInfo->pszFilename);
850 1 : poDS->TryLoadXML();
851 :
852 1 : double pts_per_area = ((double)r->getNumPoints())/(r->getBounds().x.length()*r->getBounds().y.length());
853 1 : double average_pt_spacing = sqrt(1.0 / pts_per_area) ;
854 1 : double cell_side = average_pt_spacing;
855 1 : const char * pszCellSize = CPLGetXMLValue(pxmlPCView, "CellSize", NULL);
856 1 : if (pszCellSize)
857 0 : cell_side = atof(pszCellSize);
858 1 : MaxRasterSize = MAX(poDS->reader->getBounds().x.length()/cell_side, poDS->reader->getBounds().y.length()/cell_side);
859 :
860 1 : RELEASE(r);
861 1 : RELEASE(io);
862 :
863 : // Calculate the number of levels to expose. The highest level correpsonds to a
864 : // raster size of 256 on the longest side.
865 1 : double blocksizefactor = MaxRasterSize/256.0;
866 1 : poDS->nOverviewCount = MAX(0, (int)(log(blocksizefactor)/log(RESOLUTION_RATIO) + 0.5));
867 1 : if ( poDS->nOverviewCount > 0 )
868 : {
869 : int i;
870 :
871 : poDS->papoOverviewDS = (MG4LidarDataset **)
872 1 : CPLMalloc( poDS->nOverviewCount * (sizeof(void*)) );
873 :
874 2 : for ( i = 0; i < poDS->nOverviewCount; i++ )
875 : {
876 1 : poDS->papoOverviewDS[i] = new MG4LidarDataset ();
877 2 : poDS->papoOverviewDS[i]->reader = RETAIN(poDS->reader);
878 1 : poDS->papoOverviewDS[i]->SetMetadata(poDS->GetMetadata("MG4Lidar"), "MG4Lidar");
879 1 : poDS->papoOverviewDS[i]->poXMLPCView = pxmlPCView;
880 1 : poDS->papoOverviewDS[i]->OpenZoomLevel( i+1 );
881 : }
882 : }
883 :
884 : /* -------------------------------------------------------------------- */
885 : /* Create object for the whole image. */
886 : /* -------------------------------------------------------------------- */
887 1 : poDS->OpenZoomLevel( 0 );
888 :
889 : CPLDebug( "MG4Lidar",
890 : "Opened image: width %d, height %d, bands %d",
891 1 : poDS->nRasterXSize, poDS->nRasterYSize, poDS->nBands );
892 :
893 1 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
894 : {
895 0 : delete poDS;
896 0 : return NULL;
897 : }
898 :
899 1 : if (! ((poDS->nBands == 1) || (poDS->nBands == 3)))
900 : {
901 : CPLDebug( "MG4Lidar",
902 0 : "Inappropriate number of bands (%d)", poDS->nBands );
903 0 : delete poDS;
904 0 : return(NULL);
905 : }
906 :
907 1 : return( poDS );
908 : }
909 :
910 : /************************************************************************/
911 : /* GDALRegister_MG4Lidar() */
912 : /************************************************************************/
913 :
914 582 : void GDALRegister_MG4Lidar()
915 :
916 : {
917 : GDALDriver *poDriver;
918 :
919 582 : if (! GDAL_CHECK_VERSION("MG4Lidar driver"))
920 0 : return;
921 :
922 582 : if( GDALGetDriverByName( "MG4Lidar" ) == NULL )
923 : {
924 561 : poDriver = new GDALDriver();
925 :
926 561 : poDriver->SetDescription( "MG4Lidar" );
927 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
928 561 : "MrSID Generation 4 / Lidar (.sid)" );
929 : // To do: update this help file in gdal.org
930 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
931 561 : "frmt_mrsid_lidar.html" );
932 :
933 561 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "view" );
934 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
935 561 : "Float64" );
936 :
937 561 : poDriver->pfnOpen = MG4LidarDataset::Open;
938 :
939 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
940 : }
941 : }
|
