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

Revision 3309, 135.2 kB (checked in by robert, 10 years ago)

Added protection to prevent crashes on calls when no data is set up

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