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