瓦片纹理申请和重投影
本文跟踪了如何从目标坐标系计算的瓦片，去查找源坐标系的瓦片进行加载。
主要思路：
- 基于目标坐标系的level级别，计算当前距离需要显示的瓦片TileKey
- 将此瓦片的范围转换到源数据所在坐标系（服务器获取的瓦片所对应的坐标系）
- 计算相交关系，找到目标坐标系的范围在源坐标系所覆盖的瓦片（可能对应多张）
- 申请源数据瓦片，将瓦片缓存，等待读取完毕，执行瓦片拼接为一张大图
- 将目标坐标系的范围变换到源坐标系后的范围与源坐标系的范围，按照距离（米）进行映射
- 重新采样，获取一张新的image，作为目标瓦片显示的新瓦片。
坐标系设置
osgEarth：内部使用两类坐标系：数据源坐标系：如瓦片所使用的坐标系；其二：目标坐标系，即显示坐标系
TerrainLayer----》创建TileSource，关联数据源坐标系
TerrainNode 绑定目标坐标系，同时将更新Terrain下的所有Layer，设置targetProfile
// Tell all the layers about the profile
for (TerrainLayerVector::iterator i = layers.begin(); i != layers.end(); ++i)
{
if (i->get()->getEnabled())
{
i->get()->setTargetProfileHint(_profile.get());
}
}
瓦片申请：
-
从RexTerrainEngineNode::dirtyTerrain开始，处理瓦片
-
使用目标标系，根据初始瓦片数目，计算第一级对应的x、y方向的瓦片数目
getAllKeysAtLOD -
根据瓦片，创建tileNode:
tileNode->create( keys[i], 0L, _engineContext.get() ); -
同步加载瓦片：
tileNode->loadSync();
-
RexTerrainEngineNode::dirtyTerrain()
{
................
std::vector<TileKey> keys;
getMap()->getProfile()->getAllKeysAtLOD( *_terrainOptions.firstLOD(), keys );
for( unsigned i=0; i<keys.size(); ++i )
{
TileNode* tileNode = new TileNode();
..................
// Next, build the surface geometry for the node.
tileNode->create( keys[i], 0L, _engineContext.get() );
// Add it to the scene graph
_terrain->addChild( tileNode );
// And load the tile's data synchronously (only for root tiles)
tileNode->loadSync();
}
...
}
TileNode:LoadSync调用LoadTileData::invoke,创建tileModel
// invoke runs in the background pager thread.
void
LoadTileData::invoke(ProgressCallback* progress)
{
// Assemble all the components necessary to display this tile
_dataModel = engine->createTileModel(
map.get(),
tilenode->getKey(),
_filter,
_enableCancel? progress : 0L);
}
createTileModel内部调用TerrainTileModelFactory::createTileModel来实现创建
TerrainTileModel*
TerrainTileModelFactory::createTileModel(const Map* map,
const TileKey& key,
const CreateTileModelFilter& filter,
const TerrainEngineRequirements* requirements,
ProgressCallback* progress)
{
// Make a new model:
osg::ref_ptr<TerrainTileModel> model = new TerrainTileModel(
key,
map->getDataModelRevision() );
// assemble all the components:
addColorLayers(model.get(), map, requirements, key, filter, progress);
addPatchLayers(model.get(), map, key, filter, progress);
if ( requirements == 0L || requirements->elevationTexturesRequired() )
{
unsigned border = requirements->elevationBorderRequired() ? 1u : 0u;
addElevation( model.get(), map, key, filter, border, progress );
}
// done.
return model.release();
}
在由addColorLayers来加载图层，如果imageLayer在逻辑范围内，执行创建GeoImage对象
void
TerrainTileModelFactory::addColorLayers(TerrainTileModel* model,
const Map* map,
const TerrainEngineRequirements* reqs,
const TileKey& key,
const CreateTileModelFilter& filter,
ProgressCallback* progress)
{
LayerVector layers;
map->getLayers(layers);
for (LayerVector::const_iterator i = layers.begin(); i != layers.end(); ++i)
{
Layer* layer = i->get();
if (layer->getRenderType() != layer->RENDERTYPE_TERRAIN_SURFACE)
continue;
if (!layer->getEnabled())
continue;
if (!filter.accept(layer))
continue;
ImageLayer* imageLayer = dynamic_cast<ImageLayer*>(layer);
if (imageLayer)
{
osg::Texture* tex = 0L;
osg::Matrixf textureMatrix;
if (imageLayer->isKeyInLegalRange(key) && imageLayer->mayHaveData(key))
{
if (imageLayer->useCreateTexture())
{
tex = imageLayer->createTexture( key, progress, textureMatrix );
}
else
{
GeoImage geoImage = imageLayer->createImage( key, progress );
if ( geoImage.valid() )
{
if ( imageLayer->isCoverage() )
tex = createCoverageTexture(geoImage.getImage(), imageLayer);
else
tex = createImageTexture(geoImage.getImage(), imageLayer);
}
}
}
// if this is the first LOD, and the engine requires that the first LOD
// be populated, make an empty texture if we didn't get one.
if (tex == 0L &&
_options.firstLOD() == key.getLOD() &&
reqs && reqs->fullDataAtFirstLodRequired())
{
tex = _emptyTexture.get();
}
if (tex)
{
tex->setName(model->getKey().str());
TerrainTileImageLayerModel* layerModel = new TerrainTileImageLayerModel();
layerModel->setImageLayer(imageLayer);
layerModel->setTexture(tex);
layerModel->setMatrix(new osg::RefMatrixf(textureMatrix));
model->colorLayers().push_back(layerModel);
if (imageLayer->isShared())
{
model->sharedLayers().push_back(layerModel);
}
if (imageLayer->isDynamic())
{
model->setRequiresUpdateTraverse(true);
}
}
}
else // non-image kind of TILE layer:
{
TerrainTileColorLayerModel* colorModel = new TerrainTileColorLayerModel();
colorModel->setLayer(layer);
model->colorLayers().push_back(colorModel);
}
}
}
范围判断：isKeyInLegalRange
TerrainLayer中保存的Profile是源数据的坐标系（瓦片图片所使用的坐标系）
bool TerrainLayer::isKeyInLegalRange(const TileKey& key) const
{
// We must use the equivalent lod b/c the input key can be in any profile.
// 获取目标坐标系和瓦片坐标系在接近的一个lod级别
// getProfile(): 源坐标系
// key.getProfile()： 目标坐标系
unsigned localLOD = getProfile() ?
getProfile()->getEquivalentLOD(key.getProfile(), key.getLOD()) :
key.getLOD();
// 级别范围筛选
// First check the key against the min/max level limits, it they are set.
if ((options().maxLevel().isSet() && localLOD > options().maxLevel().value()) ||
(options().minLevel().isSet() && localLOD < options().minLevel().value()))
{
return false;
}
// Next check the maxDataLevel if that is set.
if (options().maxDataLevel().isSet() && localLOD > options().maxDataLevel().get())
{
return false;
}
// 分辨率大小筛选
// Next, check against resolution limits (based on the source tile size).
if (options().minResolution().isSet() || options().maxResolution().isSet())
{
const Profile* profile = getProfile();
if ( profile )
{
// calculate the resolution in the layer's profile, which can
// be different that the key's profile.
// 获取目标坐标系下的宽度，和瓦片大小，计算尺寸
double resKey = key.getExtent().width() / (double)getTileSize();
double resLayer = key.getProfile()->getSRS()->transformUnits(resKey, profile->getSRS());
if (options().maxResolution().isSet() &&
options().maxResolution().value() > resLayer)
{
return false;
}
if (options().minResolution().isSet() &&
options().minResolution().value() < resLayer)
{
return false;
}
}
}
return true;
}
获取最接近的lod：
unsigned
Profile::getEquivalentLOD( const Profile* rhsProfile, unsigned rhsLOD ) const
{ // this 源坐标系
//If the profiles are equivalent, just use the incoming lod
if (rhsProfile->isHorizEquivalentTo( this ) ) // 目标坐标系
return rhsLOD;
// Special check for geodetic to mercator or vise versa, they should match up in LOD.
if ((rhsProfile->isEquivalentTo(Registry::instance()->getSphericalMercatorProfile()) && isEquivalentTo(Registry::instance()->getGlobalGeodeticProfile())) ||
(rhsProfile->isEquivalentTo(Registry::instance()->getGlobalGeodeticProfile()) && isEquivalentTo(Registry::instance()->getSphericalMercatorProfile())))
{
return rhsLOD;
}
double rhsWidth, rhsHeight; // 获取目标坐标系当前lod，单张瓦片的尺寸范围，基于范围计算得来
rhsProfile->getTileDimensions( rhsLOD, rhsWidth, rhsHeight );
// safety catch
if ( osg::equivalent(rhsWidth, 0.0) || osg::equivalent(rhsHeight, 0.0) )
{
OE_WARN << LC << "getEquivalentLOD: zero dimension" << std::endl;
return rhsLOD;
}
// 将目标坐标系下的尺寸单位，变换到源坐标系下的尺寸单位，对于gcgs2000适用吗？
const SpatialReference* rhsSRS = rhsProfile->getSRS();
double rhsTargetHeight = rhsSRS->transformUnits( rhsHeight, getSRS() );
int currLOD = 0;
int destLOD = currLOD;
double delta = DBL_MAX;
// Find the LOD that most closely matches the resolution of the incoming key.
// We use the closest (under or over) so that you can match back and forth between profiles and be sure to get the same results each time.
// 找到源坐标系瓦片同目标坐标系瓦片最匹配的lod值
while( true )
{
double prevDelta = delta;
double w, h;
getTileDimensions(currLOD, w, h); // 获取源坐标系，某一级别下的单张瓦片尺寸
delta = osg::absolute( h - rhsTargetHeight );
if (delta < prevDelta)
{
// We're getting closer so keep going
destLOD = currLOD;
}
else
{
// We are further away from the previous lod so stop.
break;
}
currLOD++;
}
return destLOD;
}
判断当前key是否在有效范围内
bool
TerrainLayer::isKeyInLegalRange(const TileKey& key) const
{
if ( !key.valid() )
{
return false;
}
// We must use the equivalent lod b/c the input key can be in any profile.
unsigned localLOD = getProfile() ?
getProfile()->getEquivalentLOD(key.getProfile(), key.getLOD()) :
key.getLOD();
// First check the key against the min/max level limits, it they are set.
if ((options().maxLevel().isSet() && localLOD > options().maxLevel().value()) ||
(options().minLevel().isSet() && localLOD < options().minLevel().value()))
{
return false;
}
// Next check the maxDataLevel if that is set.
if (options().maxDataLevel().isSet() && localLOD > options().maxDataLevel().get())
{
return false;
}
// Next, check against resolution limits (based on the source tile size).
if (options().minResolution().isSet() || options().maxResolution().isSet())
{
const Profile* profile = getProfile();
if ( profile )
{
// calculate the resolution in the layer's profile, which can
// be different that the key's profile.
double resKey = key.getExtent().width() / (double)getTileSize();
double resLayer = key.getProfile()->getSRS()->transformUnits(resKey, profile->getSRS());
if (options().maxResolution().isSet() &&
options().maxResolution().value() > resLayer)
{
return false;
}
if (options().minResolution().isSet() &&
options().minResolution().value() < resLayer)
{
return false;
}
}
}
return true;
}
GeoImage 创建
GeoImage geoImage = imageLayer->createImage( key, progress );
GeoImage ImageLayer::createImage(const TileKey& key,
ProgressCallback* progress)
{
ScopedMetric m("ImageLayer::createImage", 2,
"key", key.str().c_str(),
"name", getName().c_str());
if (getStatus().isError())
{
return GeoImage::INVALID;
}
return createImageInKeyProfile( key, progress );
}
根据目标坐标系和源坐标系是否相同，执行两个分支
if (key.getProfile()->isHorizEquivalentTo(getProfile()))
{
result = createImageImplementation(key, progress);
}
else
{
// If the profiles are different, use a compositing method to assemble the tile.
result = assembleImage( key, progress );
}
相同坐标系分支：
GeoImage
ImageLayer::createImageImplementation(const TileKey& key, ProgressCallback* progress)
{
return createImageFromTileSource(key, progress);
}
GeoImage ImageLayer::createImageFromTileSource(const TileKey& key, ProgressCallback* progress)
{
TileSource* source = getTileSource();
if ( !source )
return GeoImage::INVALID;
// If the profiles are different, use a compositing method to assemble the tile.
if ( !key.getProfile()->isHorizEquivalentTo( getProfile() ) )
{
return assembleImage( key, progress );
}
// Good to go, ask the tile source for an image:
osg::ref_ptr<TileSource::ImageOperation> op = getOrCreatePreCacheOp();
// Fail is the image is blacklisted.
if ( source->getBlacklist()->contains(key) )
{
OE_DEBUG << LC << "createImageFromTileSource: blacklisted(" << key.str() << ")" << std::endl;
return GeoImage::INVALID;
}
if (!mayHaveData(key))
{
OE_DEBUG << LC << "createImageFromTileSource: mayHaveData(" << key.str() << ") == false" << std::endl;
return GeoImage::INVALID;
}
// create an image from the tile source.
osg::ref_ptr<osg::Image> result = source->createImage( key, op.get(), progress );
// Process images with full alpha to properly support MP blending.
if (result.valid() &&
options().featherPixels() == true)
{
ImageUtils::featherAlphaRegions( result.get() );
}
// If image creation failed (but was not intentionally canceled and
// didn't time out or end for any other recoverable reason), then
// blacklist this tile for future requests.
if (result == 0L)
{
if ( progress == 0L || !progress->isCanceled() )
{
source->getBlacklist()->add( key );
}
}
if (progress && progress->isCanceled())
{
return GeoImage::INVALID;
}
return GeoImage(result.get(), key.getExtent());
}
由tileSource直接创建geoImage：,如：tms
osg::Image* TMSTileSource::createImage(const TileKey& key, ProgressCallback* progress)
{
if (_tileMap.valid() && key.getLevelOfDetail() <= _tileMap->getMaxLevel() )
{
std::string image_url = _tileMap->getURL( key, _invertY );
osg::ref_ptr<osg::Image> image;
if (!image_url.empty())
{
URI uri(image_url, _options.url()->context());
image = uri.readImage( _dbOptions.get(), progress ).getImage();
}
if (!image.valid())
{
if (image_url.empty() || !_tileMap->intersectsKey(key))
{
//We couldn't read the image from the URL or the cache, so check to see if the given key is less than the max level
//of the tilemap and create a transparent image.
if (key.getLevelOfDetail() <= _tileMap->getMaxLevel())
{
OE_DEBUG << LC << "Returning empty image " << std::endl;
return ImageUtils::createEmptyImage();
}
}
}
if (image.valid() && _options.coverage() == true)
{
image->setInternalTextureFormat(GL_R16F);
ImageUtils::markAsUnNormalized(image.get(), true);
}
return image.release();
}
return 0;
}
不同坐标系分支，执行组装图片（重投影)：
GeoImage ImageLayer::assembleImage(const TileKey& key, ProgressCallback* progress)
{
// If we got here, asset that there's a non-null layer profile.
if (!getProfile()) // 源坐标系
{
setStatus(Status::Error(Status::AssertionFailure, "assembleImage with undefined profile"));
return GeoImage::INVALID;
}
GeoImage mosaicedImage, result;
// Scale the extent if necessary to apply an "edge buffer"
GeoExtent ext = key.getExtent(); // 获取目标坐标系范围
if ( options().edgeBufferRatio().isSet() )
{
double ratio = options().edgeBufferRatio().get();
ext.scale(ratio, ratio);
}
// 获取目标瓦片在源瓦片中的相交部分
// Get a set of layer tiles that intersect the requested extent.
std::vector<TileKey> intersectingKeys; //
getProfile()->getIntersectingTiles( key, intersectingKeys );
if ( intersectingKeys.size() > 0 )
{
double dst_minx, dst_miny, dst_maxx, dst_maxy;
key.getExtent().getBounds(dst_minx, dst_miny, dst_maxx, dst_maxy);
// if we find at least one "real" tile in the mosaic, then the whole result tile is
// "real" (i.e. not a fallback tile)
bool retry = false;
ImageMosaic mosaic;
// keep track of failed tiles.
std::vector<TileKey> failedKeys;
// 遍历所有获取的瓦片，创建瓦片图像
for( std::vector<TileKey>::iterator k = intersectingKeys.begin(); k != intersectingKeys.end(); ++k )
{
GeoImage image = createImageImplementation( *k, progress );
if ( image.valid() )
{
if ( !isCoverage() )
{
ImageUtils::fixInternalFormat(image.getImage());
// Make sure all images in mosaic are based on "RGBA - unsigned byte" pixels.
// This is not the smarter choice (in some case RGB would be sufficient) but
// it ensure consistency between all images / layers.
//
// The main drawback is probably the CPU memory foot-print which would be reduced by allocating RGB instead of RGBA images.
// On GPU side, this should not change anything because of data alignements : often RGB and RGBA textures have the same memory footprint
if ( (image.getImage()->getDataType() != GL_UNSIGNED_BYTE)
|| (image.getImage()->getPixelFormat() != GL_RGBA) )
{
osg::ref_ptr<osg::Image> convertedImg = ImageUtils::convertToRGBA8(image.getImage());
if (convertedImg.valid())
{
image = GeoImage(convertedImg.get(), image.getExtent());
}
}
}
mosaic.getImages().push_back( TileImage(image.getImage(), *k) );
}
else
{
// the tile source did not return a tile, so make a note of it.
failedKeys.push_back( *k );
if (progress && progress->isCanceled())
{
retry = true;
break;
}
}
}
// Fail is: a) we got no data and the LOD is greater than zero; or
// b) the operation was canceled mid-stream.
if ( (mosaic.getImages().empty() && key.getLOD() > 0) || retry)
{
// if we didn't get any data at LOD>0, fail.
OE_DEBUG << LC << "Couldn't create image for ImageMosaic " << std::endl;
return GeoImage::INVALID;
}
// We got at least one good tile, OR we got nothing but since the LOD==0 we have to
// fall back on a lower resolution.
// So now we go through the failed keys and try to fall back on lower resolution data
// to fill in the gaps. The entire mosaic must be populated or this qualifies as a bad tile.
for(std::vector<TileKey>::iterator k = failedKeys.begin(); k != failedKeys.end(); ++k)
{
GeoImage image;
for(TileKey parentKey = k->createParentKey();
parentKey.valid() && !image.valid();
parentKey = parentKey.createParentKey())
{
image = createImageImplementation( parentKey, progress );
if ( image.valid() )
{
GeoImage cropped;
if ( !isCoverage() )
{
ImageUtils::fixInternalFormat(image.getImage());
if ( (image.getImage()->getDataType() != GL_UNSIGNED_BYTE)
|| (image.getImage()->getPixelFormat() != GL_RGBA) )
{
osg::ref_ptr<osg::Image> convertedImg = ImageUtils::convertToRGBA8(image.getImage());
if (convertedImg.valid())
{
image = GeoImage(convertedImg.get(), image.getExtent());
}
}
cropped = image.crop( k->getExtent(), false, image.getImage()->s(), image.getImage()->t() );
}
else
{
// TODO: may not work.... test; tilekey extent will <> cropped extent
cropped = image.crop( k->getExtent(), true, image.getImage()->s(), image.getImage()->t(), false );
}
// and queue it.
mosaic.getImages().push_back( TileImage(cropped.getImage(), *k) );
}
}
if ( !image.valid() )
{
// a tile completely failed, even with fallback. Eject.
OE_DEBUG << LC << "Couldn't fallback on tiles for ImageMosaic" << std::endl;
// let it go. The empty areas will be filled with alpha by ImageMosaic.
}
}
// all set. Mosaic all the images together.
double rxmin, rymin, rxmax, rymax;
mosaic.getExtents( rxmin, rymin, rxmax, rymax );
mosaicedImage = GeoImage(
mosaic.createImage(),
GeoExtent( getProfile()->getSRS(), rxmin, rymin, rxmax, rymax ) );
}
else
{
OE_DEBUG << LC << "assembleImage: no intersections (" << key.str() << ")" << std::endl;
}
// Final step: transform the mosaic into the requesting key's extent.
if ( mosaicedImage.valid() )
{
// GeoImage::reproject() will automatically crop the image to the correct extents.
// so there is no need to crop after reprojection. Also note that if the SRS's are the
// same (even though extents are different), then this operation is technically not a
// reprojection but merely a resampling.
result = mosaicedImage.reproject(
key.getProfile()->getSRS(),
&key.getExtent(),
getTileSize(), getTileSize(),
options().driver()->bilinearReprojection().get());
}
// Process images with full alpha to properly support MP blending.
if (result.valid() &&
options().featherPixels() == true &&
isCoverage() == false)
{
ImageUtils::featherAlphaRegions( result.getImage() );
}
if (progress && progress->isCanceled())
{
return GeoImage::INVALID;
}
return result;
}
profile求交：
void Profile::getIntersectingTiles(const TileKey& key, std::vector<TileKey>& out_intersectingKeys) const
{
OE_DEBUG << "GET ISECTING TILES for key " << key.str() << " -----------------" << std::endl;
//If the profiles are exactly equal, just add the given tile key.
if ( isHorizEquivalentTo( key.getProfile() ) )
{
//Clear the incoming list
out_intersectingKeys.clear();
out_intersectingKeys.push_back(key);
}
else
{
// figure out which LOD in the local profile is a best match for the LOD
// in the source LOD in terms of resolution.
unsigned localLOD = getEquivalentLOD(key.getProfile(), key.getLOD());
getIntersectingTiles(key.getExtent(), localLOD, out_intersectingKeys);
OE_DEBUG << LC << "GIT, key="<< key.str() << ", localLOD=" << localLOD
<< ", resulted in " << out_intersectingKeys.size() << " tiles" << std::endl;
}
}
求交：判断是否跨越国际日期变更线
void Profile::getIntersectingTiles(const GeoExtent& extent, unsigned localLOD, std::vector<TileKey>& out_intersectingKeys) const
{// this: 源坐标系
GeoExtent ext = extent;
// reproject into the profile's SRS if necessary:
if ( !getSRS()->isHorizEquivalentTo( extent.getSRS() ) )
{
// localize the extents and clamp them to legal values
// 将目标坐标范围和源坐标系范围进行求教
ext = clampAndTransformExtent( extent );
if ( !ext.isValid() )
return;
}
// 是否跨越国际日期变更线，是：添加连个范围进行比较
if ( ext.crossesAntimeridian() )
{
GeoExtent first, second;
if (ext.splitAcrossAntimeridian( first, second ))
{
addIntersectingTiles( first, localLOD, out_intersectingKeys );
addIntersectingTiles( second, localLOD, out_intersectingKeys );
}
}
else
{
addIntersectingTiles( ext, localLOD, out_intersectingKeys );
}
}
截取范围：
- 将目标坐标系范围转换到源坐标系范围
- 转换成功：用两个范围进行求交，返回求交结果的交集
- 失败：
- 将目标范围，转换到源坐标的经纬度坐标系，尝试在经纬度范围求交
- 将转换后的坐标系和实际的源坐标系的地理范围进行截取
- 再将截取后的范围，重新投影到源坐标系的范围，返回结果
GeoExtent Profile::clampAndTransformExtent(const GeoExtent& input, bool* out_clamped) const
{// this：源坐标系
// initialize the output flag
if ( out_clamped )
*out_clamped = false;
// begin by transforming the input extent to this profile's SRS.
// 将目标坐标系范围转换到源坐标系范围
GeoExtent inputInMySRS = input.transform(getSRS());
if (inputInMySRS.isValid()) // 转换成功，
{
// Compute the intersection of the two:
GeoExtent intersection = inputInMySRS.intersectionSameSRS(getExtent());
// expose whether clamping took place:
if (out_clamped != 0L)
{
*out_clamped = intersection != getExtent();
}
return intersection;
}
else
{
// The extent transformation failed, probably due to an out-of-bounds condition.
// Go to Plan B: attempt the operation in lat/long
// 转换失败，可能是范围越界导致，执行方案b，尝试在经纬度范围求交
const SpatialReference* geo_srs = getSRS()->getGeographicSRS();
// get the input in lat/long:
// 将目标坐标系转换为源坐标系所对应的地理坐标系（经纬度）
GeoExtent gcs_input =
input.getSRS()->isGeographic()?
input :
input.transform( geo_srs );
// bail out on a bad transform:
if ( !gcs_input.isValid() )
return GeoExtent::INVALID;
// bail out if the extent's do not intersect at all:
if ( !gcs_input.intersects(_latlong_extent, false) )
return GeoExtent::INVALID;
// clamp it to the profile's extents:
//将转换后的坐标系和实际的源坐标系的地理范围进行截取
GeoExtent clamped_gcs_input = GeoExtent(
gcs_input.getSRS(),
osg::clampBetween( gcs_input.xMin(), _latlong_extent.xMin(), _latlong_extent.xMax() ),
osg::clampBetween( gcs_input.yMin(), _latlong_extent.yMin(), _latlong_extent.yMax() ),
osg::clampBetween( gcs_input.xMax(), _latlong_extent.xMin(), _latlong_extent.xMax() ),
osg::clampBetween( gcs_input.yMax(), _latlong_extent.yMin(), _latlong_extent.yMax() ) );
if ( out_clamped )
*out_clamped = (clamped_gcs_input != gcs_input);
// 再将截取后的范围，重新投影到源坐标系的范围
// finally, transform the clamped extent into this profile's SRS and return it.
GeoExtent result =
clamped_gcs_input.getSRS()->isEquivalentTo( this->getSRS() )?
clamped_gcs_input :
clamped_gcs_input.transform( this->getSRS() );
if (result.isValid())
{
OE_DEBUG << LC << "clamp&xform: input=" << input.toString() << ", output=" << result.toString() << std::endl;
}
return result;
}
}
使用源坐标范围求交瓦片：addIntersectingTiles
void Profile::getIntersectingTiles(const TileKey& key, std::vector<TileKey>& out_intersectingKeys) const
{
OE_DEBUG << "GET ISECTING TILES for key " << key.str() << " -----------------" << std::endl;
//If the profiles are exactly equal, just add the given tile key.
if ( isHorizEquivalentTo( key.getProfile() ) )
{
//Clear the incoming list
out_intersectingKeys.clear();
out_intersectingKeys.push_back(key);
}
else
{
// figure out which LOD in the local profile is a best match for the LOD
// in the source LOD in terms of resolution.
unsigned localLOD = getEquivalentLOD(key.getProfile(), key.getLOD());
getIntersectingTiles(key.getExtent(), localLOD, out_intersectingKeys);
OE_DEBUG << LC << "GIT, key="<< key.str() << ", localLOD=" << localLOD
<< ", resulted in " << out_intersectingKeys.size() << " tiles" << std::endl;
}
}
void
Profile::addIntersectingTiles(const GeoExtent& key_ext, unsigned localLOD, std::vector<TileKey>& out_intersectingKeys) const
{
// assume a non-crossing extent here.
if ( key_ext.crossesAntimeridian() )
{
OE_WARN << "Profile::addIntersectingTiles cannot process date-line cross" << std::endl;
return;
}
int tileMinX, tileMaxX;
int tileMinY, tileMaxY;
double destTileWidth, destTileHeight;
getTileDimensions(localLOD, destTileWidth, destTileHeight);
//OE_DEBUG << std::fixed << " Source Tile: " << key.getLevelOfDetail() << " (" << keyWidth << ", " << keyHeight << ")" << std::endl;
//OE_DEBUG << std::fixed << " Dest Size: " << destLOD << " (" << destTileWidth << ", " << destTileHeight << ")" << std::endl;
double east = key_ext.xMax() - _extent.xMin();
bool xMaxOnTileBoundary = fmod(east, destTileWidth) == 0.0;
double south = _extent.yMax() - key_ext.yMin();
bool yMaxOnTileBoundary = fmod(south, destTileHeight) == 0.0;
tileMinX = (int)((key_ext.xMin() - _extent.xMin()) / destTileWidth);
// 在边界上，则-1，否则为0
tileMaxX = (int)(east / destTileWidth) - (xMaxOnTileBoundary ? 1 : 0);
tileMinY = (int)((_extent.yMax() - key_ext.yMax()) / destTileHeight);
tileMaxY = (int)(south / destTileHeight) - (yMaxOnTileBoundary ? 1 : 0);
unsigned int numWide, numHigh;
getNumTiles(localLOD, numWide, numHigh);
// bail out if the tiles are out of bounds.
if ( tileMinX >= (int)numWide || tileMinY >= (int)numHigh ||
tileMaxX < 0 || tileMaxY < 0 )
{
return;
}
tileMinX = osg::clampBetween(tileMinX, 0, (int)numWide-1);
tileMaxX = osg::clampBetween(tileMaxX, 0, (int)numWide-1);
tileMinY = osg::clampBetween(tileMinY, 0, (int)numHigh-1);
tileMaxY = osg::clampBetween(tileMaxY, 0, (int)numHigh-1);
OE_DEBUG << std::fixed << " Dest Tiles: " << tileMinX << "," << tileMinY << " => " << tileMaxX << "," << tileMaxY << std::endl;
// 获取源坐标系范围内的相交瓦片
for (int i = tileMinX; i <= tileMaxX; ++i)
{
for (int j = tileMinY; j <= tileMaxY; ++j)
{
//TODO: does not support multi-face destination keys.
out_intersectingKeys.push_back( TileKey(localLOD, i, j, this) );
}
}
}
返回上级函数，遍历intersectingkeys，申请源数据瓦片createImageImplementation
进入createImageFromTileSource，创建图片
// create an image from the tile source.
osg::ref_ptr<osg::Image> result = source->createImage( key, op.get(), progress );
osg::Image*TileSource::createImage(const TileKey& key,
ImageOperation* prepOp,
ProgressCallback* progress )
{
if (getStatus().isError())
return 0L;
// Try to get it from the memcache fist
if (_memCache.valid())
{
ReadResult r = _memCache->getOrCreateDefaultBin()->readImage(key.str(), 0L);
if ( r.succeeded() )
return r.releaseImage();
}
osg::ref_ptr<osg::Image> newImage = createImage(key, progress);
// Check for cancelation. The TileSource implementation should do this
// internally but we check here once last time just in case the
// implementation does not.
if (progress && progress->isCanceled())
{
return 0L;
}
// Run the pre-caching operation if there is one:
if ( prepOp )
(*prepOp)( newImage );
// Cache to the L2 cache:
if ( newImage.valid() && _memCache.valid() )
{
_memCache->getOrCreateDefaultBin()->write(key.str(), newImage.get(), 0L);
}
return newImage.release();
}
最终，返回GeoImage，包含的crs为源图层坐标系，读取的瓦片，保存在mosaic.getImages数组中
mosaic.getExtents( rxmin, rymin, rxmax, rymax ),内部遍历其所有瓦片范围，然创建一张大图，进行合并：
// all set. Mosaic all the images together.
double rxmin, rymin, rxmax, rymax;
mosaic.getExtents( rxmin, rymin, rxmax, rymax );
mosaicedImage = GeoImage(
mosaic.createImage(),
GeoExtent( getProfile()->getSRS(), rxmin, rymin, rxmax, rymax ) );
从mosaic中创建合并的图片
osg::Image*
ImageMosaic::createImage()
{
if (_images.size() == 0)
{
OE_INFO << "ImageMosaic has no images..." << std::endl;
return 0;
}
// find the first valid tile and use its size as the mosaic tile size
TileImage* tile = 0L;
for (int i = 0; i<_images.size() && !tile; ++i)
if (_images[i]._image.valid())
tile = &_images[i];
if ( !tile )
return 0L;
unsigned int tileWidth = tile->_image->s();
unsigned int tileHeight = tile->_image->t();
//OE_NOTICE << "TileDim " << tileWidth << ", " << tileHeight << std::endl;
unsigned int minTileX = tile->_tileX;
unsigned int minTileY = tile->_tileY;
unsigned int maxTileX = tile->_tileX;
unsigned int maxTileY = tile->_tileY;
//Compute the tile size.
for (TileImageList::iterator i = _images.begin(); i != _images.end(); ++i)
{
if (i->_tileX < minTileX) minTileX = i->_tileX;
if (i->_tileY < minTileY) minTileY = i->_tileY;
if (i->_tileX > maxTileX) maxTileX = i->_tileX;
if (i->_tileY > maxTileY) maxTileY = i->_tileY;
}
unsigned int tilesWide = maxTileX - minTileX + 1;
unsigned int tilesHigh = maxTileY - minTileY + 1;
unsigned int pixelsWide = tilesWide * tileWidth;
unsigned int pixelsHigh = tilesHigh * tileHeight;
unsigned int tileDepth = tile->_image->r();
osg::ref_ptr<osg::Image> image = new osg::Image;
image->allocateImage(pixelsWide, pixelsHigh, tileDepth, tile->_image->getPixelFormat(), tile->_image->getDataType());
image->setInternalTextureFormat(tile->_image->getInternalTextureFormat());
ImageUtils::markAsNormalized(image.get(), ImageUtils::isNormalized(tile->getImage()));
//Initialize the image to be completely white!
//memset(image->data(), 0xFF, image->getImageSizeInBytes());
ImageUtils::PixelWriter write(image.get());
for (unsigned t = 0; t < pixelsHigh; ++t)
for (unsigned s = 0; s < pixelsWide; ++s)
write(osg::Vec4(1,1,1,0), s, t);
//Composite the incoming images into the master image
// 将所有瓦片拷贝到同一个瓦片中
for (TileImageList::iterator i = _images.begin(); i != _images.end(); ++i)
{
//Determine the indices in the master image for this image
int dstX = (i->_tileX - minTileX) * tileWidth;
int dstY = (maxTileY - i->_tileY) * tileHeight;
osg::Image* sourceTile = i->getImage();
if ( sourceTile )
{
ImageUtils::copyAsSubImage(sourceTile, image.get(), dstX, dstY);
}
}
return image.release();
}
对瓦片进行重投影：
result = mosaicedImage.reproject(
key.getProfile()->getSRS(),
&key.getExtent(),
getTileSize(), getTileSize(),
options().driver()->bilinearReprojection().get());
重投影分为gdal投影和自定义投影
GeoImage
GeoImage::reproject(const SpatialReference* to_srs, const GeoExtent* to_extent, unsigned int width, unsigned int height, bool useBilinearInterpolation) const
{
GeoExtent destExtent;
if (to_extent)
{
destExtent = *to_extent;
}
else
{
destExtent = getExtent().transform(to_srs);
}
osg::Image* resultImage = 0L;
bool isNormalized = ImageUtils::isNormalized(getImage());
if ( getSRS()->isUserDefined() ||
to_srs->isUserDefined() ||
getSRS()->isSphericalMercator() ||
to_srs->isSphericalMercator() ||
!isNormalized )
{
// if either of the SRS is a custom projection, we have to do a manual reprojection since
// GDAL will not recognize the SRS.
resultImage = manualReproject(getImage(), getExtent(), destExtent, useBilinearInterpolation && isNormalized, width, height);
}
else
{
// otherwise use GDAL.
resultImage = reprojectImage(
getImage(),
getSRS()->getWKT(),
getExtent().xMin(), getExtent().yMin(), getExtent().xMax(), getExtent().yMax(),
to_srs->getWKT(),
destExtent.xMin(), destExtent.yMin(), destExtent.xMax(), destExtent.yMax(),
width, height, useBilinearInterpolation);
}
return GeoImage(resultImage, destExtent);
}
手动投影：
- 定义一个和目标瓦片图片width和height一致的image
- 计算目标瓦片的pixelformeter，每像素多少米（或单位）
- 计算图片占用多少像素数目
- 使用transformExtentPoints，将这个范围的点，按照像素采样（一个像素可能占了很多的单位大小），重投影到源坐标系空间
- 将投影到源坐标系的网格，从源中读取i昂二的的每个点，将他写入到目标图片的同一位置。
- 根据当前位置和边界的距离，计算比例进行重新映射，算出目标范围所对应的像素值的索引
- 基于当前px,py取临近四个方向的点，进行线性采样
- 计算后的图，作为目标瓦片的image使用。
osg::Image* manualReproject(
const osg::Image* image,
const GeoExtent& src_extent,
const GeoExtent& dest_extent,
bool interpolate,
unsigned int width = 0,
unsigned int height = 0)
{
//TODO: Compute the optimal destination size
if (width == 0 || height == 0)
{
//If no width and height are specified, just use the minimum dimension for the image
width = osg::minimum(image->s(), image->t());
height = osg::minimum(image->s(), image->t());
}
osg::Image *result = new osg::Image();
//result->allocateImage(width, height, 1, GL_RGBA, GL_UNSIGNED_BYTE);
result->allocateImage(width, height, 1, image->getPixelFormat(), image->getDataType()); //GL_UNSIGNED_BYTE);
result->setInternalTextureFormat(image->getInternalTextureFormat());
ImageUtils::markAsUnNormalized(result, ImageUtils::isUnNormalized(image));
//Initialize the image to be completely transparent/black
memset(result->data(), 0, result->getImageSizeInBytes());
//ImageUtils::PixelReader ra(result);
ImageUtils::PixelWriter writer(result);
const double dx = dest_extent.width() / (double)width;
const double dy = dest_extent.height() / (double)height;
// offset the sample points by 1/2 a pixel so we are sampling "pixel center".
// (This is especially useful in the UnifiedCubeProfile since it nullifes the chances for
// edge ambiguity.)
// 计算需要的像素数量
unsigned int numPixels = width * height;
// Start by creating a sample grid over the destination
// extent. These will be the source coordinates. Then, reproject
// the sample grid into the source coordinate system.
double *srcPointsX = new double[numPixels * 2];
double *srcPointsY = srcPointsX + numPixels;
//将目标范围，变换到源所在坐标系的范围
dest_extent.getSRS()->transformExtentPoints(
src_extent.getSRS(),
dest_extent.xMin() + .5 * dx, dest_extent.yMin() + .5 * dy,
dest_extent.xMax() - .5 * dx, dest_extent.yMax() - .5 * dy,
srcPointsX, srcPointsY, width, height);
// Next, go through the source-SRS sample grid, read the color at each point from the source image,
// and write it to the corresponding pixel in the destination image.
// 创建源坐标系下的网格，从源中读取每个点，将他写入到目标图片的同一位置。
int pixel = 0;
ImageUtils::PixelReader ia(image);
double xfac = (image->s() - 1) / src_extent.width();
double yfac = (image->t() - 1) / src_extent.height();
for (unsigned int c = 0; c < width; ++c)
{
for (unsigned int r = 0; r < height; ++r)
{
double src_x = srcPointsX[pixel];
double src_y = srcPointsY[pixel];
if ( src_x < src_extent.xMin() || src_x > src_extent.xMax() || src_y < src_extent.yMin() || src_y > src_extent.yMax() )
{
//If the sample point is outside of the bound of the source extent, increment the pixel and keep looping through.
//OE_WARN << LC << "ERROR: sample point out of bounds: " << src_x << ", " << src_y << std::endl;
pixel++;
continue;
}
// 根据当前位置和边界的距离，计算比例进行重新映射，算出目标范围所对应的像素值的索引
float px = (src_x - src_extent.xMin()) * xfac;
float py = (src_y - src_extent.yMin()) * yfac;
int px_i = osg::clampBetween( (int)osg::round(px), 0, image->s()-1 );
int py_i = osg::clampBetween( (int)osg::round(py), 0, image->t()-1 );
osg::Vec4 color(0,0,0,0);
// TODO: consider this again later. Causes blockiness.
if ( !interpolate ) //! isSrcContiguous ) // non-contiguous space- use nearest neighbot
{
// 读取像素
color = ia(px_i, py_i);
}
else // contiguous space - use bilinear sampling
{
// 基于当前px,py取临近四个方向的点，进行线性采样
int rowMin = osg::maximum((int)floor(py), 0);
int rowMax = osg::maximum(osg::minimum((int)ceil(py), (int)(image->t()-1)), 0);
int colMin = osg::maximum((int)floor(px), 0);
int colMax = osg::maximum(osg::minimum((int)ceil(px), (int)(image->s()-1)), 0);
if (rowMin > rowMax) rowMin = rowMax;
if (colMin > colMax) colMin = colMax;
// 去除四个点
osg::Vec4 urColor = ia(colMax, rowMax);
osg::Vec4 llColor = ia(colMin, rowMin);
osg::Vec4 ulColor = ia(colMin, rowMax);
osg::Vec4 lrColor = ia(colMax, rowMin);
/*Average Interpolation*/
/*double x_rem = px - (int)px;
double y_rem = py - (int)py;
double w00 = (1.0 - y_rem) * (1.0 - x_rem);
double w01 = (1.0 - y_rem) * x_rem;
double w10 = y_rem * (1.0 - x_rem);
double w11 = y_rem * x_rem;
double wsum = w00 + w01 + w10 + w11;
wsum = 1.0/wsum;
color.r() = (w00 * llColor.r() + w01 * lrColor.r() + w10 * ulColor.r() + w11 * urColor.r()) * wsum;
color.g() = (w00 * llColor.g() + w01 * lrColor.g() + w10 * ulColor.g() + w11 * urColor.g()) * wsum;
color.b() = (w00 * llColor.b() + w01 * lrColor.b() + w10 * ulColor.b() + w11 * urColor.b()) * wsum;
color.a() = (w00 * llColor.a() + w01 * lrColor.a() + w10 * ulColor.a() + w11 * urColor.a()) * wsum;*/
/*Nearest Neighbor Interpolation*/
/*if (px_i >= 0 && px_i < image->s() &&
py_i >= 0 && py_i < image->t())
{
//OE_NOTICE << "[osgEarth::GeoData] Sampling pixel " << px << "," << py << std::endl;
color = ImageUtils::getColor(image, px_i, py_i);
}
else
{
OE_NOTICE << "[osgEarth::GeoData] Pixel out of range " << px_i << "," << py_i << " image is " << image->s() << "x" << image->t() << std::endl;
}*/
/*Bilinear interpolation*/
//Check for exact value
if ((colMax == colMin) && (rowMax == rowMin))
{
//OE_NOTICE << "[osgEarth::GeoData] Exact value" << std::endl;
color = ia(px_i, py_i);
}
else if (colMax == colMin)
{
//OE_NOTICE << "[osgEarth::GeoData] Vertically" << std::endl;
//Linear interpolate vertically
for (unsigned int i = 0; i < 4; ++i)
{
color[i] = ((float)rowMax - py) * llColor[i] + (py - (float)rowMin) * ulColor[i];
}
}
else if (rowMax == rowMin)
{
//OE_NOTICE << "[osgEarth::GeoData] Horizontally" << std::endl;
//Linear interpolate horizontally
for (unsigned int i = 0; i < 4; ++i)
{
color[i] = ((float)colMax - px) * llColor[i] + (px - (float)colMin) * lrColor[i];
}
}
else
{
//OE_NOTICE << "[osgEarth::GeoData] Bilinear" << std::endl;
//Bilinear interpolate
float col1 = colMax - px, col2 = px - colMin;
float row1 = rowMax - py, row2 = py - rowMin;
for (unsigned int i = 0; i < 4; ++i)
{
float r1 = col1 * llColor[i] + col2 * lrColor[i];
float r2 = col1 * ulColor[i] + col2 * urColor[i];
//OE_INFO << "r1, r2 = " << r1 << " , " << r2 << std::endl;
color[i] = row1 * r1 + row2 * r2;
}
}
}
writer(color, c, r);
pixel++;
}
}
delete[] srcPointsX;
return result;
}
}
gdal 图片重投影：
GDALDataset*
createDataSetFromImage(const osg::Image* image, double minX, double minY, double maxX, double maxY, const std::string &projection)
{
//Clone the incoming image
osg::ref_ptr<osg::Image> clonedImage = new osg::Image(*image);
//Flip the image
clonedImage->flipVertical();
GDALDataType gdalDataType =
image->getDataType() == GL_UNSIGNED_BYTE ? GDT_Byte :
image->getDataType() == GL_UNSIGNED_SHORT ? GDT_UInt16 :
image->getDataType() == GL_FLOAT ? GDT_Float32 :
GDT_Byte;
int numBands =
image->getPixelFormat() == GL_RGBA ? 4 :
image->getPixelFormat() == GL_RGB ? 3 :
image->getPixelFormat() == GL_LUMINANCE ? 1 : 0;
if ( numBands == 0 )
{
OE_WARN << LC << "Failure in createDataSetFromImage: unsupported pixel format\n";
return 0L;
}
int pixelBytes =
gdalDataType == GDT_Byte ? numBands :
gdalDataType == GDT_UInt16 ? 2*numBands :
4*numBands;
GDALDataset* srcDS = createMemDS(image->s(), image->t(), numBands, gdalDataType, minX, minY, maxX, maxY, projection);
if ( srcDS )
{
CPLErr err = srcDS->RasterIO(
GF_Write,
0, 0,
clonedImage->s(), clonedImage->t(),
(void*)clonedImage->data(),
clonedImage->s(),
clonedImage->t(),
gdalDataType,
numBands,
NULL,
pixelBytes,
pixelBytes * image->s(),
1);
if ( err != CE_None )
{
OE_WARN << LC << "RasterIO failed.\n";
}
srcDS->FlushCache();
}
return srcDS;
}
osg::Image*
reprojectImage(osg::Image* srcImage, const std::string srcWKT, double srcMinX, double srcMinY, double srcMaxX, double srcMaxY,
const std::string destWKT, double destMinX, double destMinY, double destMaxX, double destMaxY,
int width = 0, int height = 0, bool useBilinearInterpolation = true)
{
GDAL_SCOPED_LOCK;
osg::Timer_t start = osg::Timer::instance()->tick();
//Create a dataset from the source image
GDALDataset* srcDS = createDataSetFromImage(srcImage, srcMinX, srcMinY, srcMaxX, srcMaxY, srcWKT);
OE_DEBUG << LC << "Source image is " << srcImage->s() << "x" << srcImage->t() << " in " << srcWKT << std::endl;
if (width == 0 || height == 0)
{
double outgeotransform[6];
double extents[4];
void* transformer = GDALCreateGenImgProjTransformer(srcDS, srcWKT.c_str(), NULL, destWKT.c_str(), 1, 0, 0);
GDALSuggestedWarpOutput2(srcDS,
GDALGenImgProjTransform, transformer,
outgeotransform,
&width,
&height,
extents,
0);
GDALDestroyGenImgProjTransformer(transformer);
}
OE_DEBUG << "Creating warped output of " << width <<"x" << height << " in " << destWKT << std::endl;
int numBands = srcDS->GetRasterCount();
GDALDataType dataType = srcDS->GetRasterBand(1)->GetRasterDataType();
GDALDataset* destDS = createMemDS(width, height, numBands, dataType, destMinX, destMinY, destMaxX, destMaxY, destWKT);
if (useBilinearInterpolation == true)
{
GDALReprojectImage(srcDS, NULL,
destDS, NULL,
GRA_Bilinear,
0,0,0,0,0);
}
else
{
GDALReprojectImage(srcDS, NULL,
destDS, NULL,
GRA_NearestNeighbour,
0,0,0,0,0);
}
osg::Image* result = createImageFromDataset(destDS);
delete srcDS;
delete destDS;
osg::Timer_t end = osg::Timer::instance()->tick();
OE_DEBUG << "Reprojected image in " << osg::Timer::instance()->delta_m(start,end) << std::endl;
return result;
}
再回到addlayer函数，由创建的geoImage生成Texture：
GeoImage geoImage = imageLayer->createImage( key, progress );
if ( geoImage.valid() )
{
if ( imageLayer->isCoverage() )
tex = createCoverageTexture(geoImage.getImage(), imageLayer);
else
tex = createImageTexture(geoImage.getImage(), imageLayer);
}
纹理存在，创建网格TerrainTileImageLayerModel：
if (tex)
{
tex->setName(model->getKey().str());
TerrainTileImageLayerModel* layerModel = new TerrainTileImageLayerModel();
layerModel->setImageLayer(imageLayer);
layerModel->setTexture(tex);
layerModel->setMatrix(new osg::RefMatrixf(textureMatrix));
model->colorLayers().push_back(layerModel);
if (imageLayer->isShared())
{
model->sharedLayers().push_back(layerModel);
}
if (imageLayer->isDynamic())
{
model->setRequiresUpdateTraverse(true);
}
}
返回到：createTileModel函数,添加地形
if ( requirements == 0L || requirements->elevationTexturesRequired() )
{
unsigned border = requirements->elevationBorderRequired() ? 1u : 0u;
addElevation( model.get(), map, key, filter, border, progress );
}
1047
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
