1 #include "Geographic"
2 
3 #include <algorithm>
4 #include <cmath>
5 #include <iterator>
6 #include <vector>
7 
8 #include <osg/ClusterCullingCallback>
9 #include <osg/CoordinateSystemNode>
10 #include <osg/Math>
11 #include <osg/NodeCallback>
12 #include <osg/NodeVisitor>
13 #include <osg/Texture2D>
14 
15 #include <osgEarth/ImageUtils>
16 #include <osgEarth/Notify>
17 #include <osgEarth/VerticalSpatialReference>
18 #include <osgEarth/TaskService>
19 
20 #include "GeoPatch"
21 #include "MultiArray"
22 
23 namespace seamless
24 {
25 using namespace std;
26 using namespace osg;
27 using namespace osgEarth;
28 
29 typedef multi_array_ref<Vec3f, Vec3Array, 2> PatchArray;
30 
31 Geographic::Geographic(const Map* map,
32                        const osgEarth::Drivers::SeamlessOptions& options)
33     : PatchSet(options, new PatchOptions), _profile(new EulerProfile),
34       _eModel(new EllipsoidModel)
35 {
36     setPrecisionFactor(8);
37     setMap(map);
38     {
39         int maxLevel = 0;
40         const ElevationLayerVector& elevations = _mapf->elevationLayers();
41         for (ElevationLayerVector::const_iterator itr = elevations.begin(),
42                  end = elevations.end();
43              itr != end;
44              ++itr)
45         {
46             const TerrainLayerOptions& options
47                 = (*itr)->getTerrainLayerOptions();
48             if (options.maxLevel().isSet()
49                 && options.maxLevel().get() > maxLevel)
50                 maxLevel = options.maxLevel().get();
51         }
52         if (maxLevel > 0)
53             setMaxLevel(maxLevel);
54 
55     }
56     int serviceThreads = computeLoadingThreads(_options.loadingPolicy().get());
57     _hfService = new TaskService("Height Field Service", serviceThreads);
58     _imageService = new TaskService("Image Service", serviceThreads);
59 }
60 
61 Geographic::Geographic(const Geographic& rhs, const osg::CopyOp& copyop)
62     : PatchSet(rhs, copyop),
63       _profile(static_cast<EulerProfile*>(copyop(rhs._profile.get()))),
64       _eModel(static_cast<EllipsoidModel*>(copyop(rhs._eModel.get()))),
65       _hfService(rhs._hfService), _imageService(rhs._imageService)
66 {
67 }
68 
69 Geographic::~Geographic()
70 {
71 }
72 
73 Node* Geographic::createPatchSetGraph(const std::string& filename)
74 {
75     CoordinateSystemNode* csn = new CoordinateSystemNode;
76     // Should these values come from the map profile?
77     csn->setCoordinateSystem("EPSG:4326");
78     csn->setFormat("WKT");
79     csn->setEllipsoidModel(_eModel.get());
80     for (int face = 0; face < 6; ++face)
81     {
82         double x = 0.0, y = 0.0;
83         euler::faceToCube(x, y, face);
84         PatchOptions* poptions = static_cast<PatchOptions*>(
85             osg::clone(getPatchOptionsPrototype()));
86         poptions->setPatchSet(this);
87         poptions->setTileKey(_profile->createTileKey(x, y, 2));
88         Node* node = createPatchGroup("foobar.osgearth_engine_seamless_patch",
89                                       poptions);
90         csn->addChild(node);
91     }
92     return csn;
93 }
94 
95 namespace
96 {
97 
98 GeoHeightField
99 mergeHeightFields(const GeoExtent& targetExtent, const GeoHeightFieldVector& hfs)
100 {
101     if (hfs.size() != 4)
102     {
103         OE_FATAL << "mergeHeightFields expected 4 height fields\n";
104         return GeoHeightField();
105     }
106     // List is in tile subkey quadrant order.
107     // Assume the height fields all have the same dimensions
108     unsigned targetCols = hfs[0].getHeightField()->getNumColumns() * 2 - 1;
109     unsigned targetRows = hfs[0].getHeightField()->getNumRows() * 2 - 1;
110     HeightField* targethf = new HeightField;
111     targethf->allocate(targetCols, targetRows);
112     GeoHeightField geo(targethf, targetExtent, 0);
113     for (int i = 0; i < 4; ++i)
114     {
115         const HeightField* src = hfs[i].getHeightField();
116         unsigned targetColumn
117             = floor((hfs[i].getExtent().xMin() - targetExtent.xMin())
118                     / targetExtent.width() * (targetCols - 1) + .5);
119         unsigned targetRow
120             = floor((hfs[i].getExtent().yMin() - targetExtent.yMin())
121                     / targetExtent.height() * (targetRows - 1) + .5);
122         for (unsigned sj = 0, tj = targetRow;
123              sj < src->getNumRows() && tj < targetRows;
124              ++sj, ++tj)
125         {
126             for (unsigned si = 0, ti = targetColumn;
127              si < src->getNumColumns() && ti < targetCols;
128              ++si, ++ti)
129                 targethf->setHeight(ti, tj, src->getHeight(si, sj));
130         }
131     }
132     return geo;
133 }
134 
135 GeoImage
136 mergeImages(const GeoExtent& targetExtent, const GeoImageVector& imgs)
137 {
138     Image* targetImage = new Image;
139     const Image* proto = imgs[0].getImage();
140     targetImage->setInternalTextureFormat(proto->getInternalTextureFormat());
141     int numRows = proto->s() * 2;
142     int numCols = proto->t() * 2;
143     targetImage->allocateImage(numRows, numCols, proto->r(),
144                                proto->getPixelFormat(), proto->getDataType(),
145                                proto->getPacking());
146     for (GeoImageVector::const_iterator itr = imgs.begin(),
147              end = imgs.end();
148          itr != end;
149          ++itr)
150     {
151         const GeoExtent& srcExtent = itr->getExtent();
152         int dstx
153             = floor((srcExtent.xMin() - targetExtent.xMin()) / targetExtent.width()
154                     * numCols + .5);
155         int dsty
156             = floor((srcExtent.yMin() - targetExtent.yMin()) / targetExtent.height()
157                     * numRows + .5);
158         ImageUtils::copyAsSubImage(itr->getImage(), targetImage,
159                                    dstx, dsty);
160     }
161     return GeoImage(targetImage, targetExtent);
162 }
163 }
164 
165 // Create vertex arrays from the height field for a patch and install
166 // them in the patch.
167 void expandHeights(Geographic* gpatchset, const TileKey& key,
168                    const GeoHeightField& hf, Vec3Array* verts,
169                    Vec3Array* normals)
170 {
171     int resolution = gpatchset->getResolution();
172     const GeoExtent& patchExtent = key.getExtent();
173     double centx, centy;
174     patchExtent.getCentroid(centx, centy);
175     Vec3d patchCenter = gpatchset->toModel(centx, centy, 0);
176     const SpatialReference* srs = key.getProfile()->getSRS();
177     const SpatialReference* geoSrs = srs->getGeographicSRS();
178     // Populate cell
179     ref_ptr<Patch::Data> data = new Patch::Data;
180     int patchDim = resolution + 1;
181     double xInc = (patchExtent.xMax() - patchExtent.xMin()) / resolution;
182     double yInc = (patchExtent.yMax() - patchExtent.yMin()) / resolution;
183     const EllipsoidModel* eModel = gpatchset->getEllipsoidModel();
184     const float verticalScale = gpatchset->getVerticalScale();
185     PatchArray mverts(*verts, patchDim);
186     for (int j = 0; j < patchDim; ++j)
187     {
188         for (int i = 0; i < patchDim; i++)
189         {
190             Vec2d cubeCoord(patchExtent.xMin() + i * xInc,
191                             patchExtent.yMin() + j * yInc);
192             double lon, lat;
193             srs->transform(cubeCoord.x(), cubeCoord.y(), geoSrs, lon, lat);
194             float elevation;
195 
196             bool found = hf.getElevation(srs, cubeCoord.x(), cubeCoord.y(),
197                                          INTERP_BILINEAR, 0, elevation);
198             // Into ec coordinates
199             if (!found)
200             {
201                 OE_WARN << "Couldn't find height sample for cube coordinates "
202                         << cubeCoord.x() << ", " << cubeCoord.y()
203                         << " (lon lat " << lon << ", " << lat << ")\n";
204                 continue;
205             }
206             elevation *= verticalScale;
207             Vec3d coord;
208             eModel->convertLatLongHeightToXYZ(
209                 DegreesToRadians(lat), DegreesToRadians(lon), elevation,
210                 coord.x(), coord.y(), coord.z());
211             mverts[j][i] = coord - patchCenter;
212             if (fabs(mverts[j][i].z()) > 6000000)
213                 OE_WARN << "found huge coordinate.\n";
214         }
215     }
216     // Normals. Average the normals of the triangles around the sample
217     // point. We're not following the actual tessallation of the grid.
218     for (int j = 0; j < patchDim; ++j)
219     {
220         for (int i = 0; i < patchDim; i++)
221         {
222             const Vec3& pt = (*verts)[j * patchDim + i];
223             // A cross of points.
224             Vec3 delta[4];      // Initialized to zero vectors
225             for (int k = 0; k < 2; ++k)
226             {
227                 int gridx = i + 2 * k - 1;
228                 if (gridx < patchDim && gridx >= 0)
229                     delta[2 * k] = (*verts)[j * patchDim + gridx] - pt;
230             }
231             for (int k = 0; k < 2; ++k)
232             {
233                 int gridy = j + 2 * k - 1;
234                 if (gridy < patchDim && gridy >= 0)
235                     delta[2 * k + 1] = (*verts)[gridy * patchDim + i] - pt;
236             }
237             Vec3 normal;
238             for (int k = 1; k <= 4; ++k)
239             {
240                 int v1 = k - 1, v2 = k % 4;
241                 // If One or both of the deltas are 0, then the cross
242                 // product is 0 and won't contribute to the average.
243                 normal += delta[v1] ^ delta[v2];
244             }
245             normal.normalize();
246             (*normals)[j * patchDim + i] = normal;
247         }
248     }
249 }
250 
251 void installHeightField(GeoPatch* patch, const TileKey& key,
252                         const GeoHeightField& hf)
253 {
254     Geographic* gpatchset = patch->getGeographic();
255     int resolution = gpatchset->getResolution();
256     // Populate cell
257     int patchDim = resolution + 1;
258     Vec3Array* verts = new Vec3Array(patchDim * patchDim);
259     verts->setDataVariance(Object::DYNAMIC);
260     Vec3Array* normals = new Vec3Array(patchDim * patchDim);
261     normals->setDataVariance(Object::DYNAMIC);
262     Vec2Array* texCoords = new Vec2Array(patchDim * patchDim);
263     expandHeights(gpatchset, key, hf, verts, normals);
264     const float resinv = 1.0f / static_cast<float>(resolution);
265     for (int j = 0; j < patchDim; ++j)
266     {
267         for (int i = 0; i < patchDim; i++)
268         {
269             (*texCoords)[j * patchDim +i] = Vec2(i * resinv, j * resinv);
270         }
271     }
272     // Construct the patch and its transform.
273     ref_ptr<Patch::Data> data = new Patch::Data;
274     data->vertexData.array = verts;
275     data->vertexData.binding = Geometry::BIND_PER_VERTEX;
276     data->normalData.array = normals;
277     data->normalData.binding = Geometry::BIND_PER_VERTEX;
278     Vec4Array* colors = new Vec4Array(1);
279     (*colors)[0] = Vec4(1.0, 1.0, 1.0, 1.0);
280     data->colorData.array = colors;
281     data->colorData.binding = Geometry::BIND_OVERALL;
282     data->texCoordList
283         .push_back(Geometry::ArrayData(texCoords, Geometry::BIND_PER_VERTEX));
284     patch->setData(data);
285 }
286 
287 // Create a patch and the transform that places it in the
288 // world. Install a height field if one is given.
289 MatrixTransform* createPatchAux(Geographic* gpatchset,
290                                 const TileKey& key,
291                                 const GeoHeightField& hf)
292 {
293     GeoPatch* patch = new GeoPatch(key);
294     patch->setGeographic(gpatchset);
295     const GeoExtent& patchExtent = key.getExtent();
296     double centx, centy;
297     patchExtent.getCentroid(centx, centy);
298     Vec3d patchCenter = gpatchset->toModel(centx, centy, 0);
299     Matrixd patchMat = Matrixd::translate(patchCenter);
300     installHeightField(patch, key, hf);
301     MatrixTransform* result = new MatrixTransform;
302     result->addChild(patch);
303     result->setMatrix(patchMat);
304     return result;
305 }
306 
307 
308 namespace
309 {
310 // Get a height field from the map, or an empty one if there is no
311 // data for this tile.
312 GeoHeightField getGeoHeightField(MapFrame& mapf, const TileKey& key,
313                                  int resolution)
314 {
315     osg::ref_ptr<HeightField> hf;
316     mapf.getHeightField(key, true, hf, 0L, INTERP_BILINEAR);
317     if  (!hf)
318         hf = key.getProfile()->getVerticalSRS()
319             ->createReferenceHeightField(key.getExtent(),
320                                          resolution + 1, resolution + 1);
321     return GeoHeightField(hf, key.getExtent(),
322                           key.getProfile()->getVerticalSRS());
323 }
324 
325 // Split up patch keys that cross the Date Line. The only patches
326 // that do are the the equatorial face with center at (-180, 0),
327 // and the poles faces.
328 
329 inline bool crossesDateLine(const TileKey& key)
330 {
331     int face = EulerProfile::getFace(key);
332     const GeoExtent& keyExtent = key.getExtent();
333     return ((face == 2 || face == 4 || face == 5)
334             && keyExtent.xMax() - keyExtent.xMin() > .5);
335 }
336 
337 struct HeightFieldRequest : public TaskRequest
338 {
339     HeightFieldRequest(Geographic* gpatchset, const TileKey& key)
340 
341         : _gpatchset(gpatchset), _key(key), _mapf(gpatchset->getMapFrame())
342     {
343     }
344     void operator()(ProgressCallback* progress)
345     {
346         const Map* map = _gpatchset->getMap();
347         int resolution = _gpatchset->getResolution();
348         GeoHeightField hf;
349         if (crossesDateLine(_key))
350         {
351             GeoHeightFieldVector hfs;
352             for (int child = 0; child < 4; ++child)
353             {
354                 TileKey subCubeKey = _key.createChildKey(child);
355                 hfs.push_back(getGeoHeightField(_mapf, subCubeKey, resolution));
356             }
357             hf = mergeHeightFields(_key.getExtent(), hfs);
358         }
359         else
360         {
361             hf = getGeoHeightField(_mapf, _key, resolution);
362         }
363         int patchDim = resolution + 1;
364         Vec3Array* verts = new Vec3Array(patchDim * patchDim);
365         _result = verts;
366         _normalResult = new Vec3Array(patchDim * patchDim);
367         expandHeights(_gpatchset.get(), _key, hf,
368                       verts, _normalResult.get());
369     }
370     ref_ptr<Geographic> _gpatchset;
371     TileKey _key;
372     // vertices are in _result;
373     ref_ptr<Vec3Array> _normalResult;
374     MapFrame _mapf;
375 };
376 
377 struct ImageRequest : public TaskRequest
378 {
379     ImageRequest(Geographic* gpatchset, const TileKey& key)
380         : _gpatchset(gpatchset), _key(key), _mapf(gpatchset->getMapFrame())
381     {
382     }
383 
384     void operator()(ProgressCallback* progress)
385     {
386         GeoImage gimage;
387         const ImageLayerVector& layers = _mapf.imageLayers();
388         if (crossesDateLine(_key))
389         {
390             GeoImageVector gis;
391             if (!layers.empty())
392             {
393                 for (int child = 0; child < 4; ++child)
394                 {
395                     TileKey subCubeKey = _key.createChildKey(child);
396                     gis.push_back(layers[0]->createImage(subCubeKey));
397                 }
398             }
399             if (!gis.empty())
400                 gimage = mergeImages(_key.getExtent(), gis);
401         }
402         else
403         {
404             if (!layers.empty())
405                 gimage = layers[0]->createImage(_key);
406         }
407         _result = gimage.getImage();
408     }
409     ref_ptr<Geographic> _gpatchset;
410     const TileKey _key;
411     MapFrame _mapf;
412 };
413 
414 // Update a patch node once map data is available
415 class GeoPatchUpdateCallback : public NodeCallback
416 {
417 public:
418     GeoPatchUpdateCallback() {}
419     GeoPatchUpdateCallback(HeightFieldRequest* hfRequest,
420                            ImageRequest* imageRequest)
421         : _hfRequest(hfRequest), _imageRequest(imageRequest)
422     {
423     }
424 
425     GeoPatchUpdateCallback(const GeoPatchUpdateCallback& nc,
426                            const CopyOp& copyop)
427         : NodeCallback(nc, copyop), _hfRequest(nc._hfRequest),
428           _imageRequest(nc._imageRequest)
429     {
430     }
431 
432     META_Object(seamless, GeoPatchUpdateCallback);
433 
434     virtual void operator()(Node* node, NodeVisitor* nv);
435 
436     ref_ptr<HeightFieldRequest> _hfRequest;
437     ref_ptr<ImageRequest> _imageRequest;
438 };
439 }
440 
441 Transform* Geographic::createPatch(const std::string& filename,
442                                    PatchOptions* poptions)
443 {
444     const TileKey patchKey = poptions->getTileKey();
445     // Dummy height field until data is available.
446     const VerticalSpatialReference* vsrs
447         = patchKey.getProfile()->getVerticalSRS();
448     ref_ptr<HeightField> hf = vsrs->createReferenceHeightField(
449         patchKey.getExtent(), _resolution + 1, _resolution + 1);
450     GeoHeightField ghf(hf.get(), patchKey.getExtent(), vsrs);
451     ref_ptr<MatrixTransform> transform = createPatchAux(this, patchKey, ghf);
452     GeoPatch* patch = dynamic_cast<GeoPatch*>(transform->getChild(0));
453     ref_ptr<HeightFieldRequest> hfr = new HeightFieldRequest(this, patchKey);
454     ref_ptr<ImageRequest> ir = new ImageRequest(this, patchKey);
455     patch->setUpdateCallback(new GeoPatchUpdateCallback(hfr.get(), ir.get()));
456     _hfService->add(hfr.get());
457     _imageService->add(ir.get());
458     return transform.release();
459 }
460 
461 namespace
462 {
463 ClusterCullingCallback*
464 createClusterCullingCallback(const Matrixd& transform, const Patch* patch,
465                              const EllipsoidModel* et)
466 {
467     //This code is a very slightly modified version of the
468     //DestinationTile::createClusterCullingCallback in
469     //VirtualPlanetBuilder.
470     double globe_radius =  et->getRadiusPolar();
471     Vec3 center_position(transform.getTrans());
472     Vec3 center_normal(center_position);
473     center_normal.normalize();
474 
475     // populate the vertex/normal/texcoord arrays from the grid.
476 
477     float min_dot_product = 1.0f;
478     float max_cluster_culling_height = 0.0f;
479     float max_cluster_culling_radius = 0.0f;
480     const Vec3Array* verts = static_cast<const Vec3Array*>(
481         patch->getData()->vertexData.array.get());
482     for (Vec3Array::const_iterator itr = verts->begin(), end = verts->end();
483          itr != end;
484          ++itr)
485     {
486         Vec3d dv = *itr;
487         Vec3d v = dv + center_position;
488         double lat, lon, height;
489 
490         et->convertXYZToLatLongHeight(v.x(), v.y(), v.z(),
491                                       lat, lon, height);
492 
493         double d = sqrt(dv.x()*dv.x() + dv.y()*dv.y() + dv.z()*dv.z());
494         double theta = acos(globe_radius / (globe_radius + fabs(height)));
495         double phi = 2.0 * asin (d*0.5 / globe_radius); // d/globe_radius;
496         double beta = theta + phi;
497         double sb = sin(beta);
498         double cb = cos(beta);
499         double cutoff = osg::PI_2 - 0.1;
500 
501         //log(osg::INFO,"theta="<<theta<<"\tphi="<<phi<<" beta "<<beta);
502         if (phi<cutoff && beta<cutoff)
503         {
504             float local_dot_product = -sb;
505             float local_m = globe_radius*( 1.0/ cb - 1.0);
506             float local_radius = static_cast<float>(globe_radius * sb / cb); // beta*globe_radius;
507             min_dot_product = osg::minimum(min_dot_product, local_dot_product);
508             max_cluster_culling_height = osg::maximum(max_cluster_culling_height,local_m);
509             max_cluster_culling_radius = osg::maximum(max_cluster_culling_radius,local_radius);
510         }
511         else
512         {
513             //log(osg::INFO,"Turning off cluster culling for wrap around tile.");
514             return 0;
515         }
516     }
517 
518     osg::ClusterCullingCallback* ccc = new osg::ClusterCullingCallback;
519 
520     ccc->set(center_position + center_normal*max_cluster_culling_height ,
521              center_normal,
522              min_dot_product,
523              max_cluster_culling_radius);
524 
525     return ccc;
526 }
527 }
528 
529 Node* Geographic::createPatchGroup(const string& filename,
530                                    PatchOptions* poptions)
531 {
532     Node* result = PatchSet::createPatchGroup(filename, poptions);
533     PatchGroup* pgroup = dynamic_cast<PatchGroup*>(result);
534     // Make a cluster culling callback
535     MatrixTransform* transform
536         = dynamic_cast<MatrixTransform*>(pgroup->getChild(0));
537     Patch* patch = dynamic_cast<Patch*>(transform->getChild(0));
538     ClusterCullingCallback* ccc
539         = createClusterCullingCallback(transform->getMatrix(), patch,
540                                        _eModel.get());
541     pgroup->setCullCallback(ccc);
542     return pgroup;
543 }
544 
545 Vec3d Geographic::toModel(double cubeX, double cubeY, double elevation)
546 {
547     double faceX = cubeX, faceY = cubeY;
548     int face;
549     euler::cubeToFace(faceX, faceY, face);
550     double lat_deg, lon_deg;
551     euler::faceCoordsToLatLon(faceX, faceY, face, lat_deg, lon_deg);
552     Vec3d result;
553     _eModel->convertLatLongHeightToXYZ(
554         DegreesToRadians(lat_deg), DegreesToRadians(lon_deg), elevation,
555         result.x(), result.y(), result.z());
556     return result;
557 }
558 
559 Node* Geographic::createChild(const PatchOptions* parentOptions, int childNum)
560 {
561     PatchOptions* poptions = static_cast<PatchOptions*>(
562         osg::clone(parentOptions));
563     poptions->setPatchLevel(parentOptions->getPatchLevel() + 1);
564     poptions->setTileKey(parentOptions->getTileKey().createChildKey(childNum));
565     return createPatchGroup("foobies.osgearth_engine_seamless_patch", poptions);
566 
567 }
568 
569 // A tile can be thought of lying between edges with integer
570 // coordinates at its LOD. With x going to the right and y going down,
571 // a tile between (tile_x, tile_y) and (tile_x + 1, tile_y + 1).
572 //
573 // edge order should be counter clockwise
574 
575 struct GridEdge
576 {
577     unsigned v[2][2];
578     unsigned lod;
579 };
580 
581 struct KeyIndex
582 {
583     KeyIndex() : lod(0), x(0), y(0) {}
584     KeyIndex(unsigned lod_, unsigned x_, unsigned y_)
585         : lod(lod_), x(x_), y(y_)
586     {
587     }
588     KeyIndex(const TileKey& key)
589         : lod(key.getLevelOfDetail()), x(key.getTileX()), y(key.getTileY())
590     {
591     }
592     bool operator==(const KeyIndex& rhs) const
593     {
594         return lod == rhs.lod && x == rhs.x && y == rhs.y;
595     }
596     unsigned lod;
597     unsigned x;
598     unsigned y;
599 };
600 
601 bool containsTile(const KeyIndex& parent, const KeyIndex& child)
602 {
603     if (parent.lod > child.lod)
604         return false;
605     if (parent.lod == child.lod)
606         return parent.x == child.x && parent.y == child.y;
607     int lodDiff = child.lod - parent.lod;
608     if (child.x >> lodDiff == parent.x && child.y >> lodDiff == parent.y)
609         return true;
610     else
611         return false;
612 }
613 
614 // assume that tile is at a lower or equal LOD than neighbor
615 bool isNeighborTile(const KeyIndex& tile, const KeyIndex& neighbor)
616 {
617     int lodDiff = neighbor.lod - tile.lod;
618     int lodMult = 1 << lodDiff;
619     unsigned tx = tile.x << lodDiff;
620     unsigned ty = tile.y << lodDiff;
621     if (tx == neighbor.x + 1 || tx + lodMult == neighbor.x)
622         return neighbor.y >= ty && neighbor.y + 1 <= ty + lodMult;
623     else if (ty == neighbor.y + 1 || ty + lodMult == neighbor.y)
624         return neighbor.x >= tx && neighbor.x + 1 <= tx + lodMult;
625     return false;
626 
627 }
628 
629 // Do tiles share a corner?
630 bool adjoinsTile(const KeyIndex& tile, const KeyIndex& neighbor)
631 {
632     int lodDiff = neighbor.lod - tile.lod;
633     int lodMult = 1 << lodDiff;
634     unsigned tx = tile.x << lodDiff;
635     unsigned ty = tile.y << lodDiff;
636     if (tx == neighbor.x + 1 || tx + lodMult == neighbor.x)
637         return ty == neighbor.y + 1 || ty + lodMult == neighbor.y;
638     else
639         return false;
640 }
641 
642 PatchGroup* findFaceRoot(GeoPatch* patch, NodePath& pathList)
643 {
644     // Get the patch's key
645     Group* parent = patch->getParent(0);
646     PatchGroup* parentPatch = dynamic_cast<PatchGroup*>(parent->getParent(0));
647     if (!parentPatch)
648         return 0;
649     PatchOptions* parentOptions = parentPatch->getOptions();
650     TileKey patchKey = parentOptions->getTileKey();
651     int x = patchKey.getTileX() >> (patchKey.getLevelOfDetail() - 2);
652     int y = patchKey.getTileY() >> (patchKey.getLevelOfDetail() - 2);
653 
654     for (NodePath::iterator itr = pathList.begin(), end = pathList.end();
655          itr != end;
656          ++itr)
657     {
658         PatchGroup* pg = dynamic_cast<PatchGroup*>(*itr);
659         if (pg)
660         {
661             PatchOptions* poptions = pg->getOptions();
662             if (poptions)
663             {
664                 TileKey key = poptions->getTileKey();
665                 if (key.getLevelOfDetail() == 2 && x == key.getTileX()
666                     && y == key.getTileY())
667                 return pg;
668             }
669         }
670     }
671     return 0;
672 }
673 
674 typedef vector_ref<Vec3f, Vec3Array> EdgeRef;
675 
676 EdgeRef makeEdgeRef(GeoPatch* gpatch, int edgeno, int mult)
677 {
678     Vec3Array* verts
679         = static_cast<Vec3Array*>(gpatch->getData()->vertexData.array.get());
680     int patchDim = gpatch->getPatchSet()->getResolution() + 1;
681     int shape = (patchDim - 1) / mult + 1;
682     switch(edgeno)
683     {
684     case 0:
685         return EdgeRef(*verts, mult, shape, 0);
686     case 1:
687         return EdgeRef(*verts, patchDim * mult, shape, patchDim - 1);
688     case 2:
689         return EdgeRef(*verts, mult, shape, (patchDim - 1) * patchDim);
690     case 3:
691         return EdgeRef(*verts, patchDim * mult, shape, 0);
692     default:
693         return EdgeRef(*verts, 0, 0, 0); // shouldn't happen
694     }
695 }
696 
697 struct ShareResult
698 {
699     ShareResult()
700         : numEdges(0)
701     {
702         for (int i = 0; i < 2; ++i)
703             tile1[i] = tile2[i] = -1;
704     }
705     int numEdges;
706     int tile1[2];
707     int tile2[2];
708 };
709 
710 // tile2 is at a higher LOD than tile1
711 ShareResult tilesShareEdges(const KeyIndex& tile1, const KeyIndex& tile2)
712 {
713     ShareResult result;
714     int lodDiff = tile2.lod - tile1.lod;
715     int x = tile1.x << lodDiff;
716     int xleft = (tile1.x + 1) << lodDiff;
717     int y = tile1.y << lodDiff;
718     int ybottom = (tile1.y + 1) << lodDiff;
719     if (tile2.x >= x && tile2.x + 1 <= xleft
720         && tile2.y >= y && tile2.y + 1 <= ybottom)
721     {
722         // tile1 contains tile2; do they share edges?
723         if (x == tile2.x)
724         {
725             result.tile1[0] = 3;
726             result.tile2[0] = 3;
727             result.numEdges = 1;
728         }
729         else if (xleft == tile2.x + 1)
730         {
731             result.tile1[0] = 1;
732             result.tile2[0] = 1;
733             result.numEdges = 1;
734         }
735         if (y == tile2.y)
736         {
737             result.tile1[result.numEdges] = 2;
738             result.tile2[result.numEdges] = 2;
739             result.numEdges++;
740         }
741         else if (ybottom == tile2.y + 1)
742         {
743             result.tile1[result.numEdges] = 0;
744             result.tile2[result.numEdges] = 0;
745             result.numEdges++;
746         }
747     }
748     else
749     {
750         // Tiles can share 1 edge.
751         if (x == tile2.x + 1)
752         {
753             result.tile1[0] = 3;
754             result.tile2[0] = 1;
755             result.numEdges = 1;
756         }
757         else if (xleft == tile2.x)
758         {
759             result.tile1[0] = 1;
760             result.tile2[0] = 3;
761             result.numEdges = 1;
762         }
763         else if (y == tile2.y + 1)
764         {
765             result.tile1[0] = 2;
766             result.tile2[0] = 0;
767             result.numEdges = 1;
768         }
769         else if (ybottom == tile2.y)
770         {
771             result.tile1[0] = 0;
772             result.tile2[0] = 2;
773             result.numEdges = 1;
774         }
775     }
776     return result;
777 }
778 
779 void transferEdges(
780     GeoPatch* toPatch, const Matrixd& toMat, const KeyIndex& toIdx,
781     GeoPatch* fromPatch, const Matrixd& fromMat, const KeyIndex& fromIdx,
782     const ShareResult& shared);
783 
784 void safeCopy(Vec3f& dest, const Vec3f& src, const Matrixd& mat)
785 {
786     Vec3f tmp = src * mat;
787     if ((tmp - dest).length2() > 100000000)
788         OE_WARN << "whoops!\n";
789     dest = tmp;
790 }
791 
792 class TileUpdater : public NodeVisitor
793 {
794 public:
795     TileUpdater(GeoPatch* gpatch)
796         : NodeVisitor(NodeVisitor::TRAVERSE_ALL_CHILDREN), _gpatch(gpatch)
797     {
798         const MatrixTransform* trans
799             = static_cast<const MatrixTransform*>(_gpatch->getParent(0));
800         _tileMat = trans->getMatrix();
801         const PatchGroup* pg
802             = static_cast<const PatchGroup*>(trans->getParent(0));
803         const PatchOptions* popt = pg->getOptions();
804         _tileIndex = popt->getTileKey();
805     }
806 
807     void apply(PagedLOD& node)
808     {
809         PatchGroup* pgrp = dynamic_cast<PatchGroup*>(&node);
810         if (!pgrp)
811             return;
812         const PatchOptions* popt = pgrp->getOptions();
813         if (!popt)
814             return;
815         KeyIndex idx = popt->getTileKey();
816         if (idx == _tileIndex)
817             return;
818 
819         if (containsTile(idx, _tileIndex) || isNeighborTile(idx, _tileIndex))
820             copyTileEdges(pgrp, popt);
821         else if (adjoinsTile(idx, _tileIndex))
822             copyCorner(pgrp, popt);
823         else
824             return;
825         if (node.getNumChildren() > 1)
826             traverse(*node.getChild(1));
827     }
828 protected:
829     void copyTileEdges(PatchGroup* node, const PatchOptions* gopt)
830     {
831         // The tile to update
832         MatrixTransform* trans
833             = static_cast<MatrixTransform*>(node->getChild(0));
834         GeoPatch* tpatch = static_cast<GeoPatch*>(trans->getChild(0));
835         KeyIndex idx(gopt->getTileKey());
836         ShareResult shared = tilesShareEdges(idx, _tileIndex);
837         if (shared.numEdges != 0)
838         {
839             transferEdges(tpatch, trans->getMatrix(), idx,
840                           _gpatch, _tileMat, _tileIndex, shared);
841             tpatch->dirtyVertexData();
842         }
843     }
844     void copyCorner(PatchGroup* node, const PatchOptions* gopt)
845     {
846         // The tile to update
847         MatrixTransform* trans
848             = static_cast<MatrixTransform*>(node->getChild(0));
849         Matrixd toMat = trans->getMatrix();
850         Matrixd transferMat =  _tileMat * Matrixd::inverse(toMat);
851         GeoPatch* tpatch = static_cast<GeoPatch*>(trans->getChild(0));
852         KeyIndex tidx(gopt->getTileKey());
853         Geographic* gset = _gpatch->getGeographic();
854         int patchDim = gset->getResolution() + 1;
855         Vec3Array* verts = static_cast<Vec3Array*>(_gpatch->getData()
856                                                    ->vertexData.array.get());
857         PatchArray varray(*verts, patchDim);
858         Vec3Array* tverts = static_cast<Vec3Array*>(tpatch->getData()
859                                                     ->vertexData.array.get());
860         PatchArray tarray(*tverts, patchDim);
861         int lodDiff = _tileIndex.lod - tidx.lod;
862         int lodMult = 1 << lodDiff;
863         unsigned tx = tidx.x << lodDiff;
864         unsigned ty = tidx.y << lodDiff;
865 
866         if (_tileIndex.x < tx)
867         {
868             if (_tileIndex.y == ty + lodMult)
869                 //tarray[0][0] = varray[patchDim - 1][patchDim - 1] * transferMat;
870                 safeCopy(tarray[0][0], varray[patchDim - 1][patchDim - 1], transferMat);
871             else
872                 //tarray[patchDim - 1][0] = varray[0][patchDim - 1] * transferMat;
873                 safeCopy(tarray[patchDim - 1][0], varray[0][patchDim - 1], transferMat);
874         }
875         else
876         {
877             if (_tileIndex.y == ty + lodMult)
878                 // tarray[0][patchDim - 1] = varray[patchDim - 1][0] * transferMat;
879                 safeCopy(tarray[0][patchDim - 1], varray[patchDim - 1][0], transferMat);
880             else
881                 // tarray[patchDim - 1][patchDim - 1] = varray[0][0] * transferMat;
882                 safeCopy(tarray[patchDim - 1][patchDim - 1], varray[0][0], transferMat);
883         }
884         tpatch->dirtyVertexData();
885     }
886     GeoPatch* _gpatch;
887     KeyIndex _tileIndex;
888     Matrixd _tileMat;
889 };
890 
891 void transferEdges(
892     GeoPatch* toPatch, const Matrixd& toMat, const KeyIndex& toIdx,
893     GeoPatch* fromPatch, const Matrixd& fromMat, const KeyIndex& fromIdx,
894     const ShareResult& shared)
895 {
896     int resolution = toPatch->getPatchSet()->getResolution();
897     int patchDim = resolution + 1;
898     int lodDiff = fromIdx.lod - toIdx.lod;
899     int detailMult = 1 << lodDiff;
900     Matrixd transferMat = fromMat * Matrixd::inverse(toMat);
901     for (int i = 0; i < shared.numEdges; ++i)
902     {
903         EdgeRef toEdge = makeEdgeRef(toPatch, shared.tile1[i], 1);
904         EdgeRef fromEdge = makeEdgeRef(fromPatch, shared.tile2[i], detailMult);
905         int toStart;
906         if (shared.tile1[i] == 0 || shared.tile1[i] == 2)
907             toStart = (fromIdx.x - (toIdx.x * detailMult)) * resolution / detailMult;
908         else
909             toStart = (detailMult - 1 - (fromIdx.y - (toIdx.y * detailMult)))
910                 * resolution / detailMult;
911         for (int jt = toStart, jf = 0; jf < fromEdge.shape(); ++jt, ++jf)
912         {
913 #if 0
914             Vec3d vtx = Vec3d(fromEdge[jf]);
915             vtx = vtx * transferMat;
916             toEdge[jt] = Vec3f(vtx);
917 #endif
918             safeCopy(toEdge[jt], fromEdge[jf], transferMat);
919         }
920     }
921 }
922 
923 void GeoPatchUpdateCallback::operator()(Node* node, NodeVisitor* nv)
924 {
925     GeoPatch* patch = dynamic_cast<GeoPatch*>(node);
926     if (!patch)
927         return;
928     if (_hfRequest.valid() && _hfRequest->isCompleted())
929     {
930         Vec3Array* verts = dynamic_cast<Vec3Array*>(_hfRequest->getResult());
931         Vec3Array* norms = _hfRequest->_normalResult.get();
932         if (verts && norms)
933         {
934             Vec3Array* patchVerts = static_cast<Vec3Array*>(
935                 patch->getData()->vertexData.array.get());
936             Vec3Array* patchNorms = static_cast<Vec3Array*>(
937                 patch->getData()->normalData.array.get());
938             copy(verts->begin(), verts->end(), patchVerts->begin());
939             patchVerts->dirty();
940             copy(norms->begin(), norms->end(), patchNorms->begin());
941             patchNorms->dirty();
942         }
943         _hfRequest = 0;
944         PatchGroup* faceRoot = findFaceRoot(patch, nv->getNodePath());
945         if (faceRoot)
946         {
947             TileUpdater tileUpdater(patch);
948             faceRoot->accept(tileUpdater);
949         }
950     }
951     if (_imageRequest.valid() && _imageRequest->isCompleted())
952     {
953         Image* image = dynamic_cast<Image*>(_imageRequest->getResult());
954         if (image)
955         {
956             Texture2D* tex = new Texture2D();
957             tex->setImage(image);
958             tex->setWrap(Texture::WRAP_S, Texture::CLAMP_TO_EDGE);
959             tex->setWrap(Texture::WRAP_T, Texture::CLAMP_TO_EDGE);
960             tex->setFilter(Texture::MIN_FILTER,
961                            Texture::LINEAR_MIPMAP_LINEAR);
962             tex->setFilter(Texture::MAG_FILTER, Texture::LINEAR);
963             StateSet* ss = patch->getOrCreateStateSet();
964             ss->setTextureAttributeAndModes(0, tex, StateAttribute::ON);
965         }
966         _imageRequest = 0;
967     }
968     if (!_hfRequest.valid() && !_imageRequest.valid())
969         node->setUpdateCallback(0);
970 }
971 }
972