Context Navigation

DataSet.cpp @ 3309

Revision 3309, 135.2 kB (checked in by robert, 9 years ago)
Added protection to prevent crashes on calls when no data is set up
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`

Line
1	/ --c++-- OpenSceneGraph - Copyright (C) 1998-2003 Robert Osfield*
2	*
3	* This library is open source and may be redistributed and/or modified under
4	* the terms of the OpenSceneGraph Public License (OSGPL) version 0.0 or
5	* (at your option) any later version. The full license is in LICENSE file
6	* included with this distribution, and on the openscenegraph.org website.
7	*
8	* This library is distributed in the hope that it will be useful,
9	* but WITHOUT ANY WARRANTY; without even the implied warranty of
10	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11	* OpenSceneGraph Public License for more details.
12	*/
13
14
15	#include <osg/Geode>
16	#include <osg/ShapeDrawable>
17	#include <osg/Texture2D>
18	#include <osg/Group>
19	#include <osg/Geometry>
20	#include <osg/MatrixTransform>
21	#include <osg/ClusterCullingCallback>
22	#include <osg/Notify>
23
24	#include <osgUtil/SmoothingVisitor>
25	#include <osgUtil/TriStripVisitor>
26	#include <osgUtil/Simplifier>
27
28	#include <osgDB/ReadFile>
29	#include <osgDB/WriteFile>
30	#include <osgDB/FileNameUtils>
31
32	#include <osgTerrain/DataSet>
33
34	// GDAL includes
35	#include <gdal_priv.h>
36	#include <cpl_string.h>
37	#include <gdalwarper.h>
38	#include <ogr_spatialref.h>
39
40	// standard library includes
41	#include <sstream>
42	#include <iostream>
43
44
45	using namespace osgTerrain;
46
47
48	enum CoordinateSystemType
49	{
50	GEOCENTRIC,
51	GEOGRAPHIC,
52	PROJECTED,
53	LOCAL
54	};
55
56	CoordinateSystemType getCoordinateSystemType(const osg::CoordinateSystemNode* lhs)
57	{
58	if (!lhs) return PROJECTED;
59
60	// set up LHS SpatialReference
61	char* projection_string = strdup(lhs->getCoordinateSystem().c_str());
62	char* importString = projection_string;
63
64	OGRSpatialReference lhsSR;
65	lhsSR.importFromWkt(&importString);
66
67
68
69	osg::notify(osg::INFO)<<"getCoordinateSystemType("<<projection_string<<")"<<std::endl;
70	osg::notify(osg::INFO)<<" lhsSR.IsGeographic()="<<lhsSR.IsGeographic()<<std::endl;
71	osg::notify(osg::INFO)<<" lhsSR.IsProjected()="<<lhsSR.IsProjected()<<std::endl;
72	osg::notify(osg::INFO)<<" lhsSR.IsLocal()="<<lhsSR.IsLocal()<<std::endl;
73
74	free(projection_string);
75
76	if (strcmp(lhsSR.GetRoot()->GetValue(),"GEOCCS")==0) osg::notify(osg::INFO)<<" lhsSR. is GEOCENTRIC "<<std::endl;
77
78
79	if (strcmp(lhsSR.GetRoot()->GetValue(),"GEOCCS")==0) return GEOCENTRIC;
80	if (lhsSR.IsGeographic()) return GEOGRAPHIC;
81	if (lhsSR.IsProjected()) return PROJECTED;
82	if (lhsSR.IsLocal()) return LOCAL;
83	return PROJECTED;
84	}
85
86	double getAngularUnits(const osg::CoordinateSystemNode* lhs)
87	{
88	// set up LHS SpatialReference
89	char* projection_string = strdup(lhs->getCoordinateSystem().c_str());
90	char* importString = projection_string;
91
92	OGRSpatialReference lhsSR;
93	lhsSR.importFromWkt(&importString);
94
95	free(projection_string);
96
97	char* str;
98	double result = lhsSR.GetAngularUnits(&str);
99	osg::notify(osg::INFO)<<"lhsSR.GetAngularUnits("<<str<<") "<<result<<std::endl;
100
101	return result;
102	}
103
104	double getLinearUnits(const osg::CoordinateSystemNode* lhs)
105	{
106	// set up LHS SpatialReference
107	char* projection_string = strdup(lhs->getCoordinateSystem().c_str());
108	char* importString = projection_string;
109
110	OGRSpatialReference lhsSR;
111	lhsSR.importFromWkt(&importString);
112
113	free(projection_string);
114
115	char* str;
116	double result = lhsSR.GetLinearUnits(&str);
117	osg::notify(osg::INFO)<<"lhsSR.GetLinearUnits("<<str<<") "<<result<<std::endl;
118
119	osg::notify(osg::INFO)<<"lhsSR.IsGeographic() "<<lhsSR.IsGeographic()<<std::endl;
120	osg::notify(osg::INFO)<<"lhsSR.IsProjected() "<<lhsSR.IsProjected()<<std::endl;
121	osg::notify(osg::INFO)<<"lhsSR.IsLocal() "<<lhsSR.IsLocal()<<std::endl;
122
123	return result;
124	}
125
126	bool areCoordinateSystemEquivalent(const osg::CoordinateSystemNode* lhs,const osg::CoordinateSystemNode* rhs)
127	{
128	// if ptr's equal the return true
129	if (lhs == rhs) return true;
130
131	// if one CS is NULL then true false
132	if (!lhs \|\| !rhs)
133	{
134	osg::notify(osg::INFO)<<"areCoordinateSystemEquivalent lhs="<<lhs<<" rhs="<<rhs<<" return true"<<std::endl;
135	return false;
136	}
137
138	osg::notify(osg::INFO)<<"areCoordinateSystemEquivalent lhs="<<lhs->getCoordinateSystem()<<" rhs="<<rhs->getCoordinateSystem()<<std::endl;
139
140	// use compare on ProjectionRef strings.
141	if (lhs->getCoordinateSystem() == rhs->getCoordinateSystem()) return true;
142
143	// set up LHS SpatialReference
144	char* projection_string = strdup(lhs->getCoordinateSystem().c_str());
145	char* importString = projection_string;
146
147	OGRSpatialReference lhsSR;
148	lhsSR.importFromWkt(&importString);
149
150	free(projection_string);
151
152	// set up RHS SpatialReference
153	projection_string = strdup(rhs->getCoordinateSystem().c_str());
154	importString = projection_string;
155
156	OGRSpatialReference rhsSR;
157	rhsSR.importFromWkt(&importString);
158
159	free(projection_string);
160
161	int result = lhsSR.IsSame(&rhsSR);
162
163	#if 0
164	int result2 = lhsSR.IsSameGeogCS(&rhsSR);
165
166	osg::notify(osg::INFO)<<"areCoordinateSystemEquivalent "<<std::endl
167	<<"LHS = "<<lhs->getCoordinateSystem()<<std::endl
168	<<"RHS = "<<rhs->getCoordinateSystem()<<std::endl
169	<<"result = "<<result<<" result2 = "<<result2<<std::endl;
170	#endif
171	return result ? true : false;
172	}
173
174	DataSet::SourceData* DataSet::SourceData::readData(Source* source)
175	{
176	if (!source) return 0;
177
178
179	switch(source->getType())
180	{
181	case(Source::IMAGE):
182	case(Source::HEIGHT_FIELD):
183	{
184	GDALDataset* gdalDataSet = (GDALDataset*)GDALOpen(source->getFileName().c_str(),GA_ReadOnly);
185	if (gdalDataSet)
186	{
187	SourceData* data = new SourceData(source);
188	data->_gdalDataSet = gdalDataSet;
189
190	data->_numValuesX = gdalDataSet->GetRasterXSize();
191	data->_numValuesY = gdalDataSet->GetRasterYSize();
192	data->_numValuesZ = gdalDataSet->GetRasterCount();
193	data->_hasGCPs = gdalDataSet->GetGCPCount()!=0;
194
195	const char* pszSourceSRS = gdalDataSet->GetProjectionRef();
196	if (!pszSourceSRS \|\| strlen(pszSourceSRS)==0) pszSourceSRS = gdalDataSet->GetGCPProjection();
197
198	data->_cs = new osg::CoordinateSystemNode("WKT",pszSourceSRS);
199
200	double geoTransform[6];
201	if (gdalDataSet->GetGeoTransform(geoTransform)==CE_None)
202	{
203	data->_geoTransform.set( geoTransform[1], geoTransform[4], 0.0, 0.0,
204	geoTransform[2], geoTransform[5], 0.0, 0.0,
205	0.0, 0.0, 1.0, 0.0,
206	geoTransform[0], geoTransform[3], 0.0, 1.0);
207
208	data->computeExtents();
209
210	}
211	else if (gdalDataSet->GetGCPCount()>0 && gdalDataSet->GetGCPProjection())
212	{
213	osg::notify(osg::INFO) << " Using GCP's"<< std::endl;
214
215
216	/ -------------------------------------------------------------------- /
217	/ Create a transformation object from the source to /
218	/ destination coordinate system. /
219	/ -------------------------------------------------------------------- /
220	void *hTransformArg =
221	GDALCreateGenImgProjTransformer( gdalDataSet, pszSourceSRS,
222	NULL, pszSourceSRS,
223	TRUE, 0.0, 1 );
224
225	if ( hTransformArg == NULL )
226	{
227	osg::notify(osg::INFO)<<" failed to create transformer"<<std::endl;
228	return NULL;
229	}
230
231	/ -------------------------------------------------------------------- /
232	/ Get approximate output definition. /
233	/ -------------------------------------------------------------------- /
234	double adfDstGeoTransform[6];
235	int nPixels=0, nLines=0;
236	if( GDALSuggestedWarpOutput( gdalDataSet,
237	GDALGenImgProjTransform, hTransformArg,
238	adfDstGeoTransform, &nPixels, &nLines )
239	!= CE_None )
240	{
241	osg::notify(osg::INFO)<<" failed to create warp"<<std::endl;
242	return NULL;
243	}
244
245	GDALDestroyGenImgProjTransformer( hTransformArg );
246
247
248	data->_geoTransform.set( adfDstGeoTransform[1], adfDstGeoTransform[4], 0.0, 0.0,
249	adfDstGeoTransform[2], adfDstGeoTransform[5], 0.0, 0.0,
250	0.0, 0.0, 1.0, 0.0,
251	adfDstGeoTransform[0], adfDstGeoTransform[3], 0.0, 1.0);
252
253	data->computeExtents();
254
255	}
256	else
257	{
258	osg::notify(osg::INFO) << " No GeoTransform or GCP's - unable to compute position in space"<< std::endl;
259
260	data->_geoTransform.set( 1.0, 0.0, 0.0, 0.0,
261	0.0, 1.0, 0.0, 0.0,
262	0.0, 0.0, 1.0, 0.0,
263	0.0, 0.0, 0.0, 1.0);
264
265	data->computeExtents();
266
267	}
268	return data;
269	}
270	}
271	case(Source::MODEL):
272	{
273	osg::Node* model = osgDB::readNodeFile(source->getFileName().c_str());
274	if (model)
275	{
276	SourceData* data = new SourceData(source);
277	data->_model = model;
278	data->_extents.expandBy(model->getBound());
279	}
280
281	}
282	break;
283	}
284
285	return 0;
286	}
287
288	osg::BoundingBox DataSet::SourceData::getExtents(const osg::CoordinateSystemNode* cs) const
289	{
290	return computeSpatialProperties(cs)._extents;
291	}
292
293	const DataSet::SpatialProperties& DataSet::SourceData::computeSpatialProperties(const osg::CoordinateSystemNode* cs) const
294	{
295	// check to see it exists in the _spatialPropertiesMap first.
296	SpatialPropertiesMap::const_iterator itr = _spatialPropertiesMap.find(cs);
297	if (itr!=_spatialPropertiesMap.end())
298	{
299	return itr->second;
300	}
301
302	if (areCoordinateSystemEquivalent(_cs.get(),cs))
303	{
304	return *this;
305	}
306
307	if (_cs.valid() && cs)
308	{
309
310	if (_gdalDataSet)
311	{
312
313	//osg::notify(osg::INFO)<<"Projecting bounding volume for "<<_source->getFileName()<<std::endl;
314
315
316	// insert into the _spatialPropertiesMap for future reuse.
317	_spatialPropertiesMap[cs] = *this;
318	DataSet::SpatialProperties& sp = _spatialPropertiesMap[cs];
319
320	/ -------------------------------------------------------------------- /
321	/ Create a transformation object from the source to /
322	/ destination coordinate system. /
323	/ -------------------------------------------------------------------- /
324	void *hTransformArg =
325	GDALCreateGenImgProjTransformer( _gdalDataSet,_cs->getCoordinateSystem().c_str(),
326	NULL, cs->getCoordinateSystem().c_str(),
327	TRUE, 0.0, 1 );
328
329	if (!hTransformArg)
330	{
331	osg::notify(osg::INFO)<<" failed to create transformer"<<std::endl;
332	return sp;
333	}
334
335	double adfDstGeoTransform[6];
336	int nPixels=0, nLines=0;
337	if( GDALSuggestedWarpOutput( _gdalDataSet,
338	GDALGenImgProjTransform, hTransformArg,
339	adfDstGeoTransform, &nPixels, &nLines )
340	!= CE_None )
341	{
342	osg::notify(osg::INFO)<<" failed to create warp"<<std::endl;
343	return sp;
344	}
345
346	sp._numValuesX = nPixels;
347	sp._numValuesY = nLines;
348	sp._cs = const_cast<osg::CoordinateSystemNode*>(cs);
349	sp._geoTransform.set( adfDstGeoTransform[1], adfDstGeoTransform[4], 0.0, 0.0,
350	adfDstGeoTransform[2], adfDstGeoTransform[5], 0.0, 0.0,
351	0.0, 0.0, 1.0, 0.0,
352	adfDstGeoTransform[0], adfDstGeoTransform[3], 0.0, 1.0);
353
354	GDALDestroyGenImgProjTransformer( hTransformArg );
355
356	sp.computeExtents();
357
358	return sp;
359	}
360
361	}
362	osg::notify(osg::INFO)<<"DataSet::DataSource::assuming compatible coordinates."<<std::endl;
363	return *this;
364	}
365
366	bool DataSet::SourceData::intersects(const SpatialProperties& sp) const
367	{
368	return sp._extents.intersects(getExtents(sp._cs.get()));
369	}
370
371	void DataSet::SourceData::read(DestinationData& destination)
372	{
373	osg::notify(osg::INFO)<<"A"<<std::endl;
374
375	if (!_source) return;
376
377	osg::notify(osg::INFO)<<"B"<<std::endl;
378
379	switch (_source->getType())
380	{
381	case(Source::IMAGE):
382	osg::notify(osg::INFO)<<"B.1"<<std::endl;
383	readImage(destination);
384	break;
385	case(Source::HEIGHT_FIELD):
386	osg::notify(osg::INFO)<<"B.2"<<std::endl;
387	readHeightField(destination);
388	break;
389	case(Source::MODEL):
390	osg::notify(osg::INFO)<<"B.3"<<std::endl;
391	readModels(destination);
392	break;
393	}
394	osg::notify(osg::INFO)<<"C"<<std::endl;
395	}
396
397	void DataSet::SourceData::readImage(DestinationData& destination)
398	{
399	if (destination._image.valid())
400	{
401	osg::BoundingBox s_bb = getExtents(destination._cs.get());
402
403	osg::BoundingBox d_bb = destination._extents;
404
405	osg::BoundingBox intersect_bb(s_bb.intersect(d_bb));
406	if (!intersect_bb.valid())
407	{
408	osg::notify(osg::INFO)<<"Reading image but it does not intesection destination - ignoring"<<std::endl;
409	return;
410	}
411
412
413	int windowX = osg::maximum((int)floorf((float)_numValuesX*(intersect_bb.xMin()-s_bb.xMin())/(s_bb.xMax()-s_bb.xMin())),0);
414	int windowY = osg::maximum((int)floorf((float)_numValuesY*(intersect_bb.yMin()-s_bb.yMin())/(s_bb.yMax()-s_bb.yMin())),0);
415	int windowWidth = osg::minimum((int)ceilf((float)_numValuesX(intersect_bb.xMax()-s_bb.xMin())/(s_bb.xMax()-s_bb.xMin())),(int*)_numValuesX)-windowX;
416	int windowHeight = osg::minimum((int)ceilf((float)_numValuesY(intersect_bb.yMax()-s_bb.yMin())/(s_bb.yMax()-s_bb.yMin())),(int*)_numValuesY)-windowY;
417
418	int destX = osg::maximum((int)floorf((float)destination._image->s()*(intersect_bb.xMin()-d_bb.xMin())/(d_bb.xMax()-d_bb.xMin())),0);
419	int destY = osg::maximum((int)floorf((float)destination._image->t()*(intersect_bb.yMin()-d_bb.yMin())/(d_bb.yMax()-d_bb.yMin())),0);
420	int destWidth = osg::minimum((int)ceilf((float)destination._image->s()(intersect_bb.xMax()-d_bb.xMin())/(d_bb.xMax()-d_bb.xMin())),(int*)destination._image->s())-destX;
421	int destHeight = osg::minimum((int)ceilf((float)destination._image->t()(intersect_bb.yMax()-d_bb.yMin())/(d_bb.yMax()-d_bb.yMin())),(int*)destination._image->t())-destY;
422
423	osg::notify(osg::INFO)<<" copying from "<<windowX<<"\t"<<windowY<<"\t"<<windowWidth<<"\t"<<windowHeight<<std::endl;
424	osg::notify(osg::INFO)<<" to "<<destX<<"\t"<<destY<<"\t"<<destWidth<<"\t"<<destHeight<<std::endl;
425
426	bool hasRGB = _gdalDataSet->GetRasterCount() >= 3;
427	bool hasColorTable = _gdalDataSet->GetRasterCount() >= 1 && _gdalDataSet->GetRasterBand(1)->GetColorTable();
428	bool hasGreyScale = _gdalDataSet->GetRasterCount() == 1;
429
430	if (hasRGB \|\| hasColorTable \|\| hasGreyScale)
431	{
432	// RGB
433
434	unsigned int numBytesPerPixel = 1;
435	GDALDataType targetGDALType = GDT_Byte;
436
437	int pixelSpace=3*numBytesPerPixel;
438	int lineSpace=-(int)(destination._image->getRowSizeInBytes());
439
440	unsigned char* imageData = destination._image->data(destX,destY+destHeight-1);
441	osg::notify(osg::INFO) << "reading RGB"<<std::endl;
442
443	unsigned char* tempImage = new unsigned char[destWidthdestHeight3];
444
445
446	/ New code courtesy of Frank Warmerdam of the GDAL group /
447
448	// RGB images ... or at least we assume 3+ band images can be treated
449	// as RGB.
450	if( hasRGB )
451	{
452	GDALRasterBand bandRed, bandGreen, *bandBlue;
453
454
455	bandRed = _gdalDataSet->GetRasterBand(1);
456	bandGreen = _gdalDataSet->GetRasterBand(2);
457	bandBlue = _gdalDataSet->GetRasterBand(3);
458
459
460	bandRed->RasterIO(GF_Read, windowX,_numValuesY-(windowY+windowHeight),
461	windowWidth,windowHeight,
462	(void*)(tempImage+0),destWidth,destHeight,
463	targetGDALType,pixelSpace,pixelSpace*destWidth);
464	bandGreen->RasterIO(GF_Read,
465	windowX,_numValuesY-(windowY+windowHeight),
466	windowWidth,windowHeight,
467	(void*)(tempImage+1),destWidth,destHeight,
468	targetGDALType,pixelSpace,pixelSpace*destWidth);
469	bandBlue->RasterIO(GF_Read,
470	windowX,_numValuesY-(windowY+windowHeight),
471	windowWidth,windowHeight,
472	(void*)(tempImage+2),destWidth,destHeight,
473	targetGDALType,pixelSpace,pixelSpace*destWidth);
474	}
475
476	else if( hasColorTable )
477	{
478	// Pseudocolored image. Convert 1 band + color table to 24bit RGB.
479
480	GDALRasterBand *band;
481	GDALColorTable *ct;
482	int i;
483
484
485	band = _gdalDataSet->GetRasterBand(1);
486
487
488	band->RasterIO(GF_Read,
489	windowX,_numValuesY-(windowY+windowHeight),
490	windowWidth,windowHeight,
491	(void*)(tempImage+0),destWidth,destHeight,
492	targetGDALType,pixelSpace,pixelSpace*destWidth);
493
494
495	ct = band->GetColorTable();
496
497
498	for( i = 0; i < destWidth * destHeight; i++ )
499	{
500	GDALColorEntry sEntry;
501
502
503	// default to greyscale equilvelent.
504	sEntry.c1 = tempImage[i*3];
505	sEntry.c2 = tempImage[i*3];
506	sEntry.c3 = tempImage[i*3];
507
508
509	ct->GetColorEntryAsRGB( tempImage[i*3], &sEntry );
510
511
512	// Apply RGB back over destination image.
513	tempImage[i*3 + 0] = sEntry.c1;
514	tempImage[i*3 + 1] = sEntry.c2;
515	tempImage[i*3 + 2] = sEntry.c3;
516	}
517	}
518
519
520	else if (hasGreyScale)
521	{
522	// Greyscale image. Convert 1 band to 24bit RGB.
523	GDALRasterBand *band;
524
525
526	band = _gdalDataSet->GetRasterBand(1);
527
528
529	band->RasterIO(GF_Read,
530	windowX,_numValuesY-(windowY+windowHeight),
531	windowWidth,windowHeight,
532	(void*)(tempImage+0),destWidth,destHeight,
533	targetGDALType,pixelSpace,pixelSpace*destWidth);
534	band->RasterIO(GF_Read,
535	windowX,_numValuesY-(windowY+windowHeight),
536	windowWidth,windowHeight,
537	(void*)(tempImage+1),destWidth,destHeight,
538	targetGDALType,pixelSpace,pixelSpace*destWidth);
539	band->RasterIO(GF_Read,
540	windowX,_numValuesY-(windowY+windowHeight),
541	windowWidth,windowHeight,
542	(void*)(tempImage+2),destWidth,destHeight,
543	targetGDALType,pixelSpace,pixelSpace*destWidth);
544	}
545
546
547	// now copy into destination image
548	unsigned char* sourceRowPtr = tempImage;
549	unsigned int sourceRowDelta = pixelSpace*destWidth;
550	unsigned char* destinationRowPtr = imageData;
551	unsigned int destinationRowDelta = lineSpace;
552
553	for(int row=0;
554	row<destHeight;
555	++row, sourceRowPtr+=sourceRowDelta, destinationRowPtr+=destinationRowDelta)
556	{
557	unsigned char* sourceColumnPtr = sourceRowPtr;
558	unsigned char* destinationColumnPtr = destinationRowPtr;
559
560	for(int col=0;
561	col<destWidth;
562	++col, sourceColumnPtr+=pixelSpace, destinationColumnPtr+=pixelSpace)
563	{
564	#if 0
565	unsigned int sourceTotal = sourceColumnPtr[0]+sourceColumnPtr[1]+sourceColumnPtr[2];
566	unsigned int destinationTotal = destinationColumnPtr[0]+destinationColumnPtr[1]+destinationColumnPtr[2];
567
568	if (sourceTotal>destinationTotal)
569	{
570	// copy pixel across
571	destinationColumnPtr[0] = sourceColumnPtr[0];
572	destinationColumnPtr[1] = sourceColumnPtr[1];
573	destinationColumnPtr[2] = sourceColumnPtr[2];
574	}
575	#else
576	if (sourceColumnPtr[0]!=0 \|\| sourceColumnPtr[1]!=0 \|\| sourceColumnPtr[2]!=0)
577	{
578	// copy pixel across
579	destinationColumnPtr[0] = sourceColumnPtr[0];
580	destinationColumnPtr[1] = sourceColumnPtr[1];
581	destinationColumnPtr[2] = sourceColumnPtr[2];
582	}
583	#endif
584	}
585	}
586
587	delete [] tempImage;
588
589	}
590	else
591	{
592	osg::notify(osg::INFO)<<"Warning image not read as Red, Blue and Green bands not present"<<std::endl;
593	}
594
595
596	}
597	}
598
599	void DataSet::SourceData::readHeightField(DestinationData& destination)
600	{
601	osg::notify(osg::INFO)<<"In DataSet::SourceData::readHeightField"<<std::endl;
602	if (destination._heightField.valid())
603	{
604	osg::notify(osg::INFO)<<"Reading height field"<<std::endl;
605
606	osg::BoundingBox s_bb = getExtents(destination._cs.get());
607
608	osg::BoundingBox d_bb = destination._extents;
609
610	osg::BoundingBox intersect_bb(s_bb.intersect(d_bb));
611
612	if (!intersect_bb.valid())
613	{
614	osg::notify(osg::INFO)<<"Reading height field but it does not intesection destination - ignoring"<<std::endl;
615	return;
616	}
617
618	// compute dimensions to read from.
619	int windowX = osg::maximum((int)floorf((float)_numValuesX*(intersect_bb.xMin()-s_bb.xMin())/(s_bb.xMax()-s_bb.xMin())),0);
620	int windowY = osg::maximum((int)floorf((float)_numValuesY*(intersect_bb.yMin()-s_bb.yMin())/(s_bb.yMax()-s_bb.yMin())),0);
621	int windowWidth = osg::minimum((int)ceilf((float)_numValuesX(intersect_bb.xMax()-s_bb.xMin())/(s_bb.xMax()-s_bb.xMin())),(int*)_numValuesX)-windowX;
622	int windowHeight = osg::minimum((int)ceilf((float)_numValuesY(intersect_bb.yMax()-s_bb.yMin())/(s_bb.yMax()-s_bb.yMin())),(int*)_numValuesY)-windowY;
623
624	int destX = osg::maximum((int)floorf((float)destination._heightField->getNumColumns()*(intersect_bb.xMin()-d_bb.xMin())/(d_bb.xMax()-d_bb.xMin())),0);
625	int destY = osg::maximum((int)floorf((float)destination._heightField->getNumRows()*(intersect_bb.yMin()-d_bb.yMin())/(d_bb.yMax()-d_bb.yMin())),0);
626	int destWidth = osg::minimum((int)ceilf((float)destination._heightField->getNumColumns()(intersect_bb.xMax()-d_bb.xMin())/(d_bb.xMax()-d_bb.xMin())),(int*)destination._heightField->getNumColumns())-destX;
627	int destHeight = osg::minimum((int)ceilf((float)destination._heightField->getNumRows()(intersect_bb.yMax()-d_bb.yMin())/(d_bb.yMax()-d_bb.yMin())),(int*)destination._heightField->getNumRows())-destY;
628
629	osg::notify(osg::INFO)<<" copying from "<<windowX<<"\t"<<windowY<<"\t"<<windowWidth<<"\t"<<windowHeight<<std::endl;
630	osg::notify(osg::INFO)<<" to "<<destX<<"\t"<<destY<<"\t"<<destWidth<<"\t"<<destHeight<<std::endl;
631
632
633	// which band do we want to read from...
634	int numBands = _gdalDataSet->GetRasterCount();
635	GDALRasterBand* bandGray = 0;
636	GDALRasterBand* bandRed = 0;
637	GDALRasterBand* bandGreen = 0;
638	GDALRasterBand* bandBlue = 0;
639	GDALRasterBand* bandAlpha = 0;
640
641	for(int b=1;b<=numBands;++b)
642	{
643	GDALRasterBand* band = _gdalDataSet->GetRasterBand(b);
644	if (band->GetColorInterpretation()==GCI_GrayIndex) bandGray = band;
645	else if (band->GetColorInterpretation()==GCI_RedBand) bandRed = band;
646	else if (band->GetColorInterpretation()==GCI_GreenBand) bandGreen = band;
647	else if (band->GetColorInterpretation()==GCI_BlueBand) bandBlue = band;
648	else if (band->GetColorInterpretation()==GCI_AlphaBand) bandAlpha = band;
649	else bandGray = band;
650	}
651
652
653	GDALRasterBand* bandSelected = 0;
654	if (!bandSelected && bandGray) bandSelected = bandGray;
655	else if (!bandSelected && bandAlpha) bandSelected = bandAlpha;
656	else if (!bandSelected && bandRed) bandSelected = bandRed;
657	else if (!bandSelected && bandGreen) bandSelected = bandGreen;
658	else if (!bandSelected && bandBlue) bandSelected = bandBlue;
659
660	if (bandSelected)
661	{
662
663	bool xyInDegrees = false;
664
665	CoordinateSystemType cst = getCoordinateSystemType(_cs.get());
666	if (cst==GEOGRAPHIC)
667	{
668	xyInDegrees = true;
669	}
670
671
672	if (bandSelected->GetUnitType()) osg::notify(osg::INFO) << "bandSelected->GetUnitType()=" << bandSelected->GetUnitType()<<std::endl;
673	else osg::notify(osg::INFO) << "bandSelected->GetUnitType()= null" <<std::endl;
674
675
676	int success = 0;
677	float noDataValue = bandSelected->GetNoDataValue(&success);
678	if (success)
679	{
680	osg::notify(osg::INFO)<<"We have NoDataValue = "<<noDataValue<<std::endl;
681	}
682	else
683	{
684	osg::notify(osg::INFO)<<"We have no NoDataValue"<<std::endl;
685	noDataValue = 0.0f;
686	}
687
688	float offset = bandSelected->GetOffset(&success);
689	if (success)
690	{
691	osg::notify(osg::INFO)<<"We have Offset = "<<offset<<std::endl;
692	}
693	else
694	{
695	osg::notify(osg::INFO)<<"We have no Offset"<<std::endl;
696	offset = 0.0f;
697	}
698
699	float scale = bandSelected->GetScale(&success);
700	if (success)
701	{
702	osg::notify(osg::INFO)<<"We have Scale = "<<scale<<std::endl;
703	}
704	else
705	{
706	scale = destination._dataSet->getVerticalScale();
707	osg::notify(osg::INFO)<<"We have no Scale from file so use DataSet vertical scale of "<<scale<<std::endl;
708	//scale = (xyInDegrees /&& !destination._dataSet->getConvertFromGeographicToGeocentric()/) ? 1.0f/111319.0f : 1.0f;
709
710	}
711
712	osg::notify(osg::INFO)<<"***** getLinearUnits = "<<getLinearUnits(_cs.get())<<std::endl;
713
714	// read data into temporary array
715	float* heightData = new float [ destWidth*destHeight ];
716
717	// raad the data.
718	osg::HeightField* hf = destination._heightField.get();
719
720	osg::notify(osg::INFO)<<"reading height field"<<std::endl;
721	//bandSelected->RasterIO(GF_Read,windowX,_numValuesY-(windowY+windowHeight),windowWidth,windowHeight,floatdata,destWidth,destHeight,GDT_Float32,numBytesPerZvalue,lineSpace);
722	bandSelected->RasterIO(GF_Read,windowX,_numValuesY-(windowY+windowHeight),windowWidth,windowHeight,heightData,destWidth,destHeight,GDT_Float32,0,0);
723
724	osg::notify(osg::INFO)<<" scaling height field"<<std::endl;
725
726	float noDataValueFill = 0.0f;
727	bool ignoreNoDataValue = true;
728
729	float* heightPtr = heightData;
730
731	for(int r=destY+destHeight-1;r>=destY;--r)
732	{
733	for(int c=destX;c<destX+destWidth;++c)
734	{
735	float h = *heightPtr++;
736	if (h!=noDataValue) hf->setHeight(c,r,offset + h*scale);
737	else if (!ignoreNoDataValue) hf->setHeight(c,r,noDataValueFill);
738
739	h = hf->getHeight(c,r);
740	}
741	}
742
743	delete [] heightData;
744
745	}
746
747	}
748	}
749
750	void DataSet::SourceData::readModels(DestinationData& destination)
751	{
752	if (_model.valid())
753	{
754	osg::notify(osg::INFO)<<"Raading model"<<std::endl;
755	destination._models.push_back(_model);
756	}
757	}
758
759	void DataSet::Source::setSortValueFromSourceDataResolution()
760	{
761	if (_sourceData.valid())
762	{
763	double dx = (_sourceData->_extents.xMax()-_sourceData->_extents.xMin())/(double)(_sourceData->_numValuesX-1);
764	double dy = (_sourceData->_extents.yMax()-_sourceData->_extents.yMin())/(double)(_sourceData->_numValuesY-1);
765
766	setSortValue( sqrt( dxdx + dydy ) );
767	}
768	}
769
770	void DataSet::Source::loadSourceData()
771	{
772	osg::notify(osg::INFO)<<"DataSet::Source::loadSourceData() "<<_filename<<std::endl;
773
774	_sourceData = SourceData::readData(this);
775
776	assignCoordinateSystemAndGeoTransformAccordingToParameterPolicy();
777	}
778
779	void DataSet::Source::assignCoordinateSystemAndGeoTransformAccordingToParameterPolicy()
780	{
781	if (getCoordinateSystemPolicy()==PREFER_CONFIG_SETTINGS)
782	{
783	_sourceData->_cs = _cs;
784
785	osg::notify(osg::INFO)<<"assigning CS from Source to Data."<<std::endl;
786
787	}
788	else
789	{
790	_cs = _sourceData->_cs;
791	osg::notify(osg::INFO)<<"assigning CS from Data to Source."<<std::endl;
792	}
793
794	if (getGeoTransformPolicy()==PREFER_CONFIG_SETTINGS)
795	{
796	_sourceData->_geoTransform = _geoTransform;
797
798	osg::notify(osg::INFO)<<"assigning GeoTransform from Source to Data."<<_geoTransform<<std::endl;
799
800	}
801	else if (getGeoTransformPolicy()==PREFER_CONFIG_SETTINGS_BUT_SCALE_BY_FILE_RESOLUTION)
802	{
803
804	// scale the x and y axis.
805	double div_x = 1.0/(double)(_sourceData->_numValuesX - 1);
806	double div_y = 1.0/(double)(_sourceData->_numValuesY - 1);
807
808	_geoTransform(0,0) *= div_x;
809	_geoTransform(1,0) *= div_x;
810	_geoTransform(2,0) *= div_x;
811
812	_geoTransform(0,1) *= div_y;
813	_geoTransform(1,1) *= div_y;
814	_geoTransform(2,1) *= div_y;
815
816	_sourceData->_geoTransform = _geoTransform;
817
818	osg::notify(osg::INFO)<<"assigning GeoTransform from Source to Data."<<_geoTransform<<std::endl;
819
820	}
821	else
822	{
823	_geoTransform = _sourceData->_geoTransform;
824	osg::notify(osg::INFO)<<"assigning GeoTransform from Data to Source."<<_geoTransform<<std::endl;
825	}
826
827	_sourceData->computeExtents();
828
829	_extents = _sourceData->_extents;
830	}
831
832	bool DataSet::Source::needReproject(const osg::CoordinateSystemNode* cs) const
833	{
834	return needReproject(cs,0.0,0.0);
835	}
836
837	bool DataSet::Source::needReproject(const osg::CoordinateSystemNode* cs, double minResolution, double maxResolution) const
838	{
839	if (!_sourceData) return false;
840
841	// handle modles by using a matrix transform only.
842	if (_type==MODEL) return false;
843
844	// always need to reproject imagery with GCP's.
845	if (_sourceData->_hasGCPs)
846	{
847	osg::notify(osg::INFO)<<"Need to to reproject due to presence of GCP's"<<std::endl;
848	return true;
849	}
850
851	if (!areCoordinateSystemEquivalent(_cs.get(),cs))
852	{
853	osg::notify(osg::INFO)<<"Need to do reproject !areCoordinateSystemEquivalent(_cs.get(),cs)"<<std::endl;
854
855	return true;
856	}
857
858	if (minResolution==0.0 && maxResolution==0.0) return false;
859
860	// now check resolutions.
861	const osg::Matrixd& m = _sourceData->_geoTransform;
862	double currentResolution = sqrt(osg::square(m(0,0))+osg::square(m(1,0))+
863	osg::square(m(0,1))+osg::square(m(1,1)));
864
865	if (currentResolution<minResolution) return true;
866	if (currentResolution>maxResolution) return true;
867
868	return false;
869	}
870
871	DataSet::Source* DataSet::Source::doReproject(const std::string& filename, osg::CoordinateSystemNode* cs, double targetResolution) const
872	{
873	// return nothing when repoject is inappropriate.
874	if (!_sourceData) return 0;
875	if (_type==MODEL) return 0;
876
877	osg::notify(osg::INFO)<<"reprojecting to file "<<filename<<std::endl;
878
879	GDALDriverH hDriver = GDALGetDriverByName( "GTiff" );
880
881	if (hDriver == NULL)
882	{
883	osg::notify(osg::INFO)<<"Unable to load driver for "<<"GTiff"<<std::endl;
884	return 0;
885	}
886
887	if (GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, NULL ) == NULL )
888	{
889	osg::notify(osg::INFO)<<"GDAL driver does not support create for "<<osgDB::getFileExtension(filename)<<std::endl;
890	return 0;
891	}
892
893	/ -------------------------------------------------------------------- /
894	/ Create a transformation object from the source to /
895	/ destination coordinate system. /
896	/ -------------------------------------------------------------------- /
897	void *hTransformArg =
898	GDALCreateGenImgProjTransformer( _sourceData->_gdalDataSet,_sourceData->_cs->getCoordinateSystem().c_str(),
899	NULL, cs->getCoordinateSystem().c_str(),
900	TRUE, 0.0, 1 );
901
902	if (!hTransformArg)
903	{
904	osg::notify(osg::INFO)<<" failed to create transformer"<<std::endl;
905	return 0;
906	}
907
908	double adfDstGeoTransform[6];
909	int nPixels=0, nLines=0;
910	if( GDALSuggestedWarpOutput( _sourceData->_gdalDataSet,
911	GDALGenImgProjTransform, hTransformArg,
912	adfDstGeoTransform, &nPixels, &nLines )
913	!= CE_None )
914	{
915	osg::notify(osg::INFO)<<" failed to create warp"<<std::endl;
916	return 0;
917	}
918
919	if (targetResolution>0.0f)
920	{
921	osg::notify(osg::INFO)<<"recomputing the target transform size"<<std::endl;
922
923	double currentResolution = sqrt(osg::square(adfDstGeoTransform[1])+osg::square(adfDstGeoTransform[2])+
924	osg::square(adfDstGeoTransform[4])+osg::square(adfDstGeoTransform[5]));
925
926	osg::notify(osg::INFO)<<" default computed resolution "<<currentResolution<<" nPixels="<<nPixels<<" nLines="<<nLines<<std::endl;
927
928	double extentsPixels = sqrt(osg::square(adfDstGeoTransform[1])+osg::square(adfDstGeoTransform[2]))(double*)(nPixels-1);
929	double extentsLines = sqrt(osg::square(adfDstGeoTransform[4])+osg::square(adfDstGeoTransform[5]))(double*)(nLines-1);
930
931	double ratio = targetResolution/currentResolution;
932	adfDstGeoTransform[1] *= ratio;
933	adfDstGeoTransform[2] *= ratio;
934	adfDstGeoTransform[4] *= ratio;
935	adfDstGeoTransform[5] *= ratio;
936
937	osg::notify(osg::INFO)<<" extentsPixels="<<extentsPixels<<std::endl;
938	osg::notify(osg::INFO)<<" extentsLines="<<extentsLines<<std::endl;
939	osg::notify(osg::INFO)<<" targetResolution="<<targetResolution<<std::endl;
940
941	nPixels = (int)ceil(extentsPixels/sqrt(osg::square(adfDstGeoTransform[1])+osg::square(adfDstGeoTransform[2])))+1;
942	nLines = (int)ceil(extentsLines/sqrt(osg::square(adfDstGeoTransform[4])+osg::square(adfDstGeoTransform[5])))+1;
943
944	osg::notify(osg::INFO)<<" target computed resolution "<<targetResolution<<" nPixels="<<nPixels<<" nLines="<<nLines<<std::endl;
945
946	}
947
948
949	GDALDestroyGenImgProjTransformer( hTransformArg );
950
951	GDALDataType eDT = GDALGetRasterDataType(GDALGetRasterBand(_sourceData->_gdalDataSet,1));
952
953
954	/ --------------------------------------------------------------------- /
955	/ Create the file /
956	/ --------------------------------------------------------------------- /
957
958	GDALDatasetH hDstDS = GDALCreate( hDriver, filename.c_str(), nPixels, nLines,
959	GDALGetRasterCount(_sourceData->_gdalDataSet), eDT,
960	0 );
961
962	if( hDstDS == NULL )
963	return NULL;
964
965	/ -------------------------------------------------------------------- /
966	/ Write out the projection definition. /
967	/ -------------------------------------------------------------------- /
968	GDALSetProjection( hDstDS, cs->getCoordinateSystem().c_str() );
969	GDALSetGeoTransform( hDstDS, adfDstGeoTransform );
970
971
972	// Set up the transformer along with the new datasets.
973
974	hTransformArg =
975	GDALCreateGenImgProjTransformer( _sourceData->_gdalDataSet,_sourceData->_cs->getCoordinateSystem().c_str(),
976	hDstDS, cs->getCoordinateSystem().c_str(),
977	TRUE, 0.0, 1 );
978
979	GDALTransformerFunc pfnTransformer = GDALGenImgProjTransform;
980
981
982	osg::notify(osg::INFO)<<"Setting projection "<<cs->getCoordinateSystem()<<std::endl;
983
984	/ -------------------------------------------------------------------- /
985	/ Copy the color table, if required. /
986	/ -------------------------------------------------------------------- /
987	GDALColorTableH hCT;
988
989	hCT = GDALGetRasterColorTable( GDALGetRasterBand(_sourceData->_gdalDataSet,1) );
990	if( hCT != NULL )
991	GDALSetRasterColorTable( GDALGetRasterBand(hDstDS,1), hCT );
992
993	/ -------------------------------------------------------------------- /
994	/ Setup warp options. /
995	/ -------------------------------------------------------------------- /
996	GDALWarpOptions *psWO = GDALCreateWarpOptions();
997
998	psWO->hSrcDS = _sourceData->_gdalDataSet;
999	psWO->hDstDS = hDstDS;
1000
1001	psWO->pfnTransformer = pfnTransformer;
1002	psWO->pTransformerArg = hTransformArg;
1003
1004	psWO->pfnProgress = GDALTermProgress;
1005
1006	/ -------------------------------------------------------------------- /
1007	/ Setup band mapping. /
1008	/ -------------------------------------------------------------------- /
1009	psWO->nBandCount = GDALGetRasterCount(_sourceData->_gdalDataSet);
1010	psWO->panSrcBands = (int ) CPLMalloc(psWO->nBandCountsizeof(int));
1011	psWO->panDstBands = (int ) CPLMalloc(psWO->nBandCountsizeof(int));
1012
1013	int i;
1014	for(i = 0; i < psWO->nBandCount; i++ )
1015	{
1016	psWO->panSrcBands[i] = i+1;
1017	psWO->panDstBands[i] = i+1;
1018	}
1019
1020
1021	/ -------------------------------------------------------------------- /
1022	/ Setup no datavalue /
1023	/ -----------------------------------------------------`--------------- /
1024
1025
1026	// check to see if no value values exist in source datasets.
1027	int numNoDataValues = 0;
1028	for(i = 0; i < _sourceData->_gdalDataSet->GetRasterCount(); i++ )
1029	{
1030	int success = 0;
1031	GDALRasterBand* band = _sourceData->_gdalDataSet->GetRasterBand(i+1);
1032	band->GetNoDataValue(&success);
1033	if (success) ++numNoDataValues;
1034	}
1035
1036	if (numNoDataValues)
1037	{
1038	// no data values exist, so populate the no data arrays.
1039
1040	psWO->padfSrcNoDataReal = (double) CPLMalloc(psWO->nBandCountsizeof(double));
1041	psWO->padfSrcNoDataImag = (double) CPLMalloc(psWO->nBandCountsizeof(double));
1042
1043	psWO->padfDstNoDataReal = (double) CPLMalloc(psWO->nBandCountsizeof(double));
1044	psWO->padfDstNoDataImag = (double) CPLMalloc(psWO->nBandCountsizeof(double));
1045
1046	for(i = 0; i < _sourceData->_gdalDataSet->GetRasterCount(); i++ )
1047	{
1048	int success = 0;
1049	GDALRasterBand* band = _sourceData->_gdalDataSet->GetRasterBand(i+1);
1050	double noDataValue = band->GetNoDataValue(&success);
1051	//double new_noDataValue = noDataValue;
1052	double new_noDataValue = 0;
1053	if (success)
1054	{
1055	osg::notify(osg::INFO)<<"\tassinging no data value "<<noDataValue<<" to band "<<i+1<<std::endl;
1056
1057	psWO->padfSrcNoDataReal[i] = noDataValue;
1058	psWO->padfSrcNoDataImag[i] = 0.0;
1059	psWO->padfDstNoDataReal[i] = new_noDataValue;
1060	psWO->padfDstNoDataImag[i] = 0.0;
1061
1062	GDALRasterBandH band = GDALGetRasterBand(hDstDS,i+1);
1063	GDALSetRasterNoDataValue( band, new_noDataValue);
1064	}
1065	}
1066
1067	psWO->papszWarpOptions = (char)CPLMalloc(2sizeof(char**));
1068	psWO->papszWarpOptions[0] = strdup("INIT_DEST=NO_DATA");
1069	psWO->papszWarpOptions[1] = 0;
1070
1071	}
1072	else
1073	{
1074
1075	psWO->padfSrcNoDataReal = (double) CPLMalloc(psWO->nBandCountsizeof(double));
1076	psWO->padfSrcNoDataImag = (double) CPLMalloc(psWO->nBandCountsizeof(double));
1077
1078	psWO->padfDstNoDataReal = (double) CPLMalloc(psWO->nBandCountsizeof(double));
1079	psWO->padfDstNoDataImag = (double) CPLMalloc(psWO->nBandCountsizeof(double));
1080
1081	for(i = 0; i < _sourceData->_gdalDataSet->GetRasterCount(); i++ )
1082	{
1083	int success = 0;
1084	GDALRasterBand* band = _sourceData->_gdalDataSet->GetRasterBand(i+1);
1085	double noDataValue = band->GetNoDataValue(&success);
1086	double new_noDataValue = 0.0;
1087	if (success)
1088	{
1089	osg::notify(osg::INFO)<<"\tassinging no data value "<<noDataValue<<" to band "<<i+1<<std::endl;
1090
1091	psWO->padfSrcNoDataReal[i] = noDataValue;
1092	psWO->padfSrcNoDataImag[i] = 0.0;
1093	psWO->padfDstNoDataReal[i] = noDataValue;
1094	psWO->padfDstNoDataImag[i] = 0.0;
1095
1096	GDALRasterBandH band = GDALGetRasterBand(hDstDS,i+1);
1097	GDALSetRasterNoDataValue( band, new_noDataValue);
1098	}
1099	else
1100	{
1101	psWO->padfSrcNoDataReal[i] = 0.0;
1102	psWO->padfSrcNoDataImag[i] = 0.0;
1103	psWO->padfDstNoDataReal[i] = new_noDataValue;
1104	psWO->padfDstNoDataImag[i] = 0.0;
1105
1106	GDALRasterBandH band = GDALGetRasterBand(hDstDS,i+1);
1107	GDALSetRasterNoDataValue( band, new_noDataValue);
1108	}
1109	}
1110
1111	psWO->papszWarpOptions = (char)CPLMalloc(2sizeof(char**));
1112	psWO->papszWarpOptions[0] = strdup("INIT_DEST=NO_DATA");
1113	psWO->papszWarpOptions[1] = 0;
1114
1115	}
1116
1117
1118	/ -------------------------------------------------------------------- /
1119	/ Initialize and execute the warp. /
1120	/ -------------------------------------------------------------------- /
1121	GDALWarpOperation oWO;
1122
1123	if( oWO.Initialize( psWO ) == CE_None )
1124	{
1125	oWO.ChunkAndWarpImage( 0, 0,
1126	GDALGetRasterXSize( hDstDS ),
1127	GDALGetRasterYSize( hDstDS ) );
1128	}
1129
1130	osg::notify(osg::INFO)<<"new projection is "<<GDALGetProjectionRef(hDstDS)<<std::endl;
1131
1132	/ -------------------------------------------------------------------- /
1133	/ Cleanup. /
1134	/ -------------------------------------------------------------------- /
1135	GDALDestroyGenImgProjTransformer( hTransformArg );
1136
1137	#if 0
1138	int anOverviewList[4] = { 2, 4, 8, 16 };
1139	GDALBuildOverviews( hDstDS, "AVERAGE", 4, anOverviewList, 0, NULL,
1140	GDALTermProgress/GDALDummyProgress/, NULL );
1141	#endif
1142
1143	GDALClose( hDstDS );
1144
1145	Source* newSource = new Source;
1146	newSource->_type = _type;
1147	newSource->_filename = filename;
1148	newSource->_temporaryFile = true;
1149	newSource->_cs = cs;
1150	newSource->_numValuesX = nPixels;
1151	newSource->_numValuesY = nLines;
1152	newSource->_geoTransform.set( adfDstGeoTransform[1], adfDstGeoTransform[4], 0.0, 0.0,
1153	adfDstGeoTransform[2], adfDstGeoTransform[5], 0.0, 0.0,
1154	0.0, 0.0, 1.0, 0.0,
1155	adfDstGeoTransform[0], adfDstGeoTransform[3], 0.0, 1.0);
1156
1157	newSource->computeExtents();
1158
1159	// reload the newly created file.
1160	newSource->loadSourceData();
1161
1162	return newSource;
1163	}
1164
1165	void DataSet::Source::consolodateRequiredResolutions()
1166	{
1167	if (_requiredResolutions.size()<=1) return;
1168
1169	ResolutionList consolodated;
1170
1171	ResolutionList::iterator itr = _requiredResolutions.begin();
1172
1173	double minResX = itr->_resX;
1174	double minResY = itr->_resY;
1175	double maxResX = itr->_resX;
1176	double maxResY = itr->_resY;
1177	++itr;
1178	for(;itr!=_requiredResolutions.end();++itr)
1179	{
1180	minResX = osg::minimum(minResX,itr->_resX);
1181	minResY = osg::minimum(minResY,itr->_resY);
1182	maxResX = osg::maximum(maxResX,itr->_resX);
1183	maxResY = osg::maximum(maxResY,itr->_resY);
1184	}
1185
1186	double currResX = minResX;
1187	double currResY = minResY;
1188	while (currResX<=maxResX && currResY<=maxResY)
1189	{
1190	consolodated.push_back(ResolutionPair(currResX,currResY));
1191	currResX *= 2.0f;
1192	currResY *= 2.0f;
1193	}
1194
1195
1196	_requiredResolutions.swap(consolodated);
1197	}
1198
1199	void DataSet::Source::buildOverviews()
1200	{
1201	return;
1202
1203	if (_sourceData.valid() && _sourceData->_gdalDataSet )
1204	{
1205
1206	int anOverviewList[4] = { 2, 4, 8, 16 };
1207	GDALBuildOverviews( _sourceData->_gdalDataSet, "AVERAGE", 4, anOverviewList, 0, NULL,
1208	GDALTermProgress/GDALDummyProgress/, NULL );
1209
1210	}
1211	}
1212
1213
1214	///////////////////////////////////////////////////////////////////////////////////////////////////////
1215
1216
1217	DataSet::DestinationTile::DestinationTile():
1218	_dataSet(0),
1219	_level(0),
1220	_tileX(0),
1221	_tileY(0),
1222	_maxSourceLevel(0),
1223	_imagery_maxNumColumns(4096),
1224	_imagery_maxNumRows(4096),
1225	_imagery_maxSourceResolutionX(0.0f),
1226	_imagery_maxSourceResolutionY(0.0f),
1227	_terrain_maxNumColumns(1024),
1228	_terrain_maxNumRows(1024),
1229	_terrain_maxSourceResolutionX(0.0f),
1230	_terrain_maxSourceResolutionY(0.0f),
1231	_complete(false)
1232	{
1233	for(int i=0;i<NUMBER_OF_POSITIONS;++i)
1234	{
1235	_neighbour[i]=0;
1236	_equalized[i]=false;
1237	}
1238	}
1239
1240
1241	void DataSet::DestinationTile::computeMaximumSourceResolution(CompositeSource* sourceGraph)
1242	{
1243	for(CompositeSource::source_iterator itr(sourceGraph);itr.valid();++itr)
1244	{
1245	Source* source = itr->get();
1246	if (!source \|\| source->getMaxLevel()<_level)
1247	{
1248	// skip the contribution of this source since this destination tile exceeds its contribution level.
1249	continue;
1250	}
1251
1252	SourceData* data = source->getSourceData();
1253	if (data && source->getType()!=Source::MODEL)
1254	{
1255
1256	SpatialProperties sp = data->computeSpatialProperties(_cs.get());
1257
1258	if (!sp._extents.intersects(_extents))
1259	{
1260	// skip this source since it doesn't overlap this tile.
1261	continue;
1262	}
1263
1264
1265	if (sp._numValuesX!=0 && sp._numValuesY!=0)
1266	{
1267	_maxSourceLevel = osg::maximum((*itr)->getMaxLevel(),_maxSourceLevel);
1268
1269	float sourceResolutionX = (sp._extents.xMax()-sp._extents.xMin())/(float)sp._numValuesX;
1270	float sourceResolutionY = (sp._extents.yMax()-sp._extents.yMin())/(float)sp._numValuesY;
1271
1272	switch((*itr)->getType())
1273	{
1274	case(Source::IMAGE):
1275	if (_imagery_maxSourceResolutionX==0.0f) _imagery_maxSourceResolutionX=sourceResolutionX;
1276	else _imagery_maxSourceResolutionX=osg::minimum(_imagery_maxSourceResolutionX,sourceResolutionX);
1277	if (_imagery_maxSourceResolutionY==0.0f) _imagery_maxSourceResolutionY=sourceResolutionY;
1278	else _imagery_maxSourceResolutionY=osg::minimum(_imagery_maxSourceResolutionY,sourceResolutionY);
1279	break;
1280	case(Source::HEIGHT_FIELD):
1281	if (_terrain_maxSourceResolutionX==0.0f) _terrain_maxSourceResolutionX=sourceResolutionX;
1282	else _terrain_maxSourceResolutionX=osg::minimum(_terrain_maxSourceResolutionX,sourceResolutionX);
1283	if (_terrain_maxSourceResolutionY==0.0f) _terrain_maxSourceResolutionY=sourceResolutionY;
1284	else _terrain_maxSourceResolutionY=osg::minimum(_terrain_maxSourceResolutionY,sourceResolutionY);
1285	break;
1286	default:
1287	break;
1288	}
1289	}
1290	}
1291	}
1292	}
1293
1294
1295	bool DataSet::DestinationTile::computeImageResolution(unsigned int& numColumns, unsigned int& numRows, double& resX, double& resY)
1296	{
1297	if (_imagery_maxSourceResolutionX!=0.0f && _imagery_maxSourceResolutionY!=0.0f &&
1298	_imagery_maxNumColumns!=0 && _imagery_maxNumRows!=0)
1299	{
1300
1301	unsigned int numColumnsAtFullRes = 1+(unsigned int)ceilf((_extents.xMax()-_extents.xMin())/_imagery_maxSourceResolutionX);
1302	unsigned int numRowsAtFullRes = 1+(unsigned int)ceilf((_extents.yMax()-_extents.yMin())/_imagery_maxSourceResolutionY);
1303	unsigned int numColumnsRequired = osg::minimum(_imagery_maxNumColumns,numColumnsAtFullRes);
1304	unsigned int numRowsRequired = osg::minimum(_imagery_maxNumRows,numRowsAtFullRes);
1305
1306	numColumns = 1;
1307	numRows = 1;
1308	// round to nearest power of two above or equal to the required resolution
1309	while (numColumns<numColumnsRequired) numColumns *= 2;
1310	while (numRows<numRowsRequired) numRows *= 2;
1311
1312	resX = (_extents.xMax()-_extents.xMin())/(double)(numColumns-1);
1313	resY = (_extents.yMax()-_extents.yMin())/(double)(numRows-1);
1314
1315	return true;
1316	}
1317	return false;
1318	}
1319
1320	bool DataSet::DestinationTile::computeTerrainResolution(unsigned int& numColumns, unsigned int& numRows, double& resX, double& resY)
1321	{
1322	if (_terrain_maxSourceResolutionX!=0.0f && _terrain_maxSourceResolutionY!=0.0f &&
1323	_terrain_maxNumColumns!=0 && _terrain_maxNumRows!=0)
1324	{
1325	unsigned int numColumnsAtFullRes = 1+(unsigned int)ceilf((_extents.xMax()-_extents.xMin())/_terrain_maxSourceResolutionX);
1326	unsigned int numRowsAtFullRes = 1+(unsigned int)ceilf((_extents.yMax()-_extents.yMin())/_terrain_maxSourceResolutionY);
1327	numColumns = osg::minimum(_terrain_maxNumColumns,numColumnsAtFullRes);
1328	numRows = osg::minimum(_terrain_maxNumRows,numRowsAtFullRes);
1329
1330	resX = (_extents.xMax()-_extents.xMin())/(double)(numColumns-1);
1331	resY = (_extents.yMax()-_extents.yMin())/(double)(numRows-1);
1332
1333	return true;
1334	}
1335	return false;
1336	}
1337
1338
1339	void DataSet::DestinationTile::allocate()
1340	{
1341	unsigned int texture_numColumns, texture_numRows;
1342	double texture_dx, texture_dy;
1343	if (computeImageResolution(texture_numColumns,texture_numRows,texture_dx,texture_dy))
1344	{
1345
1346	_imagery = new DestinationData(_dataSet);
1347	_imagery->_cs = _cs;
1348	_imagery->_extents = _extents;
1349	_imagery->_geoTransform.set(texture_dx, 0.0, 0.0,0.0,
1350	0.0, -texture_dy, 0.0,0.0,
1351	0.0, 0.0, 1.0,1.0,
1352	_extents.xMin(), _extents.yMax(), 0.0,1.0);
1353	_imagery->_image = new osg::Image;
1354	_imagery->_image->allocateImage(texture_numColumns,texture_numRows,1,GL_RGB,GL_UNSIGNED_BYTE);
1355	unsigned char* data = _imagery->_image->data();
1356	unsigned int totalSize = _imagery->_image->getTotalSizeInBytesIncludingMipmaps();
1357	for(unsigned int i=0;i<totalSize;++i)
1358	{
1359	*(data++) = 0;
1360	}
1361	}
1362
1363
1364	unsigned int dem_numColumns, dem_numRows;
1365	double dem_dx, dem_dy;
1366	if (computeTerrainResolution(dem_numColumns,dem_numRows,dem_dx,dem_dy))
1367	{
1368	_terrain = new DestinationData(_dataSet);
1369	_terrain->_cs = _cs;
1370	_terrain->_extents = _extents;
1371	_terrain->_geoTransform.set(dem_dx, 0.0, 0.0,0.0,
1372	0.0, -dem_dy, 0.0,0.0,
1373	0.0, 0.0, 1.0,1.0,
1374	_extents.xMin(), _extents.yMax(), 0.0,1.0);
1375	_terrain->_heightField = new osg::HeightField;
1376	_terrain->_heightField->allocate(dem_numColumns,dem_numRows);
1377	_terrain->_heightField->setOrigin(osg::Vec3(_extents.xMin(),_extents.yMin(),0.0f));
1378	_terrain->_heightField->setXInterval(dem_dx);
1379	_terrain->_heightField->setYInterval(dem_dy);
1380
1381	//float xMax = _terrain->_heightField->getOrigin().x()+_terrain->_heightField->getXInterval()(float)(dem_numColumns-1);*
1382	//osg::notify(osg::INFO)<<"ErrorX = "<<xMax-_extents.xMax()<<std::endl;
1383
1384	//float yMax = _terrain->_heightField->getOrigin().y()+_terrain->_heightField->getYInterval()(float)(dem_numRows-1);*
1385	//osg::notify(osg::INFO)<<"ErrorY = "<<yMax-_extents.yMax()<<std::endl;
1386
1387	}
1388
1389	}
1390
1391	void DataSet::DestinationTile::computeNeighboursFromQuadMap()
1392	{
1393	if (_dataSet)
1394	{
1395	setNeighbours(_dataSet->getTile(_level,_tileX-1,_tileY),_dataSet->getTile(_level,_tileX-1,_tileY-1),
1396	_dataSet->getTile(_level,_tileX,_tileY-1),_dataSet->getTile(_level,_tileX+1,_tileY-1),
1397	_dataSet->getTile(_level,_tileX+1,_tileY),_dataSet->getTile(_level,_tileX+1,_tileY+1),
1398	_dataSet->getTile(_level,_tileX,_tileY+1),_dataSet->getTile(_level,_tileX-1,_tileY+1));
1399	}
1400	}
1401
1402	void DataSet::DestinationTile::setNeighbours(DestinationTile* left, DestinationTile* left_below,
1403	DestinationTile* below, DestinationTile* below_right,
1404	DestinationTile* right, DestinationTile* right_above,
1405	DestinationTile* above, DestinationTile* above_left)
1406	{
1407	_neighbour[LEFT] = left;
1408	_neighbour[LEFT_BELOW] = left_below;
1409	_neighbour[BELOW] = below;
1410	_neighbour[BELOW_RIGHT] = below_right;
1411	_neighbour[RIGHT] = right;
1412	_neighbour[RIGHT_ABOVE] = right_above;
1413	_neighbour[ABOVE] = above;
1414	_neighbour[ABOVE_LEFT] = above_left;
1415
1416
1417	// osg::notify(osg::INFO)<<"LEFT="<<_neighbour[LEFT]<<std::endl;
1418	// osg::notify(osg::INFO)<<"LEFT_BELOW="<<_neighbour[LEFT_BELOW]<<std::endl;
1419	// osg::notify(osg::INFO)<<"BELOW="<<_neighbour[BELOW]<<std::endl;
1420	// osg::notify(osg::INFO)<<"BELOW_RIGHT="<<_neighbour[BELOW_RIGHT]<<std::endl;
1421	// osg::notify(osg::INFO)<<"RIGHT="<<_neighbour[RIGHT]<<std::endl;
1422	// osg::notify(osg::INFO)<<"RIGHT_ABOVE="<<_neighbour[RIGHT_ABOVE]<<std::endl;
1423	// osg::notify(osg::INFO)<<"ABOVE="<<_neighbour[ABOVE]<<std::endl;
1424	// osg::notify(osg::INFO)<<"ABOVE_LEFT="<<_neighbour[ABOVE_LEFT]<<std::endl;
1425
1426	for(int i=0;i<NUMBER_OF_POSITIONS;++i)
1427	{
1428	_equalized[i]=false;
1429	}
1430	}
1431
1432	void DataSet::DestinationTile::checkNeighbouringTiles()
1433	{
1434	for(int i=0;i<NUMBER_OF_POSITIONS;++i)
1435	{
1436	if (_neighbour[i] && _neighbour[i]->_neighbour[(i+4)%NUMBER_OF_POSITIONS]!=this)
1437	{
1438	osg::notify(osg::INFO)<<"Error:: Tile "<<this<<"'s _neighbour["<<i<<"] does not point back to it."<<std::endl;
1439	}
1440	}
1441	}
1442
1443	void DataSet::DestinationTile::equalizeCorner(Position position)
1444	{
1445	// don't need to equalize if already done.
1446	if (_equalized[position]) return;
1447
1448	typedef std::pair<DestinationTile*,Position> TileCornerPair;
1449	typedef std::vector<TileCornerPair> TileCornerList;
1450
1451	TileCornerList cornersToProcess;
1452	DestinationTile* tile=0;
1453
1454	cornersToProcess.push_back(TileCornerPair(this,position));
1455
1456	tile = _neighbour[(position-1)%NUMBER_OF_POSITIONS];
1457	if (tile) cornersToProcess.push_back(TileCornerPair(tile,(Position)((position+2)%NUMBER_OF_POSITIONS)));
1458
1459	tile = _neighbour[(position)%NUMBER_OF_POSITIONS];
1460	if (tile) cornersToProcess.push_back(TileCornerPair(tile,(Position)((position+4)%NUMBER_OF_POSITIONS)));
1461
1462	tile = _neighbour[(position+1)%NUMBER_OF_POSITIONS];
1463	if (tile) cornersToProcess.push_back(TileCornerPair(tile,(Position)((position+6)%NUMBER_OF_POSITIONS)));
1464
1465	// make all these tiles as equalised upfront before we return.
1466	TileCornerList::iterator itr;
1467	for(itr=cornersToProcess.begin();
1468	itr!=cornersToProcess.end();
1469	++itr)
1470	{
1471	TileCornerPair& tcp = *itr;
1472	tcp.first->_equalized[tcp.second] = true;
1473	}
1474
1475	// if there is only one valid corner to process then there is nothing to equalize against so return.
1476	if (cornersToProcess.size()==1) return;
1477
1478	typedef std::pair<osg::Image*,Position> ImageCornerPair;
1479	typedef std::vector<ImageCornerPair> ImageCornerList;
1480
1481	typedef std::pair<osg::HeightField*,Position> HeightFieldCornerPair;
1482	typedef std::vector<HeightFieldCornerPair> HeightFieldCornerList;
1483
1484	ImageCornerList imagesToProcess;
1485	HeightFieldCornerList heightFieldsToProcess;
1486
1487	for(itr=cornersToProcess.begin();
1488	itr!=cornersToProcess.end();
1489	++itr)
1490	{
1491	TileCornerPair& tcp = *itr;
1492	if (tcp.first->_imagery.valid() && tcp.first->_imagery->_image.valid())
1493	{
1494	imagesToProcess.push_back(ImageCornerPair(tcp.first->_imagery->_image.get(),tcp.second));
1495	}
1496
1497	if (tcp.first->_terrain.valid() && tcp.first->_terrain->_heightField.valid())
1498	{
1499	heightFieldsToProcess.push_back(HeightFieldCornerPair(tcp.first->_terrain->_heightField.get(),tcp.second));
1500	}
1501	}
1502
1503	if (imagesToProcess.size()>1)
1504	{
1505	int red = 0;
1506	int green = 0;
1507	int blue = 0;
1508	// accumulate colours.
1509	ImageCornerList::iterator iitr;
1510	for(iitr=imagesToProcess.begin();
1511	iitr!=imagesToProcess.end();
1512	++iitr)
1513	{
1514	ImageCornerPair& icp = *iitr;
1515	unsigned char* data=0;
1516	switch(icp.second)
1517	{
1518	case LEFT_BELOW:
1519	data = icp.first->data(0,0);
1520	break;
1521	case BELOW_RIGHT:
1522	data = icp.first->data(icp.first->s()-1,0);
1523	break;
1524	case RIGHT_ABOVE:
1525	data = icp.first->data(icp.first->s()-1,icp.first->t()-1);
1526	break;
1527	case ABOVE_LEFT:
1528	data = icp.first->data(0,icp.first->t()-1);
1529	break;
1530	default :
1531	break;
1532	}
1533	red += *(data++);
1534	green += *(data++);
1535	blue += *(data);
1536	}
1537
1538	// divide them.
1539	red /= imagesToProcess.size();
1540	green /= imagesToProcess.size();
1541	blue /= imagesToProcess.size();
1542	// apply colour to corners.
1543	for(iitr=imagesToProcess.begin();
1544	iitr!=imagesToProcess.end();
1545	++iitr)
1546	{
1547	ImageCornerPair& icp = *iitr;
1548	unsigned char* data=0;
1549	switch(icp.second)
1550	{
1551	case LEFT_BELOW:
1552	data = icp.first->data(0,0);
1553	break;
1554	case BELOW_RIGHT:
1555	data = icp.first->data(icp.first->s()-1,0);
1556	break;
1557	case RIGHT_ABOVE:
1558	data = icp.first->data(icp.first->s()-1,icp.first->t()-1);
1559	break;
1560	case ABOVE_LEFT:
1561	data = icp.first->data(0,icp.first->t()-1);
1562	break;
1563	default :
1564	break;
1565	}
1566	*(data++) = red;
1567	*(data++) = green;
1568	*(data) = blue;
1569	}
1570
1571	}
1572
1573	if (heightFieldsToProcess.size()>1)
1574	{
1575	float height = 0;
1576	// accumulate heights.
1577	HeightFieldCornerList::iterator hitr;
1578	for(hitr=heightFieldsToProcess.begin();
1579	hitr!=heightFieldsToProcess.end();
1580	++hitr)
1581	{
1582	HeightFieldCornerPair& hfcp = *hitr;
1583	switch(hfcp.second)
1584	{
1585	case LEFT_BELOW:
1586	height += hfcp.first->getHeight(0,0);
1587	break;
1588	case BELOW_RIGHT:
1589	height += hfcp.first->getHeight(hfcp.first->getNumColumns()-1,0);
1590	break;
1591	case RIGHT_ABOVE:
1592	height += hfcp.first->getHeight(hfcp.first->getNumColumns()-1,hfcp.first->getNumRows()-1);
1593	break;
1594	case ABOVE_LEFT:
1595	height += hfcp.first->getHeight(0,hfcp.first->getNumRows()-1);
1596	break;
1597	default :
1598	break;
1599	}
1600	}
1601
1602	// divide them.
1603	height /= heightFieldsToProcess.size();
1604
1605	// apply height to corners.
1606	for(hitr=heightFieldsToProcess.begin();
1607	hitr!=heightFieldsToProcess.end();
1608	++hitr)
1609	{
1610	HeightFieldCornerPair& hfcp = *hitr;
1611	switch(hfcp.second)
1612	{
1613	case LEFT_BELOW:
1614	hfcp.first->setHeight(0,0,height);
1615	break;
1616	case BELOW_RIGHT:
1617	hfcp.first->setHeight(hfcp.first->getNumColumns()-1,0,height);
1618	break;
1619	case RIGHT_ABOVE:
1620	hfcp.first->setHeight(hfcp.first->getNumColumns()-1,hfcp.first->getNumRows()-1,height);
1621	break;
1622	case ABOVE_LEFT:
1623	hfcp.first->setHeight(0,hfcp.first->getNumRows()-1,height);
1624	break;
1625	default :
1626	break;
1627	}
1628	}
1629	}
1630
1631	}
1632
1633	const char* edgeString(DataSet::DestinationTile::Position position)
1634	{
1635	switch(position)
1636	{
1637	case DataSet::DestinationTile::LEFT: return "left";
1638	case DataSet::DestinationTile::BELOW: return "below";
1639	case DataSet::DestinationTile::RIGHT: return "right";
1640	case DataSet::DestinationTile::ABOVE: return "above";
1641	default : return "<not an edge>";
1642	}
1643	}
1644
1645	void DataSet::DestinationTile::setTileComplete(bool complete)
1646	{
1647	_complete = complete;
1648	osg::notify(osg::INFO)<<"setTileComplete("<<complete<<") for "<<_level<<"\t"<<_tileX<<"\t"<<_tileY<<std::endl;
1649	}
1650
1651	void DataSet::DestinationTile::equalizeEdge(Position position)
1652	{
1653	// don't need to equalize if already done.
1654	if (_equalized[position]) return;
1655
1656	DestinationTile* tile2 = _neighbour[position];
1657	Position position2 = (Position)((position+4)%NUMBER_OF_POSITIONS);
1658
1659	_equalized[position] = true;
1660
1661	// no neighbour of this edge so nothing to equalize.
1662	if (!tile2) return;
1663
1664	tile2->_equalized[position2]=true;
1665
1666	osg::Image* image1 = _imagery.valid()?_imagery->_image.get() : 0;
1667	osg::Image* image2 = tile2->_imagery.valid()?tile2->_imagery->_image.get() : 0;
1668
1669	//osg::notify(osg::INFO)<<"Equalizing edge "<<edgeString(position)<<" of \t"<<_level<<"\t"<<_tileX<<"\t"<<_tileY
1670	// <<" neighbour "<<tile2->_level<<"\t"<<tile2->_tileX<<"\t"<<tile2->_tileY<<std::endl;
1671
1672
1673	// if (_tileY==0) return;
1674
1675	if (image1 && image2 &&
1676	image1->getPixelFormat()==image2->getPixelFormat() &&
1677	image1->getDataType()==image2->getDataType() &&
1678	image1->getPixelFormat()==GL_RGB &&
1679	image1->getDataType()==GL_UNSIGNED_BYTE)
1680	{
1681
1682	//osg::notify(osg::INFO)<<" Equalizing image1= "<<image1<< " with image2 = "<<image2<<std::endl;
1683	//osg::notify(osg::INFO)<<" data1 = 0x"<<std::hex<<(int)image1->data()<<" with data2 = 0x"<<(int)image2->data()<<std::endl;
1684
1685	unsigned char* data1 = 0;
1686	unsigned char* data2 = 0;
1687	unsigned int delta1 = 0;
1688	unsigned int delta2 = 0;
1689	int num = 0;
1690
1691	switch(position)
1692	{
1693	case LEFT:
1694	data1 = image1->data(0,1); // LEFT hand side
1695	delta1 = image1->getRowSizeInBytes();
1696	data2 = image2->data(image2->s()-1,1); // RIGHT hand side
1697	delta2 = image2->getRowSizeInBytes();
1698	num = (image1->t()==image2->t())?image2->t()-2:0; // note miss out corners.
1699	//osg::notify(osg::INFO)<<" left "<<num<<std::endl;
1700	break;
1701	case BELOW:
1702	data1 = image1->data(1,0); // BELOW hand side
1703	delta1 = 3;
1704	data2 = image2->data(1,image2->t()-1); // ABOVE hand side
1705	delta2 = 3;
1706	num = (image1->s()==image2->s())?image2->s()-2:0; // note miss out corners.
1707	//osg::notify(osg::INFO)<<" below "<<num<<std::endl;
1708	break;
1709	case RIGHT:
1710	data1 = image1->data(image1->s()-1,1); // LEFT hand side
1711	delta1 = image1->getRowSizeInBytes();
1712	data2 = image2->data(0,1); // RIGHT hand side
1713	delta2 = image2->getRowSizeInBytes();
1714	num = (image1->t()==image2->t())?image2->t()-2:0; // note miss out corners.
1715	//osg::notify(osg::INFO)<<" right "<<num<<std::endl;
1716	break;
1717	case ABOVE:
1718	data1 = image1->data(1,image1->t()-1); // ABOVE hand side
1719	delta1 = 3;
1720	data2 = image2->data(1,0); // BELOW hand side
1721	delta2 = 3;
1722	num = (image1->s()==image2->s())?image2->s()-2:0; // note miss out corners.
1723	//osg::notify(osg::INFO)<<" above "<<num<<std::endl;
1724	break;
1725	default :
1726	//osg::notify(osg::INFO)<<" default "<<num<<std::endl;
1727	break;
1728	}
1729
1730	for(int i=0;i<num;++i)
1731	{
1732	unsigned char red = (unsigned char)((((int)data1+ (int)data2)/2));
1733	unsigned char green = (unsigned char)((((int)(data1+1))+ (int)((data2+1)))/2);
1734	unsigned char blue = (unsigned char)((((int)(data1+2))+ (int)((data2+2)))/2);
1735	#if 1
1736	*data1 = red;
1737	*(data1+1) = green;
1738	*(data1+2) = blue;
1739
1740	*data2 = red;
1741	*(data2+1) = green;
1742	*(data2+2) = blue;
1743	#endif
1744
1745	#if 0
1746	*data1 = 255;
1747	*(data1+1) = 0;
1748	*(data1+2) = 0;
1749
1750	*data2 = 0;
1751	*(data2+1) = 0;
1752	*(data2+2) = 0;
1753	#endif
1754	data1 += delta1;
1755	data2 += delta2;
1756
1757	//osg::notify(osg::INFO)<<" equalizing colour "<<(int)data1<<" "<<(int)data2<<std::endl;
1758
1759	}
1760
1761	}
1762
1763	osg::HeightField* heightField1 = _terrain.valid()?_terrain->_heightField.get():0;
1764	osg::HeightField* heightField2 = tile2->_terrain.valid()?tile2->_terrain->_heightField.get():0;
1765
1766	if (heightField1 && heightField2)
1767	{
1768	//osg::notify(osg::INFO)<<" Equalizing heightfield"<<std::endl;
1769
1770	float* data1 = 0;
1771	float* data2 = 0;
1772	unsigned int delta1 = 0;
1773	unsigned int delta2 = 0;
1774	int num = 0;
1775
1776	switch(position)
1777	{
1778	case LEFT:
1779	data1 = &(heightField1->getHeight(0,1)); // LEFT hand side
1780	delta1 = heightField1->getNumColumns();
1781	data2 = &(heightField2->getHeight(heightField2->getNumColumns()-1,1)); // RIGHT hand side
1782	delta2 = heightField2->getNumColumns();
1783	num = (heightField1->getNumRows()==heightField2->getNumRows())?heightField1->getNumRows()-2:0; // note miss out corners.
1784	break;
1785	case BELOW:
1786	data1 = &(heightField1->getHeight(1,0)); // BELOW hand side
1787	delta1 = 1;
1788	data2 = &(heightField2->getHeight(1,heightField2->getNumRows()-1)); // ABOVE hand side
1789	delta2 = 1;
1790	num = (heightField1->getNumColumns()==heightField2->getNumColumns())?heightField1->getNumColumns()-2:0; // note miss out corners.
1791	break;
1792	case RIGHT:
1793	data1 = &(heightField1->getHeight(heightField1->getNumColumns()-1,1)); // LEFT hand side
1794	delta1 = heightField1->getNumColumns();
1795	data2 = &(heightField2->getHeight(0,1)); // LEFT hand side
1796	delta2 = heightField2->getNumColumns();
1797	num = (heightField1->getNumRows()==heightField2->getNumRows())?heightField1->getNumRows()-2:0; // note miss out corners.
1798	break;
1799	case ABOVE:
1800	data1 = &(heightField1->getHeight(1,heightField1->getNumRows()-1)); // ABOVE hand side
1801	delta1 = 1;
1802	data2 = &(heightField2->getHeight(1,0)); // BELOW hand side
1803	delta2 = 1;
1804	num = (heightField1->getNumColumns()==heightField2->getNumColumns())?heightField1->getNumColumns()-2:0; // note miss out corners.
1805	break;
1806	default :
1807	break;
1808	}
1809
1810	for(int i=0;i<num;++i)
1811	{
1812	float z = (data1 + data2)/2.0f;
1813
1814	//osg::notify(osg::INFO)<<" equalizing height"<<std::endl;
1815
1816	*data1 = z;
1817	*data2 = z;
1818
1819	data1 += delta1;
1820	data2 += delta2;
1821	}
1822
1823	}
1824
1825	}
1826
1827	void DataSet::DestinationTile::equalizeBoundaries()
1828	{
1829	osg::notify(osg::INFO)<<"DataSet::DestinationTile::equalizeBoundaries()"<<std::endl;
1830
1831	equalizeCorner(LEFT_BELOW);
1832	equalizeCorner(BELOW_RIGHT);
1833	equalizeCorner(RIGHT_ABOVE);
1834	equalizeCorner(ABOVE_LEFT);
1835
1836	equalizeEdge(LEFT);
1837	equalizeEdge(BELOW);
1838	equalizeEdge(RIGHT);
1839	equalizeEdge(ABOVE);
1840	}
1841
1842
1843	void DataSet::DestinationTile::optimizeResolution()
1844	{
1845	if (_terrain.valid() && _terrain->_heightField.valid())
1846	{
1847	osg::HeightField* hf = _terrain->_heightField.get();
1848
1849	// compute min max of height field
1850	float minHeight = hf->getHeight(0,0);
1851	float maxHeight = minHeight;
1852	for(unsigned int r=0;r<hf->getNumRows();++r)
1853	{
1854	for(unsigned int c=0;c<hf->getNumColumns();++c)
1855	{
1856	float h = hf->getHeight(c,r);
1857	if (h<minHeight) minHeight = h;
1858	if (h>maxHeight) maxHeight = h;
1859	}
1860	}
1861
1862	if (minHeight==maxHeight)
1863	{
1864	osg::notify(osg::INFO)<<"*** We have a flat tile *** "<<std::endl;
1865
1866	unsigned int minimumSize = 8;
1867
1868	unsigned int numColumns = minimumSize;
1869	unsigned int numRows = minimumSize;
1870
1871	float ratio_y_over_x = (_extents.yMax()-_extents.yMin())/(_extents.xMax()-_extents.xMin());
1872	if (ratio_y_over_x > 1.2) numRows = (unsigned int)ceilf((float)numRows*ratio_y_over_x);
1873	else if (ratio_y_over_x < 0.8) numColumns = (unsigned int)ceilf((float)numColumns/ratio_y_over_x);
1874
1875
1876	hf->allocate(numColumns,numRows);
1877	hf->setOrigin(osg::Vec3(_extents.xMin(),_extents.yMin(),0.0f));
1878	hf->setXInterval((_extents.xMax()-_extents.xMin())/(float)(numColumns-1));
1879	hf->setYInterval((_extents.yMax()-_extents.yMin())/(float)(numRows-1));
1880
1881	for(unsigned int r=0;r<numRows;++r)
1882	{
1883	for(unsigned int c=0;c<numColumns;++c)
1884	{
1885	hf->setHeight(c,r,minHeight);
1886	}
1887	}
1888	}
1889	}
1890	}
1891
1892	osg::Node* DataSet::DestinationTile::createScene()
1893	{
1894	if (_dataSet->getGeometryType()==HEIGHT_FIELD)
1895	{
1896	return createHeightField();
1897	}
1898	else
1899	{
1900	return createPolygonal();
1901	}
1902	}
1903
1904	osg::StateSet* DataSet::DestinationTile::createStateSet()
1905	{
1906	if (!_imagery.valid() \|\| !_imagery->_image.valid()) return 0;
1907
1908	std::string imageName(_name+".rgb");
1909	_imagery->_image->setFileName(imageName.c_str());
1910
1911	osg::StateSet* stateset = new osg::StateSet;
1912	osg::Image* image = _imagery->_image.get();
1913
1914	osg::Texture2D* texture = new osg::Texture2D;
1915	texture->setImage(image);
1916	texture->setWrap(osg::Texture::WRAP_S,osg::Texture::CLAMP_TO_EDGE);
1917	texture->setWrap(osg::Texture::WRAP_T,osg::Texture::CLAMP_TO_EDGE);
1918	switch (_dataSet->getMipMappingMode())
1919	{
1920	case(DataSet::NO_MIP_MAPPING):
1921	{
1922	texture->setFilter(osg::Texture::MIN_FILTER,osg::Texture::LINEAR);
1923	texture->setFilter(osg::Texture::MAG_FILTER,osg::Texture::LINEAR);
1924	}
1925	break;
1926	case(DataSet::MIP_MAPPING_HARDWARE):
1927	{
1928	texture->setFilter(osg::Texture::MIN_FILTER,osg::Texture::LINEAR_MIPMAP_LINEAR);
1929	texture->setFilter(osg::Texture::MAG_FILTER,osg::Texture::LINEAR);
1930	}
1931	break;
1932	case(DataSet::MIP_MAPPING_IMAGERY):
1933	{
1934	texture->setFilter(osg::Texture::MIN_FILTER,osg::Texture::LINEAR_MIPMAP_LINEAR);
1935	texture->setFilter(osg::Texture::MAG_FILTER,osg::Texture::LINEAR);
1936	}
1937	break;
1938	}
1939
1940	texture->setMaxAnisotropy(_dataSet->getMaxAnisotropy());
1941	stateset->setTextureAttributeAndModes(0,texture,osg::StateAttribute::ON);
1942
1943	bool inlineImageFile = _dataSet->getDestinationTileExtension()==".ive";
1944	bool compressedImageSupported = inlineImageFile;
1945	bool mipmapImageSupported = inlineImageFile;
1946
1947	if (compressedImageSupported &&
1948	_dataSet->getTextureType()==COMPRESSED_TEXTURE &&
1949	(image->getPixelFormat()==GL_RGB \|\| image->getPixelFormat()==GL_RGBA))
1950	{
1951	texture->setInternalFormatMode(osg::Texture::USE_S3TC_DXT3_COMPRESSION);
1952
1953	osg::ref_ptr<osg::State> state = new osg::State;
1954
1955	// force the mip mapping off temporay if we intend the graphics hardware to do the mipmapping.
1956	if (_dataSet->getMipMappingMode()==DataSet::MIP_MAPPING_HARDWARE)
1957	texture->setFilter(osg::Texture::MIN_FILTER,osg::Texture::LINEAR);
1958
1959	// get OpenGL driver to create texture from image.
1960	texture->apply(*state);
1961
1962	image->readImageFromCurrentTexture(0,true);
1963
1964	// restore the mip mapping mode.
1965	if (_dataSet->getMipMappingMode()==DataSet::MIP_MAPPING_HARDWARE)
1966	texture->setFilter(osg::Texture::MIN_FILTER,osg::Texture::LINEAR_MIPMAP_LINEAR);
1967
1968	texture->setInternalFormatMode(osg::Texture::USE_IMAGE_DATA_FORMAT);
1969
1970	osg::notify(osg::INFO)<<">>>>>>>>>>>>>>>compressed image.<<<<<<<<<<<<<<"<<std::endl;
1971
1972	}
1973	else
1974	{
1975	if (_dataSet->getTextureType()==RGB_16_BIT && image->getPixelFormat()==GL_RGB)
1976	{
1977	image->scaleImage(image->s(),image->t(),image->r(),GL_UNSIGNED_SHORT_5_6_5);
1978	}
1979
1980	if (mipmapImageSupported && _dataSet->getMipMappingMode()==DataSet::MIP_MAPPING_IMAGERY)
1981	{
1982	osg::notify(osg::NOTICE)<<"Non compressress mipmapped not yet supported yet"<<std::endl;
1983	}
1984	}
1985
1986	if (!inlineImageFile)
1987	{
1988	osg::notify(osg::INFO)<<"Writing out imagery to "<<imageName<<std::endl;
1989	osgDB::writeImageFile(*_imagery->_image,_imagery->_image->getFileName().c_str());
1990	}
1991
1992	return stateset;
1993	}
1994
1995	osg::Node* DataSet::DestinationTile::createHeightField()
1996	{
1997	osg::ShapeDrawable* shapeDrawable = 0;
1998
1999	if (_terrain.valid() && _terrain->_heightField.valid())
2000	{
2001	osg::notify(osg::INFO)<<"--- Have terrain build tile ----"<<std::endl;
2002
2003	osg::HeightField* hf = _terrain->_heightField.get();
2004
2005	shapeDrawable = new osg::ShapeDrawable(hf);
2006
2007	hf->setSkirtHeight(shapeDrawable->getBound().radius()*0.01f);
2008	}
2009	else
2010	{
2011	osg::notify(osg::INFO)<<" No terrain to build tile from use flat terrain fallback "<<std::endl;
2012	// create a dummy height field to file in the gap
2013	osg::HeightField* hf = new osg::HeightField;
2014	hf->allocate(2,2);
2015	hf->setOrigin(osg::Vec3(_extents.xMin(),_extents.yMin(),0.0f));
2016	hf->setXInterval(_extents.xMax()-_extents.xMin());
2017	hf->setYInterval(_extents.yMax()-_extents.yMin());
2018
2019	shapeDrawable = new osg::ShapeDrawable(hf);
2020
2021	hf->setSkirtHeight(shapeDrawable->getBound().radius()*0.01f);
2022	}
2023
2024	if (!shapeDrawable) return 0;
2025
2026	osg::StateSet* stateset = createStateSet();
2027	if (stateset)
2028	{
2029	shapeDrawable->setStateSet(stateset);
2030	}
2031	else
2032	{
2033	shapeDrawable->setColor(_dataSet->getDefaultColor());
2034	}
2035
2036	osg::Geode* geode = new osg::Geode;
2037	geode->addDrawable(shapeDrawable);
2038
2039	return geode;
2040
2041	}
2042
2043
2044	static osg::Vec3 computeLocalPosition(const osg::Matrixd& worldToLocal, double X, double Y, double Z)
2045	{
2046	return osg::Vec3(XworldToLocal(0,0) + YworldToLocal(1,0) + Z*worldToLocal(2,0) + worldToLocal(3,0),
2047	XworldToLocal(0,1) + YworldToLocal(1,1) + Z*worldToLocal(2,1) + worldToLocal(3,1),
2048	XworldToLocal(0,2) + YworldToLocal(1,2) + Z*worldToLocal(2,2) + worldToLocal(3,2));
2049	}
2050
2051	static inline osg::Vec3 computeLocalSkirtVector(const osg::EllipsoidModel* et, const osg::HeightField* grid, unsigned int i, unsigned int j, float length, bool useLocalToTileTransform, const osg::Matrixd& localToWorld)
2052	{
2053	// no local to tile transform + mapping from lat+longs to XYZ so we need to use
2054	// a rotatated skirt vector - use the gravity vector.
2055	double longitude = grid->getOrigin().x()+grid->getXInterval()((double*)(i));
2056	double latitude = grid->getOrigin().y()+grid->getYInterval()((double*)(j));
2057	double midZ = grid->getOrigin().z();
2058	double X,Y,Z;
2059	et->convertLatLongHeightToXYZ(osg::DegreesToRadians(latitude),osg::DegreesToRadians(longitude),midZ,X,Y,Z);
2060	osg::Vec3 gravitationVector = et->computeLocalUpVector(X,Y,Z);
2061	gravitationVector.normalize();
2062
2063	if (useLocalToTileTransform) gravitationVector = osg::Matrixd::transform3x3(localToWorld,gravitationVector);
2064
2065	return gravitationVector * -length;
2066	}
2067
2068	osg::Node* DataSet::DestinationTile::createPolygonal()
2069	{
2070	osg::notify(osg::INFO)<<"--------- DataSet::DestinationTile::createDrawableGeometry() ------------- "<<std::endl;
2071
2072	const osg::EllipsoidModel* et = _dataSet->getEllipsoidModel();
2073	bool mapLatLongsToXYZ = _dataSet->getConvertFromGeographicToGeocentric() && et;
2074	bool useLocalToTileTransform = _dataSet->getUseLocalTileTransform();
2075
2076	osg::ref_ptr<osg::HeightField> grid = 0;
2077
2078	if (_terrain.valid() && _terrain->_heightField.valid())
2079	{
2080	osg::notify(osg::INFO)<<"--- Have terrain build tile ----"<<std::endl;
2081	grid = _terrain->_heightField.get();
2082	}
2083	else
2084	{
2085	unsigned int minimumSize = 8;
2086	unsigned int numColumns = minimumSize;
2087	unsigned int numRows = minimumSize;
2088
2089	if (mapLatLongsToXYZ)
2090	{
2091	float longitude_range = (_extents.xMax()-_extents.xMin());
2092	float latitude_range = (_extents.yMax()-_extents.yMin());
2093
2094	if (longitude_range>45.0) numColumns = (unsigned int)ceilf((float)numColumns*sqrtf(longitude_range/45.0));
2095	if (latitude_range>45.0) numRows = (unsigned int)ceilf((float)numRows*sqrtf(latitude_range/45.0));
2096
2097	osg::notify(osg::INFO)<<"numColumns = "<<numColumns<<" numRows="<<numRows<<std::endl;
2098	}
2099	else
2100	{
2101	float ratio_y_over_x = (_extents.yMax()-_extents.yMin())/(_extents.xMax()-_extents.xMin());
2102	if (ratio_y_over_x > 1.2) numRows = (unsigned int)ceilf((float)numRows*ratio_y_over_x);
2103	else if (ratio_y_over_x < 0.8) numColumns = (unsigned int)ceilf((float)numColumns/ratio_y_over_x);
2104	}
2105
2106	grid = new osg::HeightField;
2107	grid->allocate(numColumns,numRows);
2108	grid->setOrigin(osg::Vec3(_extents.xMin(),_extents.yMin(),0.0f));
2109	grid->setXInterval((_extents.xMax()-_extents.xMin())/(float)(numColumns-1));
2110	grid->setYInterval((_extents.yMax()-_extents.yMin())/(float)(numRows-1));
2111	}
2112
2113	if (!grid)
2114	{
2115	osg::notify(osg::INFO)<<" No terrain to build tile from use flat terrain fallback "<<std::endl;
2116
2117	return 0;
2118	}
2119
2120	bool createSkirt = true;
2121
2122	// compute sizes.
2123	unsigned int numColumns = grid->getNumColumns();
2124	unsigned int numRows = grid->getNumRows();
2125	unsigned int numVerticesInBody = numColumns*numRows;
2126	unsigned int numVerticesInSkirt = createSkirt ? numColumns2 + numRows2 - 4 : 0;
2127	unsigned int numVertices = numVerticesInBody+numVerticesInSkirt;
2128
2129
2130	// create the geometry.
2131	osg::Geometry* geometry = new osg::Geometry;
2132
2133	osg::Vec3Array& v = (new* osg::Vec3Array(numVertices));
2134	osg::Vec2Array& t = (new* osg::Vec2Array(numVertices));
2135	osg::UByte4Array& color = (new* osg::UByte4Array(1));
2136
2137	color[0].set(255,255,255,255);
2138
2139	osg::ref_ptr<osg::Vec3Array> n = new osg::Vec3Array(numVertices); // must use ref_ptr so the array isn't removed when smooothvisitor is used
2140
2141	float skirtRatio = _dataSet->getSkirtRatio();
2142	osg::Matrixd localToWorld;
2143	osg::Matrixd worldToLocal;
2144	osg::Vec3 skirtVector(0.0f,0.0f,0.0f);
2145
2146
2147	osg::Vec3 center_position(0.0f,0.0f,0.0f);
2148	osg::Vec3 center_normal(0.0f,0.0f,1.0f);
2149	osg::Vec3 transformed_center_normal(0.0f,0.0f,1.0f);
2150	double globe_radius = et ? et->getRadiusPolar() : 1.0;
2151	float skirtLength = _extents.radius()*skirtRatio;
2152
2153	bool useClusterCullingCallback = mapLatLongsToXYZ;
2154
2155	if (useLocalToTileTransform)
2156	{
2157	if (mapLatLongsToXYZ)
2158	{
2159	double midLong = grid->getOrigin().x()+grid->getXInterval()((double)(numColumns-1))0.5;
2160	double midLat = grid->getOrigin().y()+grid->getYInterval()((double)(numRows-1))0.5;
2161	double midZ = grid->getOrigin().z();
2162	et->computeLocalToWorldTransformFromLatLongHeight(osg::DegreesToRadians(midLat),osg::DegreesToRadians(midLong),midZ,localToWorld);
2163
2164	double minLong = grid->getOrigin().x();
2165	double minLat = grid->getOrigin().y();
2166
2167	double minX,minY,minZ;
2168	et->convertLatLongHeightToXYZ(osg::DegreesToRadians(minLat),osg::DegreesToRadians(minLong),midZ, minX,minY,minZ);
2169
2170	double midX,midY;
2171	et->convertLatLongHeightToXYZ(osg::DegreesToRadians(midLat),osg::DegreesToRadians(midLong),midZ, midX,midY,midZ);
2172
2173	double length = sqrt((midX-minX)(midX-minX) + (midY-minY)(midY-minY));
2174
2175	skirtLength = length*skirtRatio;
2176	skirtVector.set(0.0f,0.0f,-skirtLength);
2177
2178	center_normal.set(midX,midY,midZ);
2179	center_normal.normalize();
2180
2181	worldToLocal.invert(localToWorld);
2182
2183	center_position = computeLocalPosition(worldToLocal,midX,midY,midZ);
2184	transformed_center_normal = osg::Matrixd::transform3x3(localToWorld,center_normal);
2185
2186	}
2187	else
2188	{
2189	double midX = grid->getOrigin().x()+grid->getXInterval()((double)(numColumns-1))0.5;
2190	double midY = grid->getOrigin().y()+grid->getYInterval()((double)(numRows-1))0.5;
2191	double midZ = grid->getOrigin().z();
2192	localToWorld.makeTranslate(midX,midY,midZ);
2193	worldToLocal.invert(localToWorld);
2194
2195	skirtVector.set(0.0f,0.0f,-skirtLength);
2196	}
2197
2198	}
2199	else if (mapLatLongsToXYZ)
2200	{
2201	// no local to tile transform + mapping from lat+longs to XYZ so we need to use
2202	// a rotatated skirt vector - use the gravity vector.
2203	double midLong = grid->getOrigin().x()+grid->getXInterval()((double)(numColumns-1))0.5;
2204	double midLat = grid->getOrigin().y()+grid->getYInterval()((double)(numRows-1))0.5;
2205	double midZ = grid->getOrigin().z();
2206	double X,Y,Z;
2207	et->convertLatLongHeightToXYZ(osg::DegreesToRadians(midLat),osg::DegreesToRadians(midLong),midZ,X,Y,Z);
2208	osg::Vec3 gravitationVector = et->computeLocalUpVector(X,Y,Z);
2209	gravitationVector.normalize();
2210	skirtVector = gravitationVector * skirtLength;
2211	}
2212	else
2213	{
2214	skirtVector.set(0.0f,0.0f,-skirtLength);
2215	}
2216
2217
2218
2219	unsigned int vi=0;
2220	unsigned int r,c;
2221
2222	// populate the vertex/normal/texcoord arrays from the grid.
2223	double orig_X = grid->getOrigin().x();
2224	double delta_X = grid->getXInterval();
2225	double orig_Y = grid->getOrigin().y();
2226	double delta_Y = grid->getYInterval();
2227	double orig_Z = grid->getOrigin().z();
2228
2229
2230	float min_dot_product = 1.0f;
2231	float max_cluster_culling_height = 0.0f;
2232	float max_cluster_culling_radius = 0.0f;
2233
2234	for(r=0;r<numRows;++r)
2235	{
2236	for(c=0;c<numColumns;++c)
2237	{
2238	double X = orig_X + delta_X(double*)c;
2239	double Y = orig_Y + delta_Y(double*)r;
2240	double Z = orig_Z + grid->getHeight(c,r);
2241	double height = Z;
2242
2243	if (mapLatLongsToXYZ)
2244	{
2245	et->convertLatLongHeightToXYZ(osg::DegreesToRadians(Y),osg::DegreesToRadians(X),Z,
2246	X,Y,Z);
2247	}
2248
2249	if (useLocalToTileTransform)
2250	{
2251	v[vi] = computeLocalPosition(worldToLocal,X,Y,Z);
2252	}
2253	else
2254	{
2255	v[vi].set(X,Y,Z);
2256	}
2257
2258
2259	if (useClusterCullingCallback)
2260	{
2261	osg::Vec3 dv = v[vi] - center_position;
2262	double d = sqrt(dv.x()dv.x() + dv.y()dv.y() + dv.z()*dv.z());
2263	double theta = acos( globe_radius/ (globe_radius + fabs(height)) );
2264	double phi = 2.0 * asin (d0.5/globe_radius); // d/globe_radius;*
2265	double beta = theta+phi;
2266	double cutoff = osg::PI_2 - 0.1;
2267	//osg::notify(osg::INFO)<<"theta="<<theta<<"\tphi="<<phi<<" beta "<<beta<<std::endl;
2268	if (phi<cutoff && beta<cutoff)
2269	{
2270
2271	float local_dot_product = -sin(theta + phi);
2272	float local_m = globe_radius*( 1.0/ cos(theta+phi) - 1.0);
2273	min_dot_product = osg::minimum(min_dot_product, local_dot_product);
2274	max_cluster_culling_height = osg::maximum(max_cluster_culling_height,local_m);
2275	max_cluster_culling_radius = osg::maximum(max_cluster_culling_radius,static_cast<float>(beta*globe_radius));
2276	}
2277	else
2278	{
2279	//osg::notify(osg::INFO)<<"Turning off cluster culling for wrap around tile."<<std::endl;
2280	useClusterCullingCallback = false;
2281	}
2282	}
2283
2284	// note normal will need rotating.
2285	if (n.valid()) (*(n.get()))[vi] = grid->getNormal(c,r);
2286
2287	t[vi].x() = (c==numColumns-1)? 1.0f : (float)(c)/(float)(numColumns-1);
2288	t[vi].y() = (r==numRows-1)? 1.0f : (float)(r)/(float)(numRows-1);
2289
2290	++vi;
2291
2292	}
2293	}
2294
2295
2296
2297	//geometry->setUseDisplayList(false);
2298	geometry->setVertexArray(&v);
2299
2300	if (n.valid())
2301	{
2302	geometry->setNormalArray(n.get());
2303	geometry->setNormalBinding(osg::Geometry::BIND_PER_VERTEX);
2304	}
2305
2306	geometry->setColorArray(&color);
2307	geometry->setColorBinding(osg::Geometry::BIND_OVERALL);
2308	geometry->setTexCoordArray(0,&t);
2309
2310	osg::DrawElementsUShort& drawElements = (new* osg::DrawElementsUShort(GL_TRIANGLES,23(numColumns-1)*(numRows-1)));
2311	geometry->addPrimitiveSet(&drawElements);
2312	int ei=0;
2313	for(r=0;r<numRows-1;++r)
2314	{
2315	for(c=0;c<numColumns-1;++c)
2316	{
2317	unsigned short i00 = (r)*numColumns+c;
2318	unsigned short i10 = (r)*numColumns+c+1;
2319	unsigned short i01 = (r+1)*numColumns+c;
2320	unsigned short i11 = (r+1)*numColumns+c+1;
2321
2322	float diff_00_11 = fabsf(v[i00].z()-v[i11].z());
2323	float diff_01_10 = fabsf(v[i01].z()-v[i10].z());
2324	if (diff_00_11<diff_01_10)
2325	{
2326	// diagonal between 00 and 11
2327	drawElements[ei++] = i00;
2328	drawElements[ei++] = i10;
2329	drawElements[ei++] = i11;
2330
2331	drawElements[ei++] = i00;
2332	drawElements[ei++] = i11;
2333	drawElements[ei++] = i01;
2334	}
2335	else
2336	{
2337	// diagonal between 01 and 10
2338	drawElements[ei++] = i01;
2339	drawElements[ei++] = i00;
2340	drawElements[ei++] = i10;
2341
2342	drawElements[ei++] = i01;
2343	drawElements[ei++] = i10;
2344	drawElements[ei++] = i11;
2345	}
2346	}
2347	}
2348
2349	#if 1
2350	osgUtil::SmoothingVisitor sv;
2351	sv.smooth(geometry); // this will replace the normal vector with a new one*
2352
2353	// now we have to reassign the normals back to the orignal pointer.
2354	n = geometry->getNormalArray();
2355	if (n.valid() && n->size()!=numVertices) n->resize(numVertices);
2356	#endif
2357
2358
2359	if (useClusterCullingCallback)
2360	{
2361	// set up cluster cullling,
2362	osg::ClusterCullingCallback* ccc = new osg::ClusterCullingCallback;
2363
2364	ccc->set(center_position + transformed_center_normal*max_cluster_culling_height ,
2365	transformed_center_normal,
2366	min_dot_product,
2367	max_cluster_culling_radius);
2368	geometry->setCullCallback(ccc);
2369	}
2370
2371	osgUtil::Simplifier::IndexList pointsToProtectDuringSimplification;
2372
2373	if (numVerticesInSkirt>0)
2374	{
2375	osg::DrawElementsUShort& skirtDrawElements = (new* osg::DrawElementsUShort(GL_QUAD_STRIP,2*numVerticesInSkirt+2));
2376	geometry->addPrimitiveSet(&skirtDrawElements);
2377	int ei=0;
2378	int firstSkirtVertexIndex = vi;
2379	// create bottom skirt vertices
2380	r=0;
2381	for(c=0;c<numColumns-1;++c)
2382	{
2383	// assign indices to primitive set
2384	skirtDrawElements[ei++] = (r)*numColumns+c;
2385	skirtDrawElements[ei++] = vi;
2386
2387	// mark these points as protected to prevent them from being removed during simplification
2388	pointsToProtectDuringSimplification.push_back((r)*numColumns+c);
2389	pointsToProtectDuringSimplification.push_back(vi);
2390
2391	osg::Vec3 localSkirtVector = !mapLatLongsToXYZ ?
2392	skirtVector :
2393	computeLocalSkirtVector(et, grid.get(), c, r, skirtLength, useLocalToTileTransform, localToWorld);
2394
2395	// add in the new point on the bottom of the skirt
2396	v[vi] = v[(r)*numColumns+c]+localSkirtVector;
2397	if (n.valid()) (n)[vi] = (n)[r*numColumns+c];
2398	t[vi++] = t[(r)*numColumns+c];
2399	}
2400	// create right skirt vertices
2401	c=numColumns-1;
2402	for(r=0;r<numRows-1;++r)
2403	{
2404	// assign indices to primitive set
2405	skirtDrawElements[ei++] = (r)*numColumns+c;
2406	skirtDrawElements[ei++] = vi;
2407
2408	// mark these points as protected to prevent them from being removed during simplification
2409	pointsToProtectDuringSimplification.push_back((r)*numColumns+c);
2410	pointsToProtectDuringSimplification.push_back(vi);
2411
2412	osg::Vec3 localSkirtVector = !mapLatLongsToXYZ ?
2413	skirtVector :
2414	computeLocalSkirtVector(et, grid.get(), c, r, skirtLength, useLocalToTileTransform, localToWorld);
2415
2416	// add in the new point on the bottom of the skirt
2417	v[vi] = v[(r)*numColumns+c]+localSkirtVector;
2418	if (n.valid()) (n)[vi] = (n)[(r)*numColumns+c];
2419	t[vi++] = t[(r)*numColumns+c];
2420	}
2421	// create top skirt vertices
2422	r=numRows-1;
2423	for(c=numColumns-1;c>0;--c)
2424	{
2425	// assign indices to primitive set
2426	skirtDrawElements[ei++] = (r)*numColumns+c;
2427	skirtDrawElements[ei++] = vi;
2428
2429	// mark these points as protected to prevent them from being removed during simplification
2430	pointsToProtectDuringSimplification.push_back((r)*numColumns+c);
2431	pointsToProtectDuringSimplification.push_back(vi);
2432
2433	osg::Vec3 localSkirtVector = !mapLatLongsToXYZ ?
2434	skirtVector :
2435	computeLocalSkirtVector(et, grid.get(), c, r, skirtLength, useLocalToTileTransform, localToWorld);
2436
2437	// add in the new point on the bottom of the skirt
2438	v[vi] = v[(r)*numColumns+c]+localSkirtVector;
2439	if (n.valid()) (n)[vi] = (n)[(r)*numColumns+c];
2440	t[vi++] = t[(r)*numColumns+c];
2441	}
2442	// create left skirt vertices
2443	c=0;
2444	for(r=numRows-1;r>0;--r)
2445	{
2446	// assign indices to primitive set
2447	skirtDrawElements[ei++] = (r)*numColumns+c;
2448	skirtDrawElements[ei++] = vi;
2449
2450	// mark these points as protected to prevent them from being removed during simplification
2451	pointsToProtectDuringSimplification.push_back((r)*numColumns+c);
2452	pointsToProtectDuringSimplification.push_back(vi);
2453
2454	osg::Vec3 localSkirtVector = !mapLatLongsToXYZ ?
2455	skirtVector :
2456	computeLocalSkirtVector(et, grid.get(), c, r, skirtLength, useLocalToTileTransform, localToWorld);
2457
2458	// add in the new point on the bottom of the skirt
2459	v[vi] = v[(r)*numColumns+c]+localSkirtVector;
2460	if (n.valid()) (n)[vi] = (n)[(r)*numColumns+c];
2461	t[vi++] = t[(r)*numColumns+c];
2462	}
2463	skirtDrawElements[ei++] = 0;
2464	skirtDrawElements[ei++] = firstSkirtVertexIndex;
2465	}
2466
2467	if (n.valid())
2468	{
2469	geometry->setNormalArray(n.get());
2470	geometry->setNormalBinding(osg::Geometry::BIND_PER_VERTEX);
2471	}
2472
2473
2474	osg::StateSet* stateset = createStateSet();
2475	if (stateset)
2476	{
2477	geometry->setStateSet(stateset);
2478	}
2479	else
2480	{
2481	osg::Vec4Array* colours = new osg::Vec4Array(1);
2482	(*colours)[0] = _dataSet->getDefaultColor();
2483
2484	geometry->setColorArray(colours);
2485	geometry->setColorBinding(osg::Geometry::BIND_OVERALL);
2486	}
2487
2488
2489	osg::Geode* geode = new osg::Geode;
2490	geode->addDrawable(geometry);
2491
2492	#if 1
2493	osgDB::writeNodeFile(geode,"NodeBeforeSimplification.osg"*);
2494	#endif
2495
2496	unsigned int targetMaxNumVertices = 2048;
2497	float sample_ratio = (numVertices <= targetMaxNumVertices) ? 1.0f : (float)targetMaxNumVertices/(float)numVertices;
2498
2499	osgUtil::Simplifier simplifier(sample_ratio,geometry->getBound().radius()/2000.0f);
2500
2501	simplifier.simplify(geometry, pointsToProtectDuringSimplification); // this will replace the normal vector with a new one*
2502
2503
2504	osgUtil::TriStripVisitor tsv;
2505	tsv.setMinStripSize(3);
2506	tsv.stripify(*geometry);
2507
2508
2509
2510	if (useLocalToTileTransform)
2511	{
2512	osg::MatrixTransform* mt = new osg::MatrixTransform;
2513	mt->setMatrix(localToWorld);
2514	mt->addChild(geode);
2515
2516	bool addLocalAxes = false;
2517	if (addLocalAxes)
2518	{
2519	float s = geode->getBound().radius()*0.5f;
2520	osg::MatrixTransform* scaleAxis = new osg::MatrixTransform;
2521	scaleAxis->setMatrix(osg::Matrix::scale(s,s,s));
2522	scaleAxis->addChild(osgDB::readNodeFile("axes.osg"));
2523	mt->addChild(scaleAxis);
2524	}
2525
2526	return mt;
2527	}
2528	else
2529	{
2530	return geode;
2531	}
2532	}
2533
2534	void DataSet::DestinationTile::readFrom(CompositeSource* sourceGraph)
2535	{
2536	allocate();
2537
2538	osg::notify(osg::INFO)<<"DestinationTile::readFrom() "<<std::endl;
2539	for(CompositeSource::source_iterator itr(sourceGraph);itr.valid();++itr)
2540	{
2541
2542	Source* source = itr->get();
2543	if (source && _level>=source->getMinLevel() && _level<=source->getMaxLevel())
2544	{
2545	SourceData* data = (*itr)->getSourceData();
2546	if (data)
2547	{
2548	osg::notify(osg::INFO)<<"DataSet::DestinationTile::readFrom -> SourceData::read() "<<std::endl;
2549	osg::notify(osg::INFO)<<" destination._level="<<_level<<"\t"<<source->getMinLevel()<<"\t"<<source->getMaxLevel()<<std::endl;
2550	if (_imagery.valid()) data->read(*_imagery);
2551	if (_terrain.valid()) data->read(*_terrain);
2552	}
2553	}
2554	}
2555
2556	optimizeResolution();
2557
2558	}
2559
2560	void DataSet::DestinationTile::unrefData()
2561	{
2562	_imagery = 0;
2563	_terrain = 0;
2564	_models = 0;
2565	}
2566
2567	void DataSet::DestinationTile::addRequiredResolutions(CompositeSource* sourceGraph)
2568	{
2569	for(CompositeSource::source_iterator itr(sourceGraph);itr.valid();++itr)
2570	{
2571	Source* source = itr->get();
2572	if (source && source->intersects(*this))
2573	{
2574	if (source->getType()==Source::IMAGE)
2575	{
2576	unsigned int numCols,numRows;
2577	double resX, resY;
2578	if (computeImageResolution(numCols,numRows,resX,resY))
2579	{
2580	source->addRequiredResolution(resX,resY);
2581	}
2582	}
2583
2584	if (source->getType()==Source::HEIGHT_FIELD)
2585	{
2586	unsigned int numCols,numRows;
2587	double resX, resY;
2588	if (computeTerrainResolution(numCols,numRows,resX,resY))
2589	{
2590	source->addRequiredResolution(resX,resY);
2591	}
2592	}
2593	}
2594
2595	}
2596	}
2597
2598	void DataSet::CompositeDestination::computeNeighboursFromQuadMap()
2599	{
2600	// handle leaves
2601	for(TileList::iterator titr=_tiles.begin();
2602	titr!=_tiles.end();
2603	++titr)
2604	{
2605	(*titr)->computeNeighboursFromQuadMap();
2606	}
2607
2608	// handle chilren
2609	for(ChildList::iterator citr=_children.begin();
2610	citr!=_children.end();
2611	++citr)
2612	{
2613	(*citr)->computeNeighboursFromQuadMap();
2614	}
2615	}
2616
2617	void DataSet::CompositeDestination::addRequiredResolutions(CompositeSource* sourceGraph)
2618	{
2619	// handle leaves
2620	for(TileList::iterator titr=_tiles.begin();
2621	titr!=_tiles.end();
2622	++titr)
2623	{
2624	(*titr)->addRequiredResolutions(sourceGraph);
2625	}
2626
2627	// handle chilren
2628	for(ChildList::iterator citr=_children.begin();
2629	citr!=_children.end();
2630	++citr)
2631	{
2632	(*citr)->addRequiredResolutions(sourceGraph);
2633	}
2634	}
2635
2636	void DataSet::CompositeDestination::readFrom(CompositeSource* sourceGraph)
2637	{
2638	osg::notify(osg::INFO)<<"CompositeDestination::readFrom() "<<std::endl;
2639
2640	// handle leaves
2641	for(TileList::iterator titr=_tiles.begin();
2642	titr!=_tiles.end();
2643	++titr)
2644	{
2645	(*titr)->readFrom(sourceGraph);
2646	}
2647
2648	// handle chilren
2649	for(ChildList::iterator citr=_children.begin();
2650	citr!=_children.end();
2651	++citr)
2652	{
2653	(*citr)->readFrom(sourceGraph);
2654	}
2655	}
2656
2657	void DataSet::CompositeDestination::equalizeBoundaries()
2658	{
2659	// handle leaves
2660	for(TileList::iterator titr=_tiles.begin();
2661	titr!=_tiles.end();
2662	++titr)
2663	{
2664	(*titr)->equalizeBoundaries();
2665	}
2666
2667	// handle chilren
2668	for(ChildList::iterator citr=_children.begin();
2669	citr!=_children.end();
2670	++citr)
2671	{
2672	(*citr)->equalizeBoundaries();
2673	}
2674
2675	}
2676
2677
2678	osg::Node* DataSet::CompositeDestination::createScene()
2679	{
2680	if (_children.empty() && _tiles.empty()) return 0;
2681
2682	if (_children.empty() && _tiles.size()==1) return _tiles.front()->createScene();
2683
2684	if (_tiles.empty() && _children.size()==1) return _children.front()->createScene();
2685
2686	if (_type==GROUP)
2687	{
2688	osg::Group* group = new osg::Group;
2689	for(TileList::iterator titr=_tiles.begin();
2690	titr!=_tiles.end();
2691	++titr)
2692	{
2693	osg::Node* node = (*titr)->createScene();
2694	if (node) group->addChild(node);
2695	}
2696
2697	// handle chilren
2698	for(ChildList::iterator citr=_children.begin();
2699	citr!=_children.end();
2700	++citr)
2701	{
2702	osg::Node* node = (*citr)->createScene();
2703	if (node) group->addChild(node);
2704	}
2705	return group;
2706	}
2707
2708	// must be either a LOD or a PagedLOD
2709
2710	typedef std::vector<osg::Node*> NodeList;
2711	typedef std::map<float,NodeList> RangeNodeListMap;
2712	RangeNodeListMap rangeNodeListMap;
2713
2714	// insert local tiles into range map
2715	for(TileList::iterator titr=_tiles.begin();
2716	titr!=_tiles.end();
2717	++titr)
2718	{
2719	osg::Node* node = (*titr)->createScene();
2720	if (node) rangeNodeListMap[_maxVisibleDistance].push_back(node);
2721	}
2722
2723	// handle chilren
2724	for(ChildList::iterator citr=_children.begin();
2725	citr!=_children.end();
2726	++citr)
2727	{
2728	osg::Node* node = (*citr)->createScene();
2729	if (node) rangeNodeListMap[(*citr)->_maxVisibleDistance].push_back(node);
2730	}
2731
2732	osg::LOD* lod = new osg::LOD;
2733
2734	unsigned int childNum = 0;
2735	for(RangeNodeListMap::reverse_iterator rnitr=rangeNodeListMap.rbegin();
2736	rnitr!=rangeNodeListMap.rend();
2737	++rnitr)
2738	{
2739	float maxVisibleDistance = rnitr->first;
2740	NodeList& nodeList = rnitr->second;
2741
2742	if (childNum==0)
2743	{
2744	// by deafult make the first child have a very visible distance so its always seen
2745	maxVisibleDistance = 1e10;
2746	}
2747	else
2748	{
2749	// set the previous child's minimum visible distance range
2750	lod->setRange(childNum-1,maxVisibleDistance,lod->getMaxRange(childNum-1));
2751	}
2752
2753	osg::Node* child = 0;
2754	if (nodeList.size()==1)
2755	{
2756	child = nodeList.front();
2757	}
2758	else if (nodeList.size()>1)
2759	{
2760	osg::Group* group = new osg::Group;
2761	for(NodeList::iterator itr=nodeList.begin();
2762	itr!=nodeList.end();
2763	++itr)
2764	{
2765	group->addChild(*itr);
2766	}
2767	child = group;
2768	}
2769
2770	if (child)
2771	{
2772
2773	lod->addChild(child,0,maxVisibleDistance);
2774
2775	++childNum;
2776	}
2777	}
2778	return lod;
2779	}
2780
2781	bool DataSet::CompositeDestination::areSubTilesComplete()
2782	{
2783	for(ChildList::iterator citr=_children.begin();
2784	citr!=_children.end();
2785	++citr)
2786	{
2787	for(TileList::iterator itr=(*citr)->_tiles.begin();
2788	itr!=(*citr)->_tiles.end();
2789	++itr)
2790	{
2791	if (!(*itr)->getTileComplete())
2792	{
2793	return false;
2794	}
2795	}
2796	}
2797	return true;
2798	}
2799
2800	std::string DataSet::CompositeDestination::getSubTileName()
2801	{
2802	return _name+"_subtile"+_dataSet->getDestinationTileExtension();
2803	}
2804
2805	osg::Node* DataSet::CompositeDestination::createPagedLODScene()
2806	{
2807	if (_children.empty() && _tiles.empty()) return 0;
2808
2809	if (_children.empty() && _tiles.size()==1) return _tiles.front()->createScene();
2810
2811	if (_tiles.empty() && _children.size()==1) return _children.front()->createPagedLODScene();
2812
2813	if (_type==GROUP)
2814	{
2815	osg::Group* group = new osg::Group;
2816	for(TileList::iterator titr=_tiles.begin();
2817	titr!=_tiles.end();
2818	++titr)
2819	{
2820	osg::Node* node = (*titr)->createScene();
2821	if (node) group->addChild(node);
2822	}
2823
2824	// handle chilren
2825	for(ChildList::iterator citr=_children.begin();
2826	citr!=_children.end();
2827	++citr)
2828	{
2829	osg::Node* node = (*citr)->createScene();
2830	if (node) group->addChild(node);
2831	}
2832	return group;
2833	}
2834
2835	// must be either a LOD or a PagedLOD
2836
2837	typedef std::vector<osg::Node*> NodeList;
2838
2839	// collect all the local tiles
2840	NodeList tileNodes;
2841	for(TileList::iterator titr=_tiles.begin();
2842	titr!=_tiles.end();
2843	++titr)
2844	{
2845	osg::Node* node = (*titr)->createScene();
2846	if (node) tileNodes.push_back(node);
2847	}
2848
2849	float cutOffDistance = -FLT_MAX;
2850	for(ChildList::iterator citr=_children.begin();
2851	citr!=_children.end();
2852	++citr)
2853	{
2854	cutOffDistance = osg::maximum(cutOffDistance,(*citr)->_maxVisibleDistance);
2855	}
2856
2857	osg::PagedLOD* pagedLOD = new osg::PagedLOD;
2858
2859	float farDistance = _dataSet->getMaximumVisibleDistanceOfTopLevel();
2860	if (tileNodes.size()==1)
2861	{
2862	pagedLOD->addChild(tileNodes.front());
2863	}
2864	else if (tileNodes.size()>1)
2865	{
2866	osg::Group* group = new osg::Group;
2867	for(NodeList::iterator itr=tileNodes.begin();
2868	itr != tileNodes.end();
2869	++itr)
2870	{
2871	group->addChild(*itr);
2872	}
2873	pagedLOD->addChild(group);
2874	}
2875
2876	cutOffDistance = pagedLOD->getBound().radius()*_dataSet->getRadiusToMaxVisibleDistanceRatio();
2877
2878	pagedLOD->setRange(0,cutOffDistance,farDistance);
2879
2880	pagedLOD->setFileName(1,getSubTileName());
2881	pagedLOD->setRange(1,0,cutOffDistance);
2882
2883	if (pagedLOD->getNumChildren()>0)
2884	pagedLOD->setCenter(pagedLOD->getBound().center());
2885
2886	return pagedLOD;
2887	}
2888
2889	osg::Node* DataSet::CompositeDestination::createSubTileScene()
2890	{
2891	if (_type==GROUP \|\|
2892	_children.empty() \|\|
2893	_tiles.empty()) return 0;
2894
2895	// handle chilren
2896	typedef std::vector<osg::Node*> NodeList;
2897	NodeList nodeList;
2898	for(ChildList::iterator citr=_children.begin();
2899	citr!=_children.end();
2900	++citr)
2901	{
2902	osg::Node* node = (*citr)->createPagedLODScene();
2903	if (node) nodeList.push_back(node);
2904	}
2905
2906	if (nodeList.size()==1)
2907	{
2908	return nodeList.front();
2909	}
2910	else if (nodeList.size()>1)
2911	{
2912	osg::Group* group = new osg::Group;
2913	for(NodeList::iterator itr=nodeList.begin();
2914	itr!=nodeList.end();
2915	++itr)
2916	{
2917	group->addChild(*itr);
2918	}
2919	return group;
2920	}
2921	else
2922	{
2923	return 0;
2924	}
2925	}
2926
2927	void DataSet::CompositeDestination::unrefSubTileData()
2928	{
2929	for(CompositeDestination::ChildList::iterator citr=_children.begin();
2930	citr!=_children.end();
2931	++citr)
2932	{
2933	(*citr)->unrefLocalData();
2934	}
2935	}
2936
2937	void DataSet::CompositeDestination::unrefLocalData()
2938	{
2939	for(CompositeDestination::TileList::iterator titr=_tiles.begin();
2940	titr!=_tiles.end();
2941	++titr)
2942	{
2943	DestinationTile* tile = titr->get();
2944	osg::notify(osg::INFO)<<" unref tile level="<<tile->_level<<" X="<<tile->_tileX<<" Y="<<tile->_tileY<<std::endl;
2945	tile->unrefData();
2946	}
2947	}
2948
2949	///////////////////////////////////////////////////////////////////////////////
2950
2951	DataSet::DataSet()
2952	{
2953	init();
2954
2955	_maximumTileImageSize = 256;
2956	_maximumTileTerrainSize = 64;
2957
2958	_maximumVisiableDistanceOfTopLevel = 1e10;
2959
2960	_radiusToMaxVisibleDistanceRatio = 7.0f;
2961	_verticalScale = 1.0f;
2962	_skirtRatio = 0.02f;
2963
2964	_convertFromGeographicToGeocentric = false;
2965
2966	_defaultColor.set(0.5f,0.5f,1.0f,1.0f);
2967	_databaseType = PagedLOD_DATABASE;
2968	_geometryType = POLYGONAL;
2969	_textureType = COMPRESSED_TEXTURE;
2970	_maxAnisotropy = 1.0;
2971	_mipMappingMode = MIP_MAPPING_IMAGERY;
2972
2973	_useLocalTileTransform = true;
2974
2975	_decorateWithCoordinateSystemNode = true;
2976
2977	setEllipsoidModel(new osg::EllipsoidModel());
2978	}
2979
2980	void DataSet::init()
2981	{
2982	static bool s_initialized = false;
2983	if (!s_initialized)
2984	{
2985	s_initialized = true;
2986	GDALAllRegister();
2987	}
2988	}
2989
2990	void DataSet::addSource(Source* source)
2991	{
2992	if (!source) return;
2993
2994	if (!_sourceGraph.valid()) _sourceGraph = new CompositeSource;
2995
2996	_sourceGraph->_sourceList.push_back(source);
2997	}
2998
2999	void DataSet::addSource(CompositeSource* composite)
3000	{
3001	if (!composite) return;
3002
3003	if (!_sourceGraph.valid()) _sourceGraph = new CompositeSource;
3004
3005	_sourceGraph->_children.push_back(composite);
3006	}
3007
3008	void DataSet::loadSources()
3009	{
3010	for(CompositeSource::source_iterator itr(_sourceGraph.get());itr.valid();++itr)
3011	{
3012	(*itr)->loadSourceData();
3013	}
3014	}
3015
3016	DataSet::CompositeDestination* DataSet::createDestinationGraph(CompositeDestination* parent,
3017	osg::CoordinateSystemNode* cs,
3018	const osg::BoundingBox& extents,
3019	unsigned int maxImageSize,
3020	unsigned int maxTerrainSize,
3021	unsigned int currentLevel,
3022	unsigned int currentX,
3023	unsigned int currentY,
3024	unsigned int maxNumLevels)
3025	{
3026
3027	DataSet::CompositeDestination* destinationGraph = new DataSet::CompositeDestination(cs,extents);
3028
3029	destinationGraph->_maxVisibleDistance = extents.radius()*getRadiusToMaxVisibleDistanceRatio();
3030
3031	// first create the topmost tile
3032
3033	// create the name
3034	std::ostringstream os;
3035	os << _tileBasename << "_L"<<currentLevel<<"_X"<<currentX<<"_Y"<<currentY;
3036
3037	destinationGraph->_parent = parent;
3038	destinationGraph->_name = os.str();
3039	destinationGraph->_level = currentLevel;
3040	destinationGraph->_tileX = currentX;
3041	destinationGraph->_tileY = currentY;
3042	destinationGraph->_dataSet = this;
3043
3044
3045	DestinationTile* tile = new DestinationTile;
3046	tile->_name = destinationGraph->_name;
3047	tile->_level = currentLevel;
3048	tile->_tileX = currentX;
3049	tile->_tileY = currentY;
3050	tile->_dataSet = this;
3051	tile->_cs = cs;
3052	tile->_extents = extents;
3053	tile->setMaximumImagerySize(maxImageSize,maxImageSize);
3054	tile->setMaximumTerrainSize(maxTerrainSize,maxTerrainSize);
3055	tile->computeMaximumSourceResolution(_sourceGraph.get());
3056
3057	insertTileToQuadMap(destinationGraph);
3058
3059	if (currentLevel>=maxNumLevels-1 \|\| currentLevel>=tile->_maxSourceLevel)
3060	{
3061	// bottom level can't divide any further.
3062	destinationGraph->_tiles.push_back(tile);
3063	}
3064	else
3065	{
3066	destinationGraph->_type = LOD;
3067	destinationGraph->_tiles.push_back(tile);
3068
3069	bool needToDivideX = false;
3070	bool needToDivideY = false;
3071
3072	unsigned int texture_numColumns;
3073	unsigned int texture_numRows;
3074	double texture_dx;
3075	double texture_dy;
3076	if (tile->computeImageResolution(texture_numColumns,texture_numRows,texture_dx,texture_dy))
3077	{
3078	if (texture_dx>tile->_imagery_maxSourceResolutionX) needToDivideX = true;
3079	if (texture_dy>tile->_imagery_maxSourceResolutionY) needToDivideY = true;
3080	}
3081
3082	unsigned int dem_numColumns;
3083	unsigned int dem_numRows;
3084	double dem_dx;
3085	double dem_dy;
3086	if (tile->computeTerrainResolution(dem_numColumns,dem_numRows,dem_dx,dem_dy))
3087	{
3088	if (dem_dx>tile->_terrain_maxSourceResolutionX) needToDivideX = true;
3089	if (dem_dy>tile->_terrain_maxSourceResolutionY) needToDivideY = true;
3090	}
3091
3092	float xCenter = (extents.xMin()+extents.xMax())*0.5f;
3093	float yCenter = (extents.yMin()+extents.yMax())*0.5f;
3094
3095	unsigned int newLevel = currentLevel+1;
3096	unsigned int newX = currentX*2;
3097	unsigned int newY = currentY*2;
3098
3099	if (needToDivideX && needToDivideY)
3100	{
3101	float aspectRatio = (extents.yMax()- extents.yMin())/(extents.xMax()- extents.xMin());
3102
3103	if (aspectRatio>1.414) needToDivideX = false;
3104	else if (aspectRatio<.707) needToDivideY = false;
3105	}
3106
3107	if (needToDivideX && needToDivideY)
3108	{
3109	osg::notify(osg::INFO)<<"Need to Divide X + Y for level "<<currentLevel<<std::endl;
3110	// create four tiles.
3111	osg::BoundingBox bottom_left(extents.xMin(),extents.yMin(),extents.zMin(),xCenter,yCenter,extents.zMax());
3112	osg::BoundingBox bottom_right(xCenter,extents.yMin(),extents.zMin(),extents.xMax(),yCenter,extents.zMax());
3113	osg::BoundingBox top_left(extents.xMin(),yCenter,extents.zMin(),xCenter,extents.yMax(),extents.zMax());
3114	osg::BoundingBox top_right(xCenter,yCenter,extents.zMin(),extents.xMax(),extents.yMax(),extents.zMax());
3115
3116	destinationGraph->_children.push_back(createDestinationGraph(destinationGraph,
3117	cs,
3118	bottom_left,
3119	maxImageSize,
3120	maxTerrainSize,
3121	newLevel,
3122	newX,
3123	newY,
3124	maxNumLevels));
3125
3126	destinationGraph->_children.push_back(createDestinationGraph(destinationGraph,
3127	cs,
3128	bottom_right,
3129	maxImageSize,
3130	maxTerrainSize,
3131	newLevel,
3132	newX+1,
3133	newY,
3134	maxNumLevels));
3135
3136	destinationGraph->_children.push_back(createDestinationGraph(destinationGraph,
3137	cs,
3138	top_left,
3139	maxImageSize,
3140	maxTerrainSize,
3141	newLevel,
3142	newX,
3143	newY+1,
3144	maxNumLevels));
3145
3146	destinationGraph->_children.push_back(createDestinationGraph(destinationGraph,
3147	cs,
3148	top_right,
3149	maxImageSize,
3150	maxTerrainSize,
3151	newLevel,
3152	newX+1,
3153	newY+1,
3154	maxNumLevels));
3155
3156	// Set all there max distance to the same value to ensure the same LOD bining.
3157	float cutOffDistance = destinationGraph->_maxVisibleDistance*0.5f;
3158
3159	for(CompositeDestination::ChildList::iterator citr=destinationGraph->_children.begin();
3160	citr!=destinationGraph->_children.end();
3161	++citr)
3162	{
3163	(*citr)->_maxVisibleDistance = cutOffDistance;
3164	}
3165
3166	}
3167	else if (needToDivideX)
3168	{
3169	osg::notify(osg::INFO)<<"Need to Divide X only"<<std::endl;
3170
3171	// create two tiles.
3172	osg::BoundingBox left(extents.xMin(),extents.yMin(),extents.zMin(),xCenter,extents.yMax(),extents.zMax());
3173	osg::BoundingBox right(xCenter,extents.yMin(),extents.zMin(),extents.xMax(),extents.yMax(),extents.zMax());
3174
3175	destinationGraph->_children.push_back(createDestinationGraph(destinationGraph,
3176	cs,
3177	left,
3178	maxImageSize,
3179	maxTerrainSize,
3180	newLevel,
3181	newX,
3182	newY,
3183	maxNumLevels));
3184
3185	destinationGraph->_children.push_back(createDestinationGraph(destinationGraph,
3186	cs,
3187	right,
3188	maxImageSize,
3189	maxTerrainSize,
3190	newLevel,
3191	newX+1,
3192	newY,
3193	maxNumLevels));
3194
3195
3196	// Set all there max distance to the same value to ensure the same LOD bining.
3197	float cutOffDistance = destinationGraph->_maxVisibleDistance*0.5f;
3198
3199	for(CompositeDestination::ChildList::iterator citr=destinationGraph->_children.begin();
3200	citr!=destinationGraph->_children.end();
3201	++citr)
3202	{
3203	(*citr)->_maxVisibleDistance = cutOffDistance;
3204	}
3205
3206	}
3207	else if (needToDivideY)
3208	{
3209	osg::notify(osg::INFO)<<"Need to Divide Y only"<<std::endl;
3210
3211	// create two tiles.
3212	osg::BoundingBox top(extents.xMin(),yCenter,extents.zMin(),extents.xMax(),extents.yMax(),extents.zMax());
3213	osg::BoundingBox bottom(extents.xMin(),extents.yMin(),extents.zMin(),extents.xMax(),yCenter,extents.zMax());
3214
3215	destinationGraph->_children.push_back(createDestinationGraph(destinationGraph,
3216	cs,
3217	bottom,
3218	maxImageSize,
3219	maxTerrainSize,
3220	newLevel,
3221	newX,
3222	newY,
3223	maxNumLevels));
3224
3225	destinationGraph->_children.push_back(createDestinationGraph(destinationGraph,
3226	cs,
3227	top,
3228	maxImageSize,
3229	maxTerrainSize,
3230	newLevel,
3231	newX,
3232	newY+1,
3233	maxNumLevels));
3234
3235	// Set all there max distance to the same value to ensure the same LOD bining.
3236	float cutOffDistance = destinationGraph->_maxVisibleDistance*0.5f;
3237
3238	for(CompositeDestination::ChildList::iterator citr=destinationGraph->_children.begin();
3239	citr!=destinationGraph->_children.end();
3240	++citr)
3241	{
3242	(*citr)->_maxVisibleDistance = cutOffDistance;
3243	}
3244
3245	}
3246	else
3247	{
3248	osg::notify(osg::INFO)<<"No Need to Divide"<<std::endl;
3249	}
3250	}
3251
3252	return destinationGraph;
3253	}
3254
3255	void DataSet::computeDestinationGraphFromSources(unsigned int numLevels)
3256	{
3257	if (!_sourceGraph) return;
3258
3259	// ensure we have a valid coordinate system
3260	if (!_destinationCoordinateSystem)
3261	{
3262	for(CompositeSource::source_iterator itr(_sourceGraph.get());itr.valid();++itr)
3263	{
3264	SourceData* sd = (*itr)->getSourceData();
3265	if (sd)
3266	{
3267	if (sd->_cs.valid())
3268	{
3269	_destinationCoordinateSystem = sd->_cs;
3270	osg::notify(osg::INFO)<<"Setting coordinate system to "<<_destinationCoordinateSystem->getCoordinateSystem()<<std::endl;
3271	break;
3272	}
3273	}
3274	}
3275	}
3276
3277	if (!_intermediateCoordinateSystem)
3278	{
3279	CoordinateSystemType cst = getCoordinateSystemType(_destinationCoordinateSystem.get());
3280
3281	osg::notify(osg::INFO)<<"new DataSet::createDestination()"<<std::endl;
3282	if (cst!=GEOGRAPHIC && getConvertFromGeographicToGeocentric())
3283	{
3284	// need to use the geocentric coordinate system as a base for creating an geographic intermediate
3285	// coordinate system.
3286	OGRSpatialReference oSRS;
3287
3288	char *pszWKT = NULL;
3289	oSRS.SetWellKnownGeogCS( "WGS84" );
3290	oSRS.exportToWkt( &pszWKT );
3291
3292	setIntermediateCoordinateSystem(pszWKT);
3293	}
3294	else
3295	{
3296	_intermediateCoordinateSystem = _destinationCoordinateSystem;
3297	}
3298	}
3299
3300
3301	CoordinateSystemType destinateCoordSytemType = getCoordinateSystemType(_destinationCoordinateSystem.get());
3302	if (destinateCoordSytemType==GEOGRAPHIC && !getConvertFromGeographicToGeocentric())
3303	{
3304	// convert elevation into degrees.
3305	setVerticalScale(1.0f/111319.0f);
3306	}
3307
3308	// get the extents of the sources and
3309	osg::BoundingBox extents(_extents);
3310	if (!extents.valid())
3311	{
3312	for(CompositeSource::source_iterator itr(_sourceGraph.get());itr.valid();++itr)
3313	{
3314	SourceData* sd = (*itr)->getSourceData();
3315	if (sd)
3316	{
3317	osg::BoundingBox local_extents(sd->getExtents(_intermediateCoordinateSystem.get()));
3318	osg::notify(osg::INFO)<<"local_extents = xMin()"<<local_extents.xMin()<<" "<<local_extents.xMax()<<std::endl;
3319	osg::notify(osg::INFO)<<" yMin()"<<local_extents.yMin()<<" "<<local_extents.yMax()<<std::endl;
3320	extents.expandBy(local_extents);
3321	}
3322	}
3323	}
3324
3325	osg::notify(osg::INFO)<<"extents = xMin()"<<extents.xMin()<<" "<<extents.xMax()<<std::endl;
3326	osg::notify(osg::INFO)<<" yMin()"<<extents.yMin()<<" "<<extents.yMax()<<std::endl;
3327
3328	// then create the destinate graph accordingly.
3329	_destinationGraph = createDestinationGraph(0,
3330	_intermediateCoordinateSystem.get(),
3331	extents,
3332	_maximumTileImageSize,
3333	_maximumTileTerrainSize,
3334	0,
3335	0,
3336	0,
3337	numLevels);
3338
3339	// now traverse the destination graph to build neighbours.
3340	_destinationGraph->computeNeighboursFromQuadMap();
3341
3342	}
3343
3344	template<class T>
3345	struct DerefLessFunctor
3346	{
3347	bool operator () (const T& lhs, const T& rhs)
3348	{
3349	return (lhs.valid() && rhs.valid()) && (lhs->getSortValue() > rhs->getSortValue());
3350	}
3351	};
3352
3353	void DataSet::CompositeSource::setSortValueFromSourceDataResolution()
3354	{
3355	for(SourceList::iterator sitr=_sourceList.begin();sitr!=_sourceList.end();++sitr)
3356	{
3357	(*sitr)->setSortValueFromSourceDataResolution();
3358	}
3359
3360
3361	for(ChildList::iterator citr=_children.begin();citr!=_children.end();++citr)
3362	{
3363	(*citr)->setSortValueFromSourceDataResolution();
3364	}
3365	}
3366
3367	void DataSet::CompositeSource::sort()
3368	{
3369	// sort the sources.
3370	std::sort(_sourceList.begin(),_sourceList.end(),DerefLessFunctor< osg::ref_ptr<Source> >());
3371
3372	// sort the composite sources internal data
3373	for(ChildList::iterator itr=_children.begin();itr!=_children.end();++itr)
3374	{
3375	if (itr->valid()) (*itr)->sort();
3376	}
3377	}
3378
3379	void DataSet::updateSourcesForDestinationGraphNeeds()
3380	{
3381	if (!_destinationGraph \|\| !_sourceGraph) return;
3382
3383
3384	std::string temporyFilePrefix("temporaryfile_");
3385
3386	// compute the resolutions of the source that are required.
3387	{
3388	_destinationGraph->addRequiredResolutions(_sourceGraph.get());
3389
3390
3391	for(CompositeSource::source_iterator sitr(_sourceGraph.get());sitr.valid();++sitr)
3392	{
3393	Source* source = sitr->get();
3394	if (source)
3395	{
3396	osg::notify(osg::INFO)<<"Source File "<<source->getFileName()<<std::endl;
3397
3398
3399	const Source::ResolutionList& resolutions = source->getRequiredResolutions();
3400	osg::notify(osg::INFO)<<" resolutions.size() "<<resolutions.size()<<std::endl;
3401	osg::notify(osg::INFO)<<" { "<<std::endl;
3402	Source::ResolutionList::const_iterator itr;
3403	for(itr=resolutions.begin();
3404	itr!=resolutions.end();
3405	++itr)
3406	{
3407	osg::notify(osg::INFO)<<" resX="<<itr->_resX<<" resY="<<itr->_resY<<std::endl;
3408	}
3409	osg::notify(osg::INFO)<<" } "<<std::endl;
3410
3411	source->consolodateRequiredResolutions();
3412
3413	osg::notify(osg::INFO)<<" consolodated resolutions.size() "<<resolutions.size()<<std::endl;
3414	osg::notify(osg::INFO)<<" consolodated { "<<std::endl;
3415	for(itr=resolutions.begin();
3416	itr!=resolutions.end();
3417	++itr)
3418	{
3419	osg::notify(osg::INFO)<<" resX="<<itr->_resX<<" resY="<<itr->_resY<<std::endl;
3420	}
3421	osg::notify(osg::INFO)<<" } "<<std::endl;
3422	}
3423
3424	}
3425
3426
3427	}
3428
3429	// do standardisation of coordinates systems.
3430	// do any reprojection if required.
3431	{
3432	for(CompositeSource::source_iterator itr(_sourceGraph.get());itr.valid();++itr)
3433	{
3434	Source* source = itr->get();
3435	if (source && source->needReproject(_intermediateCoordinateSystem.get()))
3436	{
3437	// do the reprojection to a tempory file.
3438	std::string newFileName = temporyFilePrefix + osgDB::getStrippedName(source->getFileName()) + ".tif";
3439
3440	Source* newSource = source->doReproject(newFileName,_intermediateCoordinateSystem.get());
3441
3442	// replace old source by new one.
3443	if (newSource) *itr = newSource;
3444	else
3445	{
3446	osg::notify(osg::WARN)<<"Failed to reproject"<<source->getFileName()<<std::endl;
3447	*itr = 0;
3448	}
3449	}
3450	}
3451	}
3452
3453	// do sampling of data to required values.
3454	{
3455	for(CompositeSource::source_iterator itr(_sourceGraph.get());itr.valid();++itr)
3456	{
3457	Source* source = itr->get();
3458	if (source) source->buildOverviews();
3459	}
3460	}
3461
3462	// sort the sources so that the lowest res tiles are drawn first.
3463	{
3464	for(CompositeSource::source_iterator itr(_sourceGraph.get());itr.valid();++itr)
3465	{
3466	Source* source = itr->get();
3467	if (source)
3468	{
3469	source->setSortValueFromSourceDataResolution();
3470	osg::notify(osg::INFO)<<"sort "<<source->getFileName()<<" value "<<source->getSortValue()<<std::endl;
3471	}
3472
3473	}
3474
3475	// sort them so highest sortValue is first.
3476
3477	_sourceGraph->sort();
3478	}
3479
3480	osg::notify(osg::INFO)<<"Using source_lod_iterator itr"<<std::endl;
3481	for(CompositeSource::source_lod_iterator csitr(_sourceGraph.get(),CompositeSource::LODSourceAdvancer(0.0));csitr.valid();++csitr)
3482	{
3483	Source* source = csitr->get();
3484	if (source)
3485	{
3486	osg::notify(osg::INFO)<<" LOD "<<(*csitr)->getFileName()<<std::endl;
3487	}
3488	}
3489	osg::notify(osg::INFO)<<"End of Using Source Iterator itr"<<std::endl;
3490
3491	}
3492
3493	void DataSet::populateDestinationGraphFromSources()
3494	{
3495	if (!_destinationGraph \|\| !_sourceGraph) return;
3496
3497	osg::notify(osg::NOTICE)<<std::endl<<"started DataSet::populateDestinationGraphFromSources)"<<std::endl;
3498
3499	if (_databaseType==LOD_DATABASE)
3500	{
3501
3502	// for each DestinationTile populate it.
3503	_destinationGraph->readFrom(_sourceGraph.get());
3504
3505	// for each DestinationTile equalize the boundaries so they all fit each other without gaps.
3506	_destinationGraph->equalizeBoundaries();
3507
3508
3509	}
3510	else
3511	{
3512	// for each level
3513	// compute x and y range
3514	// from top row down to bottom row equalize boundairies a write out
3515	}
3516	osg::notify(osg::NOTICE)<<"completed DataSet::populateDestinationGraphFromSources)"<<std::endl;
3517	}
3518
3519
3520	void DataSet::_readRow(Row& row)
3521	{
3522	osg::notify(osg::NOTICE)<<"_readRow "<<row.size()<<std::endl;
3523	for(Row::iterator citr=row.begin();
3524	citr!=row.end();
3525	++citr)
3526	{
3527	CompositeDestination* cd = citr->second;
3528	for(CompositeDestination::TileList::iterator titr=cd->_tiles.begin();
3529	titr!=cd->_tiles.end();
3530	++titr)
3531	{
3532	DestinationTile* tile = titr->get();
3533	osg::notify(osg::NOTICE)<<" reading tile level="<<tile->_level<<" X="<<tile->_tileX<<" Y="<<tile->_tileY<<std::endl;
3534	tile->readFrom(_sourceGraph.get());
3535	}
3536	}
3537	}
3538
3539	void DataSet::_equalizeRow(Row& row)
3540	{
3541	osg::notify(osg::NOTICE)<<"_equalizeRow "<<row.size()<<std::endl;
3542	for(Row::iterator citr=row.begin();
3543	citr!=row.end();
3544	++citr)
3545	{
3546	CompositeDestination* cd = citr->second;
3547	for(CompositeDestination::TileList::iterator titr=cd->_tiles.begin();
3548	titr!=cd->_tiles.end();
3549	++titr)
3550	{
3551	DestinationTile* tile = titr->get();
3552	osg::notify(osg::NOTICE)<<" equalizing tile level="<<tile->_level<<" X="<<tile->_tileX<<" Y="<<tile->_tileY<<std::endl;
3553	tile->equalizeBoundaries();
3554	tile->setTileComplete(true);
3555	}
3556	}
3557	}
3558
3559	void DataSet::_writeRow(Row& row)
3560	{
3561	osg::notify(osg::NOTICE)<<"_writeRow "<<row.size()<<std::endl;
3562	for(Row::iterator citr=row.begin();
3563	citr!=row.end();
3564	++citr)
3565	{
3566	CompositeDestination* cd = citr->second;
3567	CompositeDestination* parent = cd->_parent;
3568
3569	if (parent)
3570	{
3571	if (!parent->getSubTilesGenerated() && parent->areSubTilesComplete())
3572	{
3573	osg::ref_ptr<osg::Node> node = parent->createSubTileScene();
3574	std::string filename = parent->getSubTileName();
3575	if (node.valid())
3576	{
3577	osg::notify(osg::NOTICE)<<" writeSubTile filename="<<filename<<std::endl;
3578	osgDB::writeNodeFile(*node,filename);
3579	parent->setSubTilesGenerated(true);
3580	parent->unrefSubTileData();
3581	}
3582	else
3583	{
3584	osg::notify(osg::WARN)<<" failed to writeSubTile node for tile, filename="<<filename<<std::endl;
3585	}
3586	}
3587	}
3588	else
3589	{
3590	osg::ref_ptr<osg::Node> node = cd->createPagedLODScene();
3591
3592	if (_decorateWithCoordinateSystemNode)
3593	{
3594	node = decorateWithCoordinateSystemNode(node.get());
3595	}
3596
3597
3598	if (!_comment.empty())
3599	{
3600	node->addDescription(_comment);
3601	}
3602
3603	//std::string filename = cd->_name + _tileExtension;
3604	std::string filename = _tileBasename+_tileExtension;
3605
3606	if (node.valid())
3607	{
3608	osg::notify(osg::NOTICE)<<" writeNodeFile = "<<cd->_level<<" X="<<cd->_tileX<<" Y="<<cd->_tileY<<" filename="<<filename<<std::endl;
3609	osgDB::writeNodeFile(*node,filename);
3610	}
3611	else
3612	{
3613	osg::notify(osg::WARN)<<" faild to write node for tile = "<<cd->_level<<" X="<<cd->_tileX<<" Y="<<cd->_tileY<<" filename="<<filename<<std::endl;
3614	}
3615	}
3616	}
3617	}
3618
3619	void DataSet::createDestination(unsigned int numLevels)
3620	{
3621	osg::notify(osg::NOTICE)<<std::endl<<"started DataSet::createDestination("<<numLevels<<")"<<std::endl;
3622
3623	computeDestinationGraphFromSources(numLevels);
3624
3625	updateSourcesForDestinationGraphNeeds();
3626
3627	osg::notify(osg::NOTICE)<<"completed DataSet::createDestination("<<numLevels<<")"<<std::endl;
3628
3629	}
3630
3631	osg::Node* DataSet::decorateWithCoordinateSystemNode(osg::Node* subgraph)
3632	{
3633	// don't decorate if no coord system is set.
3634	if (_destinationCoordinateSystem->getCoordinateSystem().empty())
3635	return subgraph;
3636
3637	osg::CoordinateSystemNode* csn = new osg::CoordinateSystemNode(
3638	_destinationCoordinateSystem->getFormat(),
3639	_destinationCoordinateSystem->getCoordinateSystem());
3640
3641	// set the ellipsoid model if geocentric coords are used.
3642	if (getConvertFromGeographicToGeocentric()) csn->setEllipsoidModel(getEllipsoidModel());
3643
3644	// add the a subgraph.
3645	csn->addChild(subgraph);
3646
3647	return csn;
3648	}
3649
3650
3651	void DataSet::writeDestination()
3652	{
3653	if (_destinationGraph.valid())
3654	{
3655	std::string filename = _tileBasename+_tileExtension;
3656	osg::notify(osg::NOTICE)<<std::endl<<"started DataSet::writeDestination("<<filename<<")"<<std::endl;
3657
3658	if (_databaseType==LOD_DATABASE)
3659	{
3660	populateDestinationGraphFromSources();
3661	_rootNode = _destinationGraph->createScene();
3662
3663	if (_decorateWithCoordinateSystemNode)
3664	{
3665	_rootNode = decorateWithCoordinateSystemNode(_rootNode.get());
3666	}
3667
3668	if (!_comment.empty())
3669	{
3670	_rootNode->addDescription(_comment);
3671	}
3672
3673	osgDB::writeNodeFile(*_rootNode,filename);
3674	}
3675	else // _databaseType==PagedLOD_DATABASE
3676	{
3677	// for each level build read and write the rows.
3678	for(QuadMap::iterator qitr=_quadMap.begin();
3679	qitr!=_quadMap.end();
3680	++qitr)
3681	{
3682	Level& level = qitr->second;
3683
3684	// skip is level is empty.
3685	if (level.empty()) continue;
3686
3687	osg::notify(osg::INFO)<<"New level"<<std::endl;
3688
3689	Level::iterator prev_itr = level.begin();
3690	_readRow(prev_itr->second);
3691	Level::iterator curr_itr = prev_itr;
3692	++curr_itr;
3693	for(;
3694	curr_itr!=level.end();
3695	++curr_itr)
3696	{
3697	_readRow(curr_itr->second);
3698
3699	_equalizeRow(prev_itr->second);
3700	_writeRow(prev_itr->second);
3701
3702	prev_itr = curr_itr;
3703	}
3704
3705	_equalizeRow(prev_itr->second);
3706	_writeRow(prev_itr->second);
3707	}
3708	}
3709	osg::notify(osg::NOTICE)<<"completed DataSet::writeDestination("<<filename<<")"<<std::endl;
3710
3711	}
3712	else
3713	{
3714	osg::notify(osg::WARN)<<"Error: no scene graph to output, no file written."<<std::endl;
3715	}
3716
3717	}

Note: See TracBrowser for help on using the browser.

Context Navigation

root/OpenSceneGraph/trunk/src/osgTerrain/DataSet.cpp @ 3309

Download in other formats: