The tool is called rgb2pct.py. But it uses some internal GDAL code. I
attached the patched versions of the sourcefiles. But be careful. The
files might have changed since the last time I used it (> 1 year)
If GDAL has to rotate your file when referencing it it has to fill the
triangles at the borders caused by the rotation. You define a
transparent color for that if you want to stitch several files. Such a
32bit color file with a transparent channel, could not be handled back
in these days. But maybe the problem is already solved by GDAL. I
haven't had a look lately.
Am 07.10.2015 um 15:53 schrieb Flo:
>
> On 10/07/15 15:28, Oliver Eichler wrote:
>> Hi Flo,
>>
>> you need either 32bit colors or all files must have the same color
>> table. If referencing does not scale/rotate the file you can apply a
>> common table right before the referencing with a normal image
>> manipulation software like Gimp. If that is not possible you can expand
>> the files to 32bit color and do the referencing. If that results into
>> files too large, you can extract referencing information, reduce the
>> colors to a common color table with Gimp, and add the information again.
>>
> Hm, file size gets much bigger, that's right if I expand the table with
> gdal_translate. Now I use compression like deflate to reduce that effect.
>
> If this is becoming a problem for me I will go into Gimp to see if I can
> do these operations from the command line in order to integrate that in
> a loop/script.
>
>> There is even a GDAL tool to reduce colors. However it does not take a
>> transparent channel into account. Years ago I patched it. But that was
>> quite quick and dirty. It never made it into
>> the official code.
>>
> How is this tool called? I don't think I have to deal with a transparent
> channel now. Thus, I could work with it without a patch.
>
> Best Regards,
> Flo.
>
>> Am 07.10.2015 um 14:55 schrieb Flo:
>>> On 10/06/15 17:45, Oliver Eichler wrote:
>>>> Hi Flo,
>>>>
>>>> the file has reference points (GCP) but no reference information as such:
>>>>
>>>> Corner Coordinates:
>>>> Upper Left ( 0.0, 0.0)
>>>> Lower Left ( 0.0, 2336.0)
>>>> Upper Right ( 3125.0, 0.0)
>>>> Lower Right ( 3125.0, 2336.0)
>>>> Center ( 1562.5, 1168.0)
>>>>
>>>>
>>>> This has to look a bit like:
>>>>
>>>> Corner Coordinates:
>>>> Upper Left ( 1193772.946, 5877433.181) ( 10d43'25.84"E, 46d47'40.14"N)
>>>> Lower Left ( 1193772.946, 5756585.800) ( 10d43'25.84"E, 46d 2'37.39"N)
>>>> Upper Right ( 1400769.640, 5877433.181) ( 12d34'59.98"E, 46d47'40.14"N)
>>>> Lower Right ( 1400769.640, 5756585.800) ( 12d34'59.98"E, 46d 2'37.39"N)
>>>> Center ( 1297271.293, 5817009.490) ( 11d39'12.91"E, 46d25'13.45"N)
>>>>
>>> Thank you. That is very valueable for me.
>>>
>>>> Try gdalwarp instead of gdal_translate.
>>>>
>>> Hit!!! It's working with gdalwrap. However I have to postprocess the
>>> output with gdal_translate to expand the color table if I build a mapset
>>> of several images. At least I haven't found an option to do so with
>>> gdalwarp.
>>>
>>> Thank you very much!
>>>
>>> Best Regards,
>>> Flo.
>>>
>>> ------------------------------------------------------------------------------
>>> Full-scale, agent-less Infrastructure Monitoring from a single dashboard
>>> Integrate with 40+ ManageEngine ITSM Solutions for complete visibility
>>> Physical-Virtual-Cloud Infrastructure monitoring from one console
>>> Real user monitoring with APM Insights and performance trend reports
>>> Learn More http://pubads.g.doubleclick.net/gampad/clk?id=247754911&iu=/4140
>>> _______________________________________________
>>> Qlandkartegt-users mailing list
>>> [email protected]
>>> https://lists.sourceforge.net/lists/listinfo/qlandkartegt-users
>>
>> ------------------------------------------------------------------------------
>> Full-scale, agent-less Infrastructure Monitoring from a single dashboard
>> Integrate with 40+ ManageEngine ITSM Solutions for complete visibility
>> Physical-Virtual-Cloud Infrastructure monitoring from one console
>> Real user monitoring with APM Insights and performance trend reports
>> Learn More http://pubads.g.doubleclick.net/gampad/clk?id=247754911&iu=/4140
>> _______________________________________________
>> Qlandkartegt-users mailing list
>> [email protected]
>> https://lists.sourceforge.net/lists/listinfo/qlandkartegt-users
>>
> ------------------------------------------------------------------------------
> Full-scale, agent-less Infrastructure Monitoring from a single dashboard
> Integrate with 40+ ManageEngine ITSM Solutions for complete visibility
> Physical-Virtual-Cloud Infrastructure monitoring from one console
> Real user monitoring with APM Insights and performance trend reports
> Learn More http://pubads.g.doubleclick.net/gampad/clk?id=247754911&iu=/4140
> _______________________________________________
> Qlandkartegt-users mailing list
> [email protected]
> https://lists.sourceforge.net/lists/listinfo/qlandkartegt-users
/******************************************************************************
* $Id: gdaldither.cpp 13342 2007-12-14 20:58:31Z rouault $
*
* Project: CIETMap Phase 2
* Purpose: Convert RGB (24bit) to a pseudo-colored approximation using
* Floyd-Steinberg dithering (error diffusion).
* Author: Frank Warmerdam, [email protected]
*
******************************************************************************
* Copyright (c) 2001, Frank Warmerdam
*
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the "Software"),
* to deal in the Software without restriction, including without limitation
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
* and/or sell copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
* DEALINGS IN THE SOFTWARE.
******************************************************************************
*
* Notes:
*
* [1] Floyd-Steinberg dither:
* I should point out that the actual fractions we used were, assuming
* you are at X, moving left to right:
*
* X 7/16
* 3/16 5/16 1/16
*
* Note that the error goes to four neighbors, not three. I think this
* will probably do better (at least for black and white) than the
* 3/8-3/8-1/4 distribution, at the cost of greater processing. I have
* seen the 3/8-3/8-1/4 distribution described as "our" algorithm before,
* but I have no idea who the credit really belongs to.
* --
* Lou Steinberg
*
*/
#include "gdal_priv.h"
#include "gdal_alg.h"
CPL_CVSID("$Id: gdaldither.cpp 13342 2007-12-14 20:58:31Z rouault $");
#define C_LEVELS 32
static void FindNearestColor( int nColors, int *panPCT, GByte *pabyColorMap );
/************************************************************************/
/* GDALDitherRGB2PCT() */
/************************************************************************/
/**
* 24bit to 8bit conversion with dithering.
*
* This functions utilizes Floyd-Steinberg dithering in the process of
* converting a 24bit RGB image into a pseudocolored 8bit image using a
* provided color table.
*
* The red, green and blue input bands do not necessarily need to come
* from the same file, but they must be the same width and height. They will
* be clipped to 8bit during reading, so non-eight bit bands are generally
* inappropriate. Likewise the hTarget band will be written with 8bit values
* and must match the width and height of the source bands.
*
* The color table cannot have more than 256 entries.
*
* @param hRed Red input band.
* @param hGreen Green input band.
* @param hBlue Blue input band.
* @param hTarget Output band.
* @param hColorTable the color table to use with the output band.
* @param pfnProgress callback for reporting algorithm progress matching the
* GDALProgressFunc() semantics. May be NULL.
* @param pProgressArg callback argument passed to pfnProgress.
*
* @return CE_None on success or CE_Failure if an error occurs.
*/
int CPL_STDCALL
GDALDitherRGB2PCT( GDALRasterBandH hRed,
GDALRasterBandH hGreen,
GDALRasterBandH hBlue,
GDALRasterBandH hTarget,
GDALColorTableH hColorTable,
GDALProgressFunc pfnProgress,
void * pProgressArg )
{
VALIDATE_POINTER1( hRed, "GDALDitherRGB2PCT", CE_Failure );
VALIDATE_POINTER1( hGreen, "GDALDitherRGB2PCT", CE_Failure );
VALIDATE_POINTER1( hBlue, "GDALDitherRGB2PCT", CE_Failure );
VALIDATE_POINTER1( hTarget, "GDALDitherRGB2PCT", CE_Failure );
VALIDATE_POINTER1( hColorTable, "GDALDitherRGB2PCT", CE_Failure );
int nXSize, nYSize;
CPLErr err = CE_None;
int hasNoData = 0;
int iNoDataR = -1;
int iNoDataG = -1;
int iNoDataB = -1;
int iNoData = GDALGetRasterNoDataValue(hTarget, &hasNoData);
if(hasNoData)
{
iNoDataR = GDALGetRasterNoDataValue(hRed, 0);
iNoDataG = GDALGetRasterNoDataValue(hGreen, 0);
iNoDataB = GDALGetRasterNoDataValue(hBlue, 0);
}
/* -------------------------------------------------------------------- */
/* Validate parameters. */
/* -------------------------------------------------------------------- */
nXSize = GDALGetRasterBandXSize( hRed );
nYSize = GDALGetRasterBandYSize( hRed );
if( GDALGetRasterBandXSize( hGreen ) != nXSize
|| GDALGetRasterBandYSize( hGreen ) != nYSize
|| GDALGetRasterBandXSize( hBlue ) != nXSize
|| GDALGetRasterBandYSize( hBlue ) != nYSize )
{
CPLError( CE_Failure, CPLE_IllegalArg,
"Green or blue band doesn't match size of red band.\n" );
return CE_Failure;
}
if( GDALGetRasterBandXSize( hTarget ) != nXSize
|| GDALGetRasterBandYSize( hTarget ) != nYSize )
{
CPLError( CE_Failure, CPLE_IllegalArg,
"GDALDitherRGB2PCT(): "
"Target band doesn't match size of source bands.\n" );
return CE_Failure;
}
if( pfnProgress == NULL )
pfnProgress = GDALDummyProgress;
/* -------------------------------------------------------------------- */
/* Setup more direct colormap. */
/* -------------------------------------------------------------------- */
int nColors, anPCT[768], iColor, maxColor;
nColors = GDALGetColorEntryCount( hColorTable );
if(hasNoData)
{
nColors -= 1;
}
maxColor = nColors - 1;
if (nColors == 0 )
{
CPLError( CE_Failure, CPLE_IllegalArg,
"GDALDitherRGB2PCT(): "
"Color table must not be empty.\n" );
return CE_Failure;
}
else if (nColors > 256)
{
CPLError( CE_Failure, CPLE_IllegalArg,
"GDALDitherRGB2PCT(): "
"Color table cannot have more than 256 entries.\n" );
return CE_Failure;
}
for( iColor = 0; iColor < nColors; iColor++ )
{
GDALColorEntry sEntry;
GDALGetColorEntryAsRGB( hColorTable, iColor, &sEntry );
anPCT[iColor ] = sEntry.c1;
anPCT[iColor+256] = sEntry.c2;
anPCT[iColor+512] = sEntry.c3;
}
/* -------------------------------------------------------------------- */
/* Build a 24bit to 8 bit color mapping. */
/* -------------------------------------------------------------------- */
GByte *pabyColorMap;
pabyColorMap = (GByte *) CPLMalloc(C_LEVELS * C_LEVELS * C_LEVELS
* sizeof(int));
FindNearestColor( nColors, anPCT, pabyColorMap );
/* -------------------------------------------------------------------- */
/* Setup various variables. */
/* -------------------------------------------------------------------- */
GByte *pabyRed, *pabyGreen, *pabyBlue, *pabyIndex;
int *panError;
pabyRed = (GByte *) VSIMalloc(nXSize);
pabyGreen = (GByte *) VSIMalloc(nXSize);
pabyBlue = (GByte *) VSIMalloc(nXSize);
pabyIndex = (GByte *) VSIMalloc(nXSize);
panError = (int *) VSICalloc(sizeof(int),(nXSize+2) * 3);
if (pabyRed == NULL ||
pabyGreen == NULL ||
pabyBlue == NULL ||
pabyIndex == NULL ||
panError == NULL)
{
CPLError( CE_Failure, CPLE_OutOfMemory,
"VSIMalloc(): Out of memory in GDALDitherRGB2PCT" );
err = CE_Failure;
goto end_and_cleanup;
}
/* ==================================================================== */
/* Loop over all scanlines of data to process. */
/* ==================================================================== */
int iScanline;
for( iScanline = 0; iScanline < nYSize; iScanline++ )
{
int nLastRedError, nLastGreenError, nLastBlueError, i;
/* -------------------------------------------------------------------- */
/* Report progress */
/* -------------------------------------------------------------------- */
if( !pfnProgress( iScanline / (double) nYSize, NULL, pProgressArg ) )
{
CPLError( CE_Failure, CPLE_UserInterrupt, "User Terminated" );
err = CE_Failure;
goto end_and_cleanup;
}
/* -------------------------------------------------------------------- */
/* Read source data. */
/* -------------------------------------------------------------------- */
GDALRasterIO( hRed, GF_Read, 0, iScanline, nXSize, 1,
pabyRed, nXSize, 1, GDT_Byte, 0, 0 );
GDALRasterIO( hGreen, GF_Read, 0, iScanline, nXSize, 1,
pabyGreen, nXSize, 1, GDT_Byte, 0, 0 );
GDALRasterIO( hBlue, GF_Read, 0, iScanline, nXSize, 1,
pabyBlue, nXSize, 1, GDT_Byte, 0, 0 );
/* -------------------------------------------------------------------- */
/* Apply the error from the previous line to this one. */
/* -------------------------------------------------------------------- */
for( i = 0; i < nXSize; i++ )
{
if(pabyRed[i] != iNoDataR)
{
pabyRed[i] = (GByte) MAX(0,MIN(maxColor,(pabyRed[i] + panError[i*3+0+3])));
}
if(pabyGreen[i] != iNoDataG)
{
pabyGreen[i] = (GByte) MAX(0,MIN(maxColor,(pabyGreen[i] + panError[i*3+1+3])));
}
if(pabyBlue[i] != iNoDataB)
{
pabyBlue[i] = (GByte) MAX(0,MIN(maxColor,(pabyBlue[i] + panError[i*3+2+3])));
}
}
memset( panError, 0, sizeof(int) * (nXSize+2) * 3 );
/* -------------------------------------------------------------------- */
/* Figure out the nearest color to the RGB value. */
/* -------------------------------------------------------------------- */
nLastRedError = 0;
nLastGreenError = 0;
nLastBlueError = 0;
for( i = 0; i < nXSize; i++ )
{
int iIndex, nError, nSixth, iRed, iGreen, iBlue;
int nRedValue, nGreenValue, nBlueValue;
if(pabyRed[i] != iNoDataR || pabyGreen[i] != iNoDataG || pabyBlue[i] != iNoDataB)
{
nRedValue = MAX(0,MIN(maxColor, pabyRed[i] + nLastRedError));
nGreenValue = MAX(0,MIN(maxColor, pabyGreen[i] + nLastGreenError));
nBlueValue = MAX(0,MIN(maxColor, pabyBlue[i] + nLastBlueError));
iRed = nRedValue * C_LEVELS / 256;
iGreen = nGreenValue * C_LEVELS / 256;
iBlue = nBlueValue * C_LEVELS / 256;
iIndex = pabyColorMap[iRed + iGreen * C_LEVELS
+ iBlue * C_LEVELS * C_LEVELS];
pabyIndex[i] = (GByte) iIndex;
}
else
{
pabyIndex[i] = (GByte) iNoData;
continue;
}
/* -------------------------------------------------------------------- */
/* Compute Red error, and carry it on to the next error line. */
/* -------------------------------------------------------------------- */
nError = nRedValue - anPCT[iIndex ];
nSixth = nError / 6;
panError[i*3 ] += nSixth;
panError[i*3+6 ] = nSixth;
panError[i*3+3 ] += nError - 5 * nSixth;
nLastRedError = 2 * nSixth;
/* -------------------------------------------------------------------- */
/* Compute Green error, and carry it on to the next error line. */
/* -------------------------------------------------------------------- */
nError = nGreenValue - anPCT[iIndex+256];
nSixth = nError / 6;
panError[i*3 +1] += nSixth;
panError[i*3+6+1] = nSixth;
panError[i*3+3+1] += nError - 5 * nSixth;
nLastGreenError = 2 * nSixth;
/* -------------------------------------------------------------------- */
/* Compute Blue error, and carry it on to the next error line. */
/* -------------------------------------------------------------------- */
nError = nBlueValue - anPCT[iIndex+512];
nSixth = nError / 6;
panError[i*3 +2] += nSixth;
panError[i*3+6+2] = nSixth;
panError[i*3+3+2] += nError - 5 * nSixth;
nLastBlueError = 2 * nSixth;
}
/* -------------------------------------------------------------------- */
/* Write results. */
/* -------------------------------------------------------------------- */
GDALRasterIO( hTarget, GF_Write, 0, iScanline, nXSize, 1,
pabyIndex, nXSize, 1, GDT_Byte, 0, 0 );
}
pfnProgress( 1.0, NULL, pProgressArg );
/* -------------------------------------------------------------------- */
/* Cleanup */
/* -------------------------------------------------------------------- */
end_and_cleanup:
CPLFree( pabyRed );
CPLFree( pabyGreen );
CPLFree( pabyBlue );
CPLFree( pabyIndex );
CPLFree( panError );
CPLFree( pabyColorMap );
return err;
}
/************************************************************************/
/* FindNearestColor() */
/* */
/* Finear near PCT color for any RGB color. */
/************************************************************************/
static void FindNearestColor( int nColors, int *panPCT, GByte *pabyColorMap )
{
int iBlue, iGreen, iRed;
int iColor;
/* -------------------------------------------------------------------- */
/* Loop over all the cells in the high density cube. */
/* -------------------------------------------------------------------- */
for( iBlue = 0; iBlue < C_LEVELS; iBlue++ )
{
for( iGreen = 0; iGreen < C_LEVELS; iGreen++ )
{
for( iRed = 0; iRed < C_LEVELS; iRed++ )
{
int nRedValue, nGreenValue, nBlueValue;
int nBestDist = 768, nBestIndex = 0;
nRedValue = (iRed * 255) / (C_LEVELS-1);
nGreenValue = (iGreen * 255) / (C_LEVELS-1);
nBlueValue = (iBlue * 255) / (C_LEVELS-1);
for( iColor = 0; iColor < nColors; iColor++ )
{
int nThisDist;
nThisDist = ABS(nRedValue - panPCT[iColor ])
+ ABS(nGreenValue - panPCT[iColor+256])
+ ABS(nBlueValue - panPCT[iColor+512]);
if( nThisDist < nBestDist )
{
nBestIndex = iColor;
nBestDist = nThisDist;
}
}
pabyColorMap[iRed + iGreen*C_LEVELS
+ iBlue*C_LEVELS*C_LEVELS] = (GByte)nBestIndex;
}
}
}
}
/******************************************************************************
* $Id: gdalmediancut.cpp 20629 2010-09-16 19:08:40Z rouault $
*
* Project: CIETMap Phase 2
* Purpose: Use median cut algorithm to generate an near-optimal PCT for a
* given RGB image. Implemented as function GDALComputeMedianCutPCT.
* Author: Frank Warmerdam, [email protected]
*
******************************************************************************
* Copyright (c) 2001, Frank Warmerdam
*
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the "Software"),
* to deal in the Software without restriction, including without limitation
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
* and/or sell copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
* DEALINGS IN THE SOFTWARE.
******************************************************************************
*
* This code was based on the tiffmedian.c code from libtiff (www.libtiff.org)
* which was based on a paper by Paul Heckbert:
*
* "Color Image Quantization for Frame Buffer Display", Paul
* Heckbert, SIGGRAPH proceedings, 1982, pp. 297-307.
*
*/
#include "gdal_priv.h"
#include "gdal_alg.h"
CPL_CVSID("$Id: gdalmediancut.cpp 20629 2010-09-16 19:08:40Z rouault $");
#define MAX_CMAP_SIZE 256
#define COLOR_DEPTH 8
#define MAX_COLOR 256
#define GMC_B_DEPTH 5 /* # bits/pixel to use */
#define GMC_B_LEN (1L<<GMC_B_DEPTH)
#define COLOR_SHIFT (COLOR_DEPTH-GMC_B_DEPTH)
typedef struct colorbox {
struct colorbox *next, *prev;
int rmin, rmax;
int gmin, gmax;
int bmin, bmax;
int total;
} Colorbox;
static void splitbox(Colorbox* ptr, int (*histogram)[GMC_B_LEN][GMC_B_LEN],
Colorbox **pfreeboxes, Colorbox **pusedboxes);
static void shrinkbox(Colorbox* box, int (*histogram)[GMC_B_LEN][GMC_B_LEN]);
static Colorbox* largest_box(Colorbox *usedboxes);
/************************************************************************/
/* GDALComputeMedianCutPCT() */
/************************************************************************/
/**
* Compute optimal PCT for RGB image.
*
* This function implements a median cut algorithm to compute an "optimal"
* pseudocolor table for representing an input RGB image. This PCT could
* then be used with GDALDitherRGB2PCT() to convert a 24bit RGB image into
* an eightbit pseudo-colored image.
*
* This code was based on the tiffmedian.c code from libtiff (www.libtiff.org)
* which was based on a paper by Paul Heckbert:
*
* \verbatim
* "Color Image Quantization for Frame Buffer Display", Paul
* Heckbert, SIGGRAPH proceedings, 1982, pp. 297-307.
* \endverbatim
*
* The red, green and blue input bands do not necessarily need to come
* from the same file, but they must be the same width and height. They will
* be clipped to 8bit during reading, so non-eight bit bands are generally
* inappropriate.
*
* @param hRed Red input band.
* @param hGreen Green input band.
* @param hBlue Blue input band.
* @param pfnIncludePixel function used to test which pixels should be included
* in the analysis. At this time this argument is ignored and all pixels are
* utilized. This should normally be NULL.
* @param nColors the desired number of colors to be returned (2-256).
* @param hColorTable the colors will be returned in this color table object.
* @param pfnProgress callback for reporting algorithm progress matching the
* GDALProgressFunc() semantics. May be NULL.
* @param pProgressArg callback argument passed to pfnProgress.
*
* @return returns CE_None on success or CE_Failure if an error occurs.
*/
extern "C" int CPL_STDCALL
GDALComputeMedianCutPCT( GDALRasterBandH hRed,
GDALRasterBandH hGreen,
GDALRasterBandH hBlue,
int (*pfnIncludePixel)(int,int,void*),
int nColors,
GDALColorTableH hColorTable,
GDALProgressFunc pfnProgress,
void * pProgressArg )
{
VALIDATE_POINTER1( hRed, "GDALComputeMedianCutPCT", CE_Failure );
VALIDATE_POINTER1( hGreen, "GDALComputeMedianCutPCT", CE_Failure );
VALIDATE_POINTER1( hBlue, "GDALComputeMedianCutPCT", CE_Failure );
int nXSize, nYSize;
CPLErr err = CE_None;
int iNoDataOkR = 0;
int iNoDataOkG = 0;
int iNoDataOkB = 0;
int iNoDataR = GDALGetRasterNoDataValue(hRed, &iNoDataOkR);
int iNoDataG = GDALGetRasterNoDataValue(hGreen, &iNoDataOkG);
int iNoDataB = GDALGetRasterNoDataValue(hBlue, &iNoDataOkB);
int hasNoData = iNoDataOkR && iNoDataOkG && iNoDataOkB;
if(hasNoData)
{ // decrease nColors as the last entry will be the no data entry
nColors -= 1;
}
/* -------------------------------------------------------------------- */
/* Validate parameters. */
/* -------------------------------------------------------------------- */
nXSize = GDALGetRasterBandXSize( hRed );
nYSize = GDALGetRasterBandYSize( hRed );
if( GDALGetRasterBandXSize( hGreen ) != nXSize
|| GDALGetRasterBandYSize( hGreen ) != nYSize
|| GDALGetRasterBandXSize( hBlue ) != nXSize
|| GDALGetRasterBandYSize( hBlue ) != nYSize )
{
CPLError( CE_Failure, CPLE_IllegalArg,
"Green or blue band doesn't match size of red band.\n" );
return CE_Failure;
}
if( pfnIncludePixel != NULL )
{
CPLError( CE_Failure, CPLE_IllegalArg,
"GDALComputeMedianCutPCT() doesn't currently support "
" pfnIncludePixel function." );
return CE_Failure;
}
if ( nColors <= 0 )
{
CPLError( CE_Failure, CPLE_IllegalArg,
"GDALComputeMedianCutPCT() : nColors must be strictly greater than 1." );
return CE_Failure;
}
if ( nColors > 256 )
{
CPLError( CE_Failure, CPLE_IllegalArg,
"GDALComputeMedianCutPCT() : nColors must be lesser than or equal to 256." );
return CE_Failure;
}
if( pfnProgress == NULL )
pfnProgress = GDALDummyProgress;
/* ==================================================================== */
/* STEP 1: crate empty boxes. */
/* ==================================================================== */
int i;
Colorbox *box_list, *ptr;
int (*histogram)[GMC_B_LEN][GMC_B_LEN];
Colorbox *freeboxes;
Colorbox *usedboxes;
histogram = (int (*)[GMC_B_LEN][GMC_B_LEN])
CPLCalloc(GMC_B_LEN * GMC_B_LEN * GMC_B_LEN,sizeof(int));
usedboxes = NULL;
box_list = freeboxes = (Colorbox *)CPLMalloc(nColors*sizeof (Colorbox));
freeboxes[0].next = &freeboxes[1];
freeboxes[0].prev = NULL;
for (i = 1; i < nColors-1; ++i) {
freeboxes[i].next = &freeboxes[i+1];
freeboxes[i].prev = &freeboxes[i-1];
}
freeboxes[nColors-1].next = NULL;
freeboxes[nColors-1].prev = &freeboxes[nColors-2];
/* ==================================================================== */
/* Build histogram. */
/* ==================================================================== */
GByte *pabyRedLine, *pabyGreenLine, *pabyBlueLine;
int iLine, iPixel;
/* -------------------------------------------------------------------- */
/* Initialize the box datastructures. */
/* -------------------------------------------------------------------- */
ptr = freeboxes;
freeboxes = ptr->next;
if (freeboxes)
freeboxes->prev = NULL;
ptr->next = usedboxes;
usedboxes = ptr;
if (ptr->next)
ptr->next->prev = ptr;
ptr->rmin = ptr->gmin = ptr->bmin = 999;
ptr->rmax = ptr->gmax = ptr->bmax = -1;
ptr->total = nXSize * nYSize;
/* -------------------------------------------------------------------- */
/* Collect histogram. */
/* -------------------------------------------------------------------- */
pabyRedLine = (GByte *) VSIMalloc(nXSize);
pabyGreenLine = (GByte *) VSIMalloc(nXSize);
pabyBlueLine = (GByte *) VSIMalloc(nXSize);
if (pabyRedLine == NULL ||
pabyGreenLine == NULL ||
pabyBlueLine == NULL)
{
CPLError( CE_Failure, CPLE_OutOfMemory,
"VSIMalloc(): Out of memory in GDALComputeMedianCutPCT" );
err = CE_Failure;
goto end_and_cleanup;
}
for( iLine = 0; iLine < nYSize; iLine++ )
{
if( !pfnProgress( iLine / (double) nYSize,
"Generating Histogram", pProgressArg ) )
{
CPLError( CE_Failure, CPLE_UserInterrupt, "User Terminated" );
err = CE_Failure;
goto end_and_cleanup;
}
GDALRasterIO( hRed, GF_Read, 0, iLine, nXSize, 1,
pabyRedLine, nXSize, 1, GDT_Byte, 0, 0 );
GDALRasterIO( hGreen, GF_Read, 0, iLine, nXSize, 1,
pabyGreenLine, nXSize, 1, GDT_Byte, 0, 0 );
GDALRasterIO( hBlue, GF_Read, 0, iLine, nXSize, 1,
pabyBlueLine, nXSize, 1, GDT_Byte, 0, 0 );
for( iPixel = 0; iPixel < nXSize; iPixel++ )
{
int nRed, nGreen, nBlue;
nRed = pabyRedLine[iPixel] >> COLOR_SHIFT;
nGreen = pabyGreenLine[iPixel] >> COLOR_SHIFT;
nBlue = pabyBlueLine[iPixel] >> COLOR_SHIFT;
if(hasNoData)
{
if(nRed == pabyRedLine[iPixel] && nGreen == pabyGreenLine[iPixel] && nBlue == pabyBlueLine[iPixel])
{ // skip no data colors
continue;
}
}
ptr->rmin = MIN(ptr->rmin, nRed);
ptr->gmin = MIN(ptr->gmin, nGreen);
ptr->bmin = MIN(ptr->bmin, nBlue);
ptr->rmax = MAX(ptr->rmax, nRed);
ptr->gmax = MAX(ptr->gmax, nGreen);
ptr->bmax = MAX(ptr->bmax, nBlue);
histogram[nRed][nGreen][nBlue]++;
}
}
if( !pfnProgress( 1.0, "Generating Histogram", pProgressArg ) )
{
CPLError( CE_Failure, CPLE_UserInterrupt, "User Terminated" );
err = CE_Failure;
goto end_and_cleanup;
}
/* ==================================================================== */
/* STEP 3: continually subdivide boxes until no more free */
/* boxes remain or until all colors assigned. */
/* ==================================================================== */
while (freeboxes != NULL) {
ptr = largest_box(usedboxes);
if (ptr != NULL)
splitbox(ptr, histogram, &freeboxes, &usedboxes);
else
freeboxes = NULL;
}
/* ==================================================================== */
/* STEP 4: assign colors to all boxes */
/* ==================================================================== */
for (i = 0, ptr = usedboxes; ptr != NULL; ++i, ptr = ptr->next)
{
GDALColorEntry sEntry;
sEntry.c1 = (GByte) (((ptr->rmin + ptr->rmax) << COLOR_SHIFT) / 2);
sEntry.c2 = (GByte) (((ptr->gmin + ptr->gmax) << COLOR_SHIFT) / 2);
sEntry.c3 = (GByte) (((ptr->bmin + ptr->bmax) << COLOR_SHIFT) / 2);
sEntry.c4 = 255;
GDALSetColorEntry( hColorTable, i, &sEntry );
}
if(hasNoData)
{ // add no data entry
GDALColorEntry sEntry;
sEntry.c1 = (GByte) (iNoDataR);
sEntry.c2 = (GByte) (iNoDataG);
sEntry.c3 = (GByte) (iNoDataB);
sEntry.c4 = 0;
GDALSetColorEntry( hColorTable, nColors, &sEntry );
}
end_and_cleanup:
CPLFree( pabyRedLine );
CPLFree( pabyGreenLine );
CPLFree( pabyBlueLine );
/* We're done with the boxes now */
CPLFree(box_list);
freeboxes = usedboxes = NULL;
CPLFree( histogram );
return err;
}
/************************************************************************/
/* largest_box() */
/************************************************************************/
static Colorbox *
largest_box(Colorbox *usedboxes)
{
Colorbox *p, *b;
int size;
b = NULL;
size = -1;
for (p = usedboxes; p != NULL; p = p->next)
if ((p->rmax > p->rmin || p->gmax > p->gmin ||
p->bmax > p->bmin) && p->total > size)
size = (b = p)->total;
return (b);
}
/************************************************************************/
/* splitbox() */
/************************************************************************/
static void
splitbox(Colorbox* ptr, int (*histogram)[GMC_B_LEN][GMC_B_LEN],
Colorbox **pfreeboxes, Colorbox **pusedboxes)
{
int hist2[GMC_B_LEN];
int first=0, last=0;
Colorbox *new_cb;
int *iptr, *histp;
int i, j;
int ir,ig,ib;
int sum, sum1, sum2;
enum { RED, GREEN, BLUE } axis;
/*
* See which axis is the largest, do a histogram along that
* axis. Split at median point. Contract both new boxes to
* fit points and return
*/
i = ptr->rmax - ptr->rmin;
if (i >= ptr->gmax - ptr->gmin && i >= ptr->bmax - ptr->bmin)
axis = RED;
else if (ptr->gmax - ptr->gmin >= ptr->bmax - ptr->bmin)
axis = GREEN;
else
axis = BLUE;
/* get histogram along longest axis */
switch (axis) {
case RED:
histp = &hist2[ptr->rmin];
for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
*histp = 0;
for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
iptr = &histogram[ir][ig][ptr->bmin];
for (ib = ptr->bmin; ib <= ptr->bmax; ++ib)
*histp += *iptr++;
}
histp++;
}
first = ptr->rmin;
last = ptr->rmax;
break;
case GREEN:
histp = &hist2[ptr->gmin];
for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
*histp = 0;
for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
iptr = &histogram[ir][ig][ptr->bmin];
for (ib = ptr->bmin; ib <= ptr->bmax; ++ib)
*histp += *iptr++;
}
histp++;
}
first = ptr->gmin;
last = ptr->gmax;
break;
case BLUE:
histp = &hist2[ptr->bmin];
for (ib = ptr->bmin; ib <= ptr->bmax; ++ib) {
*histp = 0;
for (ir = ptr->rmin; ir <= ptr->rmax; ++ir) {
iptr = &histogram[ir][ptr->gmin][ib];
for (ig = ptr->gmin; ig <= ptr->gmax; ++ig) {
*histp += *iptr;
iptr += GMC_B_LEN;
}
}
histp++;
}
first = ptr->bmin;
last = ptr->bmax;
break;
}
/* find median point */
sum2 = ptr->total / 2;
histp = &hist2[first];
sum = 0;
for (i = first; i <= last && (sum += *histp++) < sum2; ++i)
;
if (i == first)
i++;
/* Create new box, re-allocate points */
new_cb = *pfreeboxes;
*pfreeboxes = new_cb->next;
if (*pfreeboxes)
(*pfreeboxes)->prev = NULL;
if (*pusedboxes)
(*pusedboxes)->prev = new_cb;
new_cb->next = *pusedboxes;
*pusedboxes = new_cb;
histp = &hist2[first];
for (sum1 = 0, j = first; j < i; j++)
sum1 += *histp++;
for (sum2 = 0, j = i; j <= last; j++)
sum2 += *histp++;
new_cb->total = sum1;
ptr->total = sum2;
new_cb->rmin = ptr->rmin;
new_cb->rmax = ptr->rmax;
new_cb->gmin = ptr->gmin;
new_cb->gmax = ptr->gmax;
new_cb->bmin = ptr->bmin;
new_cb->bmax = ptr->bmax;
switch (axis) {
case RED:
new_cb->rmax = i-1;
ptr->rmin = i;
break;
case GREEN:
new_cb->gmax = i-1;
ptr->gmin = i;
break;
case BLUE:
new_cb->bmax = i-1;
ptr->bmin = i;
break;
}
shrinkbox(new_cb, histogram);
shrinkbox(ptr, histogram);
}
/************************************************************************/
/* shrinkbox() */
/************************************************************************/
static void
shrinkbox(Colorbox* box, int (*histogram)[GMC_B_LEN][GMC_B_LEN])
{
int *histp, ir, ig, ib;
if (box->rmax > box->rmin) {
for (ir = box->rmin; ir <= box->rmax; ++ir)
for (ig = box->gmin; ig <= box->gmax; ++ig) {
histp = &histogram[ir][ig][box->bmin];
for (ib = box->bmin; ib <= box->bmax; ++ib)
if (*histp++ != 0) {
box->rmin = ir;
goto have_rmin;
}
}
have_rmin:
if (box->rmax > box->rmin)
for (ir = box->rmax; ir >= box->rmin; --ir)
for (ig = box->gmin; ig <= box->gmax; ++ig) {
histp = &histogram[ir][ig][box->bmin];
ib = box->bmin;
for (; ib <= box->bmax; ++ib)
if (*histp++ != 0) {
box->rmax = ir;
goto have_rmax;
}
}
}
have_rmax:
if (box->gmax > box->gmin) {
for (ig = box->gmin; ig <= box->gmax; ++ig)
for (ir = box->rmin; ir <= box->rmax; ++ir) {
histp = &histogram[ir][ig][box->bmin];
for (ib = box->bmin; ib <= box->bmax; ++ib)
if (*histp++ != 0) {
box->gmin = ig;
goto have_gmin;
}
}
have_gmin:
if (box->gmax > box->gmin)
for (ig = box->gmax; ig >= box->gmin; --ig)
for (ir = box->rmin; ir <= box->rmax; ++ir) {
histp = &histogram[ir][ig][box->bmin];
ib = box->bmin;
for (; ib <= box->bmax; ++ib)
if (*histp++ != 0) {
box->gmax = ig;
goto have_gmax;
}
}
}
have_gmax:
if (box->bmax > box->bmin) {
for (ib = box->bmin; ib <= box->bmax; ++ib)
for (ir = box->rmin; ir <= box->rmax; ++ir) {
histp = &histogram[ir][box->gmin][ib];
for (ig = box->gmin; ig <= box->gmax; ++ig) {
if (*histp != 0) {
box->bmin = ib;
goto have_bmin;
}
histp += GMC_B_LEN;
}
}
have_bmin:
if (box->bmax > box->bmin)
for (ib = box->bmax; ib >= box->bmin; --ib)
for (ir = box->rmin; ir <= box->rmax; ++ir) {
histp = &histogram[ir][box->gmin][ib];
ig = box->gmin;
for (; ig <= box->gmax; ++ig) {
if (*histp != 0) {
box->bmax = ib;
goto have_bmax;
}
histp += GMC_B_LEN;
}
}
}
have_bmax:
;
}
------------------------------------------------------------------------------
Full-scale, agent-less Infrastructure Monitoring from a single dashboard
Integrate with 40+ ManageEngine ITSM Solutions for complete visibility
Physical-Virtual-Cloud Infrastructure monitoring from one console
Real user monitoring with APM Insights and performance trend reports
Learn More http://pubads.g.doubleclick.net/gampad/clk?id=247754911&iu=/4140
_______________________________________________
Qlandkartegt-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/qlandkartegt-users
