osgEarth and 16bit images?

classic Classic list List threaded Threaded
15 messages Options
MCA4213 MCA4213
Reply | Threaded
Open this post in threaded view
|

osgEarth and 16bit images?

Hi,
I have a question about how osgEarth treat the 16 bits images, even if it is an old question, but in my case when I try to read a Landsat-8 tiff image I got a white image.
<image driver="gdal" name="world-tiff" cache_enabled="false">
        <url>LC81220322015227LGN00_B3.TIF</url>
</image>


According to my understanding after testing different images, osgEarth is doing the following:
if (pixel <255) resultedColor=pixel;
else  resultedColor=255;

and that make it impossible to make any processing on the image even with the GLSL code.

So according to my basic thinking, because I am not expert, it should be just a normal division to convert a 16bit  or  a 32bit value to 255, and after that if we got a darker image we can use the GLSL code to stretch the color and got a good results.

My questions are: 
1- Is this situation (this strange conversion) done for a specific purpose or no?
2- which class or file can I find where the value is converted because I could not find it (I don't have much knowledge about how osgEarth work) ?

Thank you.
gwaldron gwaldron
Reply | Threaded
Open this post in threaded view
|

Re: osgEarth and 16bit images?

MCA,
I tested a 16-bit per channel GeoTIFF and it worked. Can you post the gdalinfo output of your image, or make a 16-bit image available for testing?
Glenn Waldron / Pelican Mapping
MCA4213 MCA4213
Reply | Threaded
Open this post in threaded view
|

Re: osgEarth and 16bit images?

Thanks, here is my  gdal info:

G:\DevTools\LibraryC\GDAL\bin\gdal\apps>gdalinfo.exe  "F:\Mes Etudes\4_China\Courses\Spring\Team Project\Images\Original\LC81220322015227LGN00_B3.TIF" -hist
Driver: GTiff/GeoTIFF
Files: F:\Mes Etudes\4_China\Courses\Spring\Team Project\Images\Original\LC81220322015227LGN00_B3.TIF
       F:\Mes Etudes\4_China\Courses\Spring\Team Project\Images\Original\LC81220322015227LGN00_MTL.txt
Size is 7681, 7821
Coordinate System is:
PROJCS["WGS 84 / UTM zone 50N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",117],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32650"]]
Origin = (493185.000000000000000,4582215.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Point
  METADATATYPE=ODL
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  493185.000, 4582215.000) (116d55' 6.54"E, 41d23'29.17"N)
Lower Left  (  493185.000, 4347585.000) (116d55'15.56"E, 39d16'39.33"N)
Upper Right (  723615.000, 4582215.000) (119d40'24.10"E, 41d21'37.45"N)
Lower Right (  723615.000, 4347585.000) (119d35'28.78"E, 39d14'55.60"N)
Center      (  608400.000, 4464900.000) (118d16'33.68"E, 40d19'39.36"N)
Band 1 Block=7681x1 Type=UInt16, ColorInterp=Gray
0...10...20...30...40...50...60...70...80...90...100 - done.
  256 buckets from -55.1078 to 28160.1:
  18448347 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 505 2694 6699 10542 23622 94679 311819 699552 1083479 1363365 1588793 1866149 2135746 2410293 2604210 2695373 2581752 2409911 2174837 1901064 1642073 1393807 1214538 1072700 988775 920034 851472 777082 703687 638574 565214 508073 459520 421918 391550 356334 328864 302983 280121 251344 224709 195758 168862 144589 119282 98480 81966 68341 58294 48192 40988 35730 30917 26584 23079 20591 17987 16054 14281 12659 11477 10297 9597 8576 7687 7127 6474 5998 5492 4949 4602 4255 3839 3339 3206 3021 2714 2517 2256 2137 1998 1725 1605 1529 1379 1277 1205 1205 1084 1015 860 843 781 685 703 610 609 603 513 522 485 434 384 394 373 394 311 305 302 289 259 268 251 256 210 220 215 198 203 206 164 175 171 167 155 160 151 146 141 120 127 133 139 114 113 121 114 103 118 90 99 69 93 91 101 93 90 90 89 105 81 84 64 77 65 58 69 54 41 62 62 56 64 50 53 50 51 48 53 46 42 51 45 31 32 43 30 33 44 34 34 34 29 30 28 40 30 31 39 22 33 38 38 38 23 27 1662
  Metadata:
    STATISTICS_MAXIMUM=28105
    STATISTICS_MEAN=5966.1452175149
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=4072.1814879297

and this is the histogram:




I tried with the green band of the Landsat-8 image in my case, my results was:


I don't know how to make the image available for you, but I can upload it to a google drive or mediafire if you want.
Thank you for you help.
dranders dranders
Reply | Threaded
Open this post in threaded view
|

Re: osgEarth and 16bit images?

In reply to this post by MCA4213
This is an issue with how the GDAL driver reads the images.  It reads into an array assuming they're stored as GDT_Bytes (the 8-bit representation), but GDALRasterBand::RasterIO does not automatically scale larger values, it simply clamps them to the range [0, 255].

From: http://www.gdal.org/gdal_tutorial.html

"The pData is the memory buffer the data is read into, or written from. It's real type must be whatever is passed as eBufType, such as GDT_Float32, or GDT_Byte. The RasterIO() call will take care of converting between the buffer's data type and the data type of the band. Note that when converting floating point data to integer RasterIO() rounds down, and when converting source values outside the legal range of the output the nearest legal value is used. This implies, for instance, that 16bit data read into a GDT_Byte buffer will map all values greater than 255 to 255, the data is not scaled!"

For Int16 images, you can either use the gdal tools to pre-process and scale your values down with reasonable limits to an 8-bit image, or you can modify the osgEarth GDAL driver to do that for you.  From working with a lot of these non-Byte image formats, I'll add it's sometimes hard to programmatically determine reasonable scaling limits.
dranders dranders
Reply | Threaded
Open this post in threaded view
|

Re: osgEarth and 16bit images?

One other option is to read your higher-bit data as an elevation layer, and then use the color ramp driver to convert it to grayscale or false-color.  When creating height fields, the osgEarth GDAL driver looks for the first GDT_Float32 band, but if that fails it just uses the first band.  The conversion between Int16 and Float32 should work fine without clamping.

 This isn't really a solution if you're trying to view multi-band imagery, but I've been using it pretty successfully.  My only quibble is it'd be nice to be able to adjust the color ramp dynamically without removing and re-adding the layer.
MCA4213 MCA4213
Reply | Threaded
Open this post in threaded view
|

Re: osgEarth and 16bit images?

Thank you very much dranders, your explanation for the problem was exactly the issue, and was very helpful for me.

I don't know if it is the same case for the RGB 16bit image, but for the gray uint16 all the value greater than 255 will be forced to 255 because,as dranders motioned before, of the GDT_Bytes in the RasterIO in the:
ReaderWriterGDAL.cpp --> osg::Image* createImage( const TileKey& key,  ProgressCallback*  progress)
-->  else if (bandGray) --> if (!*_options.interpolateImagery() || _options.interpolation() == INTERP_NEAREST)

Here osgEarth read all the images into GDT_Bytes even for the non 8 bit images (this is the problem).

So to test my idea, I just comment all the code in this section and replace it with the code from the next else, I put:
for (int r = 0; r < tileSize; ++r)
{
        double geoY   = ymin + (dy * (double)r);

        for (int c = 0; c < tileSize; ++c)
        {
                double geoX = xmin + (dx * (double)c);
                float  color = getInterpolatedValue(bandGray,geoX,geoY,false) / 255 ;
                //std::cout << "color=" <<color<< std::endl;
                *(image->data(c,r) + 0) = (unsigned char)color;
                *(image->data(c,r) + 1) = (unsigned char)color;
                *(image->data(c,r) + 2) = (unsigned char)color;
                if (bandAlpha != NULL)
                        *(image->data(c,r) + 3) = (unsigned char)getInterpolatedValue(bandAlpha,geoX,geoY,false);
                else
                        *(image->data(c,r) + 3) = 255;
        }
}
In the getInterpolatedValue function a float32 is used, so I just divide the value on 255 because the image was uint16.
The results was an image a little bit darker (my max value was 11000) but with a full detail, that I enhanced after that using the GLSL (stretching).

Of course it is not a solution, because it should be done for the different types, and also it can be enhanced if we do a stretching (we can got the mean max from the hist using gdal--> but take time)  to get a better results.

So this is just my point of view of this problem, I am not an expert, so what do you think about this test?

Thank you.
gwaldron gwaldron
Reply | Threaded
Open this post in threaded view
|

Re: osgEarth and 16bit images?

Is the real issue here the fact that this is not visible-spectrum RGB data? That is requires some interpretation in order to draw it as a raster layer? Could it be possible to use a VRT to interpret the data and "transform" it into visible spectrum?
Glenn Waldron / Pelican Mapping
MCA4213 MCA4213
Reply | Threaded
Open this post in threaded view
|

Re: osgEarth and 16bit images?

First thanks for the reply gwaldron.

The issue is that RasterIO is not able to scale the data with GDT_Bytes, let take the example of a 16bit image containing 4 classes of values (100, 1000, 5000, 10000), the gdal function will return 2 classes (100, 255, 255, 255: any value greater than 255 will equal to 255), so the use of this function for all kinds is wrong (according to the gdal doc).
if you have a 16bit gray image and you get it displayed, it means that in your image all the values are between 0 and 255, or most of them, because of that it looks like it is working (I am not sure).

I think it is better to use GDT_Float32 to get a float or a double values of the pixel as it is done for the elevation layers (I think the code above that I tested, refers to how to read the elev), then divide the value according to the type of the image (manual scale):
if( image is 8bit  ) float  color = getInterpolatedValue(bandGray,geoX,geoY,false);
if( image is 16bit ) float  color = getInterpolatedValue(bandGray,geoX,geoY,false)/255;
// for the same  example  (100, 1000, 5000, 10000) we will get (0, 4, 20, 40)
if( image is 32bit ) float  color = getInterpolatedValue(bandGray,geoX,geoY,false)/65025;

I test with 16 and 32 bit gray image and it worked, I got the image displayed with the detail inside (not a white image).

I don't know about vrt, so I can not help for that until I read on it, but using GDT_Float32  instead of the GDT_Bytes will fix the issue (just a proposition).
Thank you.
gwaldron gwaldron
Reply | Threaded
Open this post in threaded view
|

Re: osgEarth and 16bit images?

MCA,
I understand all that. What I am asking is this: while scaling (by dividing the value by 255 for 16-bit and by 255*255 for 32-bit) will certainly give you values in the 0..255 range, are those actually always the values you want? Is that kind of linear scaling appropriate in all circumstances, or is it just a "best guess" and is a "best guess" good enough?
Glenn Waldron / Pelican Mapping
dranders dranders
Reply | Threaded
Open this post in threaded view
|

Re: osgEarth and 16bit images?

In reply to this post by MCA4213
Glenn isn't wrong, you could use a VRT with an appropriate scale and offset to convert all your 16-bit values into an 8-bit range when read.  From an osgEarth::ImageLayer perspective, it's not like your screen can render higher-order bits anyway, so why support them?

MCA, your pseudo-code might get the job done for some images, but it's also very slow!  getInterpolatedValue performs bi-lienar interpolation and reads data out one pixel at a time.  Also, a lot of satellite imagery is actually only 12-bits per channel, so you're going to end up with lousy dynamic-range if you assume the full 16-bit range.

For that reason, I would probably have some additional options for controlling what the GDAL driver should do with non GDT_Byte images.  I would have a flag for enabling value scaling, and then a minValueRange and maxValueRange.  The code would look something like this:

uint8_t* imageBuffer = new uint8_t[target_width * target_height]
DataType* tempBuffer = new DataType[nRows * nCols];
rasterBand->RasterIO( GF_Read, off_x, off_y, width, height, tempBuffer, target_width, target_height, GDT_DataType, 0, 0);
scaleBufferValues( tempBuffer, imageBuffer, target_width*target_height, minValueRange, maxValueRange );

with function (probably templated on DataType):

void scaleBufferValues( DataType* tempBuffer, uint8_t* imageBuffer, size_t numElem, DataType minValue, DataType maxValue )
{
  for ( size_t i=0; i < numElem; ++i )
  {
    imageBuffer[i] = (uint8_t) ((tempBuffer[i] - minValue)/(maxValue-minValue));
   }
}

The beauty of OpenSource is you can change the GDAL driver yourself to implement this, or create your own forked version!

Dean
MCA4213 MCA4213
Reply | Threaded
Open this post in threaded view
|

Re: osgEarth and 16bit images?

You are right, dividing on 255 or 255*255 will results on loosing some data.
If we take the Landsat-8 images they have 14bit values, so dividing on 255 will results in a little bit black image with loosing some details. I used before the min max after that using the GLSL code :
<glsl>
                float NMIN = 0.0;
                float NMAX = 255.0;
                float OMIN = 40.0; //from histogram /255 (after scaling)
                float OMAX = 175.0; //from histogram /255
                float space = ( NMAX - NMIN ) / ( OMAX - OMIN ) ;
                float bins  = ( OMAX - OMIN );
                float value = color.r * 255.0 - OMIN;
                color.rgb = vec3(( NMIN + ( value * space ))/255.0 );
        </glsl>
BUT I already loose the detail when I scaled the image by dividing on 255.


Your proposition, dranders, seems to be good, because the function that I used read the pixels one by one so it is slow, and you are proposing a parametric scale.
 
The reason that I divided on 255 and not using min max here, is because the only way that I know to get the mean max is by calling the hist in GDAL, and it takes some seconds to get the results, so I think it will slow the GDAL plugin ( I didn't compare the times, I just test gdal_info -hist).
If there is a better way to get the mean and the max faster, then scaling according this values will be certainly the best choice that give us an image with the max possible details.
Thank you.
dranders dranders
Reply | Threaded
Open this post in threaded view
|

Re: osgEarth and 16bit images?

MCA4213 wrote
The reason that I divided on 255 and not using min max here, is because the only way that I know to get the mean max is by calling the hist in GDAL, and it takes some seconds to get the results, so I think it will slow the GDAL plugin ( I didn't compare the times, I just test gdal_info -hist).
If there is a better way to get the mean and the max faster, then scaling according this values will be certainly the best choice that give us an image with the max possible details.
Thank you.
You only need to call RasterBand::GetHistogram or ComputeStatistics when initializing the GDALTileSource, not during every call to createImage.  So it's a one time hit when you load the driver, and shouldn't slow down the viewer.  And you only need to use those methods if you aren't directly specifying scaling ranges through the options!
jasonbeverage jasonbeverage
Reply | Threaded
Open this post in threaded view
|

Re: osgEarth and 16bit images?

Just to add a bit more to this, you can use gdal_translate along with the "-scale" option to generate a VRT that you can use.

You can either specify explicit source and destination data ranges or if you omit then entirely then GDAL will do it's best to come up with something sensible based on the data itself.


For example,

gdal_translate -of VRT -scale input.tif output.vrt

If you are happy with the results you can also output a TIF instead of a VRT.  VRT is just really convenient b/c you don't need to actually process very large images just to see a small change like scaling data.

Jason


On Fri, Jul 22, 2016 at 9:06 AM dranders [via osgEarth] <[hidden email]> wrote:
MCA4213 wrote
The reason that I divided on 255 and not using min max here, is because the only way that I know to get the mean max is by calling the hist in GDAL, and it takes some seconds to get the results, so I think it will slow the GDAL plugin ( I didn't compare the times, I just test gdal_info -hist).
If there is a better way to get the mean and the max faster, then scaling according this values will be certainly the best choice that give us an image with the max possible details.
Thank you.
You only need to call RasterBand::GetHistogram or ComputeStatistics when initializing the GDALTileSource, not during every call to createImage.  So it's a one time hit when you load the driver, and shouldn't slow down the viewer.  And you only need to use those methods if you aren't directly specifying scaling ranges through the options!



If you reply to this email, your message will be added to the discussion below:
To start a new topic under osgEarth, email [hidden email]
To unsubscribe from osgEarth, click here.
NAML
MCA4213 MCA4213
Reply | Threaded
Open this post in threaded view
|

Re: osgEarth and 16bit images?

Thanks for all your replies,

Sorry I am a little bit out this days because of the vacation, but lest discuss some points:

jasonbeverage proposition: first thanks for the advise.
Using the gdal_translate will make osgEarth a non-real time rendering toolkit for a large numbers of images that others software can deal with them easily, I use it before in one of my application for viewing images (browsing a hundred of satellite images), and I need to wait 1 minutes or more to get only 1 image converted each time: really a loose of time especially if the image is shared on a small network without any services.
--> using this conversion is only useful if you have 1 image, if not, if you are creating an app to deal with a lot of images,  it will be really a very bad choice.

dranders proposition:
I think you are right, I will find the initialization part and calculate the mean/max there, then I will tray to change the GDT_Byte to the GDT_DataType and then I will create a scale function as you mentioned before.



I still have some questions, because I don't know if I should activate the scale on the all  parts of creatImage function, so:
* what is gray coverage?? (the non-coverage I think is the normal image).
* should the coverage image be scaled?
* should the scale be done on Nearest interpolation and Exact point? (I think it should be done on both)

I don't see a reason to maintain the plugin as it is now (not able to deal with all images), and as I said before osgEarth is a real time rendering toolkit, so it must be able to do this kind of processing on the fly (we can add an option --scale to activate this).
I will tray to find some free time soon then I will publish my changes to discuss them.

Thank you very much.
MCA4213 MCA4213
Reply | Threaded
Open this post in threaded view
|

Re: osgEarth and 16bit images?

In reply to this post by jasonbeverage
Hi every one,
I did some modification for the GDAL driver so now it can display the 8bit/16bit and 32 bit images on the fly.
I am not expert so I just convert the world.tif to 16bit and 32bit image for testing. also I tested the Landsat-8 image that  the values range is on 12bit.

The principle that I took is as it was explained by "dranders", with some enhancement to work in these cases, so:
1- we need to calculate the min/max in the initialization function;
2- get the type of image (in create image function for RGB and Gray)
2- create a temp buffer according to the type of image;

GDALDataType gdalDataType = bandGray->GetRasterDataType();
void* tempGray,*tempAlpha;
if (gdalDataType == GDT_Byte)
{
        tempGray = new unsigned char[target_width * target_height];
        tempAlpha = new unsigned char[target_width * target_height];
}
else
        if (gdalDataType == GDT_Int16 || gdalDataType == GDT_UInt16)
        {
                tempGray = new int16_t[target_width * target_height];
                tempAlpha = new int16_t[target_width * target_height];
        }
        else
                if (gdalDataType == GDT_Int32 || gdalDataType == GDT_UInt32)
                {
                        tempGray = new int32_t[target_width * target_height];
                        tempAlpha = new int32_t[target_width * target_height];
                }
4- get the value using RasterIO but according to each type:
bandGray->RasterIO(GF_Read, off_x, off_y, width, height, tempGray, target_width, target_height, gdalDataType, 0, 0);

5- scaling the results:
void scaleBufferValues(void* tempBuffer, uint8_t* &imageBuffer, size_t numElem, int minValue, int maxValue,GDALDataType gdalDataType)
{
        if (gdalDataType==GDT_Byte)
                for (size_t i = 0; i < numElem; ++i)
                {
                        imageBuffer[i] = (uint8_t)((((unsigned char*)(tempBuffer))[i] - minValue)* 255 / (maxValue - minValue));
                }
        else
                if (gdalDataType == GDT_Int16 || gdalDataType == GDT_UInt16)
                for (size_t i = 0; i < numElem; ++i)
                {
                        imageBuffer[i] = (uint8_t)((((int16_t*)(tempBuffer))[i] - minValue) * 255 / (maxValue - minValue));
                }
                else
                        if (gdalDataType == GDT_Int32 || gdalDataType == GDT_UInt32)
                                for (size_t i = 0; i < numElem; ++i)
                                {
                                        imageBuffer[i] = (uint8_t)((double)(((int32_t*)(tempBuffer))[i] - minValue) * (double)255 / (double)(maxValue - minValue));
                                }
}


with this code no need to use gdal_translate anymore, it can be done on the fly for gray images and RGB images.
The final file is uploaded here:
ReaderWriterGDAL.cpp

Any suggestions, or remarks???