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

Revision 3292, 135.1 kB (checked in by robert, 10 years ago)

Removed computeMipMaps call

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