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