|
osgEarth 2.1.1
|
|
osgEarth 2.1.1
|
#include <osgEarth/GeoData>#include <osgEarth/ImageUtils>#include <osgEarth/Registry>#include <osgEarth/Cube>#include <osg/Notify>#include <osg/Timer>#include <gdal_priv.h>#include <gdalwarper.h>#include <ogr_spatialref.h>#include <memory.h>#include <sstream>#include <iomanip>
Include dependency graph for GeoData.cpp:Go to the source code of this file.
Defines | |
| #define | LC "[GeoData] " |
| #define | LC "[Geoid] " |
Functions | |
| static osg::Image * | createImageFromDataset (GDALDataset *ds) |
| static GDALDataset * | createMemDS (int width, int height, double minX, double minY, double maxX, double maxY, const std::string &projection) |
| static GDALDataset * | createDataSetFromImage (const osg::Image *image, double minX, double minY, double maxX, double maxY, const std::string &projection) |
| static 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) |
| static osg::Image * | manualReproject (const osg::Image *image, const GeoExtent &src_extent, const GeoExtent &dest_extent, unsigned int width=0, unsigned int height=0) |
| #define LC "[GeoData] " |
Definition at line 1130 of file GeoData.cpp.
| #define LC "[Geoid] " |
Definition at line 1130 of file GeoData.cpp.
| static GDALDataset* createDataSetFromImage | ( | const osg::Image * | image, |
| double | minX, | ||
| double | minY, | ||
| double | maxX, | ||
| double | maxY, | ||
| const std::string & | projection | ||
| ) | [static] |
Definition at line 660 of file GeoData.cpp.
{
//Clone the incoming image
osg::ref_ptr<osg::Image> clonedImage = new osg::Image(*image);
//Flip the image
clonedImage->flipVertical();
GDALDataset* srcDS = createMemDS(image->s(), image->t(), minX, minY, maxX, maxY, projection);
//Write the image data into the memory dataset
//If the image is already RGBA, just read all 4 bands in one call
if (image->getPixelFormat() == GL_RGBA)
{
srcDS->RasterIO(GF_Write, 0, 0, clonedImage->s(), clonedImage->t(), (void*)clonedImage->data(), clonedImage->s(), clonedImage->t(), GDT_Byte, 4, NULL, 4, 4 * image->s(), 1);
}
else if (image->getPixelFormat() == GL_RGB)
{
//OE_NOTICE << "[osgEarth::GeoData] Reprojecting RGB " << std::endl;
//Read the read, green and blue bands
srcDS->RasterIO(GF_Write, 0, 0, clonedImage->s(), clonedImage->t(), (void*)clonedImage->data(), clonedImage->s(), clonedImage->t(), GDT_Byte, 3, NULL, 3, 3 * image->s(), 1);
//Initialize the alpha values to 255.
unsigned char *alpha = new unsigned char[clonedImage->s() * clonedImage->t()];
memset(alpha, 255, clonedImage->s() * clonedImage->t());
GDALRasterBand* alphaBand = srcDS->GetRasterBand(4);
alphaBand->RasterIO(GF_Write, 0, 0, clonedImage->s(), clonedImage->t(), alpha, clonedImage->s(),clonedImage->t(), GDT_Byte, 0, 0);
delete[] alpha;
}
else
{
OE_WARN << LC << "createDataSetFromImage: unsupported pixel format " << std::hex << image->getPixelFormat() << std::endl;
}
srcDS->FlushCache();
return srcDS;
}
Here is the call graph for this function: Here is the caller graph for this function:| static osg::Image* createImageFromDataset | ( | GDALDataset * | ds | ) | [static] |
Definition at line 608 of file GeoData.cpp.
{
// called internally -- GDAL lock not required
//Allocate the image
osg::Image *image = new osg::Image;
image->allocateImage(ds->GetRasterXSize(), ds->GetRasterYSize(), 1, GL_RGBA, GL_UNSIGNED_BYTE);
ds->RasterIO(GF_Read, 0, 0, image->s(), image->t(), (void*)image->data(), image->s(), image->t(), GDT_Byte, 4, NULL, 4, 4 * image->s(), 1);
ds->FlushCache();
image->flipVertical();
return image;
}
Here is the caller graph for this function:| static GDALDataset* createMemDS | ( | int | width, |
| int | height, | ||
| double | minX, | ||
| double | minY, | ||
| double | maxX, | ||
| double | maxY, | ||
| const std::string & | projection | ||
| ) | [static] |
Definition at line 625 of file GeoData.cpp.
{
//Get the MEM driver
GDALDriver* memDriver = (GDALDriver*)GDALGetDriverByName("MEM");
if (!memDriver)
{
OE_NOTICE << "[osgEarth::GeoData] Could not get MEM driver" << std::endl;
}
//Create the in memory dataset.
GDALDataset* ds = memDriver->Create("", width, height, 4, GDT_Byte, 0);
//Initialize the color interpretation
ds->GetRasterBand(1)->SetColorInterpretation(GCI_RedBand);
ds->GetRasterBand(2)->SetColorInterpretation(GCI_GreenBand);
ds->GetRasterBand(3)->SetColorInterpretation(GCI_BlueBand);
ds->GetRasterBand(4)->SetColorInterpretation(GCI_AlphaBand);
//Initialize the geotransform
double geotransform[6];
double x_units_per_pixel = (maxX - minX) / (double)width;
double y_units_per_pixel = (maxY - minY) / (double)height;
geotransform[0] = minX;
geotransform[1] = x_units_per_pixel;
geotransform[2] = 0;
geotransform[3] = maxY;
geotransform[4] = 0;
geotransform[5] = -y_units_per_pixel;
ds->SetGeoTransform(geotransform);
ds->SetProjection(projection.c_str());
return ds;
}
Here is the caller graph for this function:| static osg::Image* manualReproject | ( | const osg::Image * | image, |
| const GeoExtent & | src_extent, | ||
| const GeoExtent & | dest_extent, | ||
| unsigned int | width = 0, |
||
| unsigned int | height = 0 |
||
| ) | [static] |
Definition at line 752 of file GeoData.cpp.
{
//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());
}
// need to know this in order to choose the right interpolation algorithm
const bool isSrcContiguous = src_extent.getSRS()->isContiguous();
osg::Image *result = new osg::Image();
result->allocateImage(width, height, 1, GL_RGBA, GL_UNSIGNED_BYTE);
//Initialize the image to be completely transparent
memset(result->data(), 0, result->getImageSizeInBytes());
ImageUtils::PixelReader ra(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, 0, true);
// 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);
if ( ! isSrcContiguous ) // non-contiguous space- use nearest neighbot
{
color = ia(px_i, py_i);
}
else // contiguous space - use bilinear sampling
{
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;
}
}
}
unsigned char* rgba = const_cast<unsigned char*>(ra.data(c,r,0));
rgba[0] = (unsigned char)(color.r() * 255);
rgba[1] = (unsigned char)(color.g() * 255);
rgba[2] = (unsigned char)(color.b() * 255);
rgba[3] = (unsigned char)(color.a() * 255);
pixel++;
}
}
delete[] srcPointsX;
return result;
}
Here is the call graph for this function: Here is the caller graph for this function:| static 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 |
||
| ) | [static] |
Definition at line 701 of file GeoData.cpp.
{
//OE_NOTICE << "Reprojecting..." << std::endl;
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 << "Source image is " << srcImage->s() << "x" << srcImage->t() << 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 << std::endl;
GDALDataset* destDS = createMemDS(width, height, destMinX, destMinY, destMaxX, destMaxY, destWKT);
GDALReprojectImage(srcDS, NULL,
destDS, NULL,
//GDALResampleAlg::GRA_NearestNeighbour,
GRA_Bilinear,
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;
}
Here is the call graph for this function: Here is the caller graph for this function:
1.7.3
