1 : /******************************************************************************
2 : * $Id: gdaldither.cpp 13342 2007-12-14 20:58:31Z rouault $
3 : *
4 : * Project: CIETMap Phase 2
5 : * Purpose: Convert RGB (24bit) to a pseudo-colored approximation using
6 : * Floyd-Steinberg dithering (error diffusion).
7 : * Author: Frank Warmerdam, warmerdam@pobox.com
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 2001, Frank Warmerdam
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : ******************************************************************************
30 : *
31 : * Notes:
32 : *
33 : * [1] Floyd-Steinberg dither:
34 : * I should point out that the actual fractions we used were, assuming
35 : * you are at X, moving left to right:
36 : *
37 : * X 7/16
38 : * 3/16 5/16 1/16
39 : *
40 : * Note that the error goes to four neighbors, not three. I think this
41 : * will probably do better (at least for black and white) than the
42 : * 3/8-3/8-1/4 distribution, at the cost of greater processing. I have
43 : * seen the 3/8-3/8-1/4 distribution described as "our" algorithm before,
44 : * but I have no idea who the credit really belongs to.
45 : * --
46 : * Lou Steinberg
47 : *
48 : */
49 :
50 : #include "gdal_priv.h"
51 : #include "gdal_alg.h"
52 :
53 : CPL_CVSID("$Id: gdaldither.cpp 13342 2007-12-14 20:58:31Z rouault $");
54 :
55 : #define C_LEVELS 32
56 :
57 : static void FindNearestColor( int nColors, int *panPCT, GByte *pabyColorMap );
58 :
59 : /************************************************************************/
60 : /* GDALDitherRGB2PCT() */
61 : /************************************************************************/
62 :
63 : /**
64 : * 24bit to 8bit conversion with dithering.
65 : *
66 : * This functions utilizes Floyd-Steinberg dithering in the process of
67 : * converting a 24bit RGB image into a pseudocolored 8bit image using a
68 : * provided color table.
69 : *
70 : * The red, green and blue input bands do not necessarily need to come
71 : * from the same file, but they must be the same width and height. They will
72 : * be clipped to 8bit during reading, so non-eight bit bands are generally
73 : * inappropriate. Likewise the hTarget band will be written with 8bit values
74 : * and must match the width and height of the source bands.
75 : *
76 : * The color table cannot have more than 256 entries.
77 : *
78 : * @param hRed Red input band.
79 : * @param hGreen Green input band.
80 : * @param hBlue Blue input band.
81 : * @param hTarget Output band.
82 : * @param hColorTable the color table to use with the output band.
83 : * @param pfnProgress callback for reporting algorithm progress matching the
84 : * GDALProgressFunc() semantics. May be NULL.
85 : * @param pProgressArg callback argument passed to pfnProgress.
86 : *
87 : * @return CE_None on success or CE_Failure if an error occurs.
88 : */
89 :
90 : int CPL_STDCALL
91 1 : GDALDitherRGB2PCT( GDALRasterBandH hRed,
92 : GDALRasterBandH hGreen,
93 : GDALRasterBandH hBlue,
94 : GDALRasterBandH hTarget,
95 : GDALColorTableH hColorTable,
96 : GDALProgressFunc pfnProgress,
97 : void * pProgressArg )
98 :
99 : {
100 1 : VALIDATE_POINTER1( hRed, "GDALDitherRGB2PCT", CE_Failure );
101 1 : VALIDATE_POINTER1( hGreen, "GDALDitherRGB2PCT", CE_Failure );
102 1 : VALIDATE_POINTER1( hBlue, "GDALDitherRGB2PCT", CE_Failure );
103 1 : VALIDATE_POINTER1( hTarget, "GDALDitherRGB2PCT", CE_Failure );
104 1 : VALIDATE_POINTER1( hColorTable, "GDALDitherRGB2PCT", CE_Failure );
105 :
106 : int nXSize, nYSize;
107 1 : CPLErr err = CE_None;
108 :
109 : /* -------------------------------------------------------------------- */
110 : /* Validate parameters. */
111 : /* -------------------------------------------------------------------- */
112 1 : nXSize = GDALGetRasterBandXSize( hRed );
113 1 : nYSize = GDALGetRasterBandYSize( hRed );
114 :
115 1 : if( GDALGetRasterBandXSize( hGreen ) != nXSize
116 : || GDALGetRasterBandYSize( hGreen ) != nYSize
117 : || GDALGetRasterBandXSize( hBlue ) != nXSize
118 : || GDALGetRasterBandYSize( hBlue ) != nYSize )
119 : {
120 : CPLError( CE_Failure, CPLE_IllegalArg,
121 0 : "Green or blue band doesn't match size of red band.\n" );
122 :
123 0 : return CE_Failure;
124 : }
125 :
126 1 : if( GDALGetRasterBandXSize( hTarget ) != nXSize
127 : || GDALGetRasterBandYSize( hTarget ) != nYSize )
128 : {
129 : CPLError( CE_Failure, CPLE_IllegalArg,
130 : "GDALDitherRGB2PCT(): "
131 0 : "Target band doesn't match size of source bands.\n" );
132 :
133 0 : return CE_Failure;
134 : }
135 :
136 1 : if( pfnProgress == NULL )
137 1 : pfnProgress = GDALDummyProgress;
138 :
139 : /* -------------------------------------------------------------------- */
140 : /* Setup more direct colormap. */
141 : /* -------------------------------------------------------------------- */
142 : int nColors, anPCT[768], iColor;
143 :
144 1 : nColors = GDALGetColorEntryCount( hColorTable );
145 :
146 1 : if (nColors == 0 )
147 : {
148 : CPLError( CE_Failure, CPLE_IllegalArg,
149 : "GDALDitherRGB2PCT(): "
150 0 : "Color table must not be empty.\n" );
151 :
152 0 : return CE_Failure;
153 : }
154 1 : else if (nColors > 256)
155 : {
156 : CPLError( CE_Failure, CPLE_IllegalArg,
157 : "GDALDitherRGB2PCT(): "
158 0 : "Color table cannot have more than 256 entries.\n" );
159 :
160 0 : return CE_Failure;
161 : }
162 :
163 9 : for( iColor = 0; iColor < nColors; iColor++ )
164 : {
165 : GDALColorEntry sEntry;
166 :
167 8 : GDALGetColorEntryAsRGB( hColorTable, iColor, &sEntry );
168 :
169 8 : anPCT[iColor ] = sEntry.c1;
170 8 : anPCT[iColor+256] = sEntry.c2;
171 8 : anPCT[iColor+512] = sEntry.c3;
172 : }
173 :
174 : /* -------------------------------------------------------------------- */
175 : /* Build a 24bit to 8 bit color mapping. */
176 : /* -------------------------------------------------------------------- */
177 : GByte *pabyColorMap;
178 :
179 : pabyColorMap = (GByte *) CPLMalloc(C_LEVELS * C_LEVELS * C_LEVELS
180 1 : * sizeof(int));
181 :
182 1 : FindNearestColor( nColors, anPCT, pabyColorMap );
183 :
184 : /* -------------------------------------------------------------------- */
185 : /* Setup various variables. */
186 : /* -------------------------------------------------------------------- */
187 : GByte *pabyRed, *pabyGreen, *pabyBlue, *pabyIndex;
188 : int *panError;
189 :
190 1 : pabyRed = (GByte *) VSIMalloc(nXSize);
191 1 : pabyGreen = (GByte *) VSIMalloc(nXSize);
192 1 : pabyBlue = (GByte *) VSIMalloc(nXSize);
193 :
194 1 : pabyIndex = (GByte *) VSIMalloc(nXSize);
195 :
196 1 : panError = (int *) VSICalloc(sizeof(int),(nXSize+2) * 3);
197 :
198 1 : if (pabyRed == NULL ||
199 : pabyGreen == NULL ||
200 : pabyBlue == NULL ||
201 : pabyIndex == NULL ||
202 : panError == NULL)
203 : {
204 : CPLError( CE_Failure, CPLE_OutOfMemory,
205 0 : "VSIMalloc(): Out of memory in GDALDitherRGB2PCT" );
206 0 : err = CE_Failure;
207 0 : goto end_and_cleanup;
208 : }
209 :
210 : /* ==================================================================== */
211 : /* Loop over all scanlines of data to process. */
212 : /* ==================================================================== */
213 : int iScanline;
214 :
215 51 : for( iScanline = 0; iScanline < nYSize; iScanline++ )
216 : {
217 : int nLastRedError, nLastGreenError, nLastBlueError, i;
218 :
219 : /* -------------------------------------------------------------------- */
220 : /* Report progress */
221 : /* -------------------------------------------------------------------- */
222 50 : if( !pfnProgress( iScanline / (double) nYSize, NULL, pProgressArg ) )
223 : {
224 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User Terminated" );
225 0 : err = CE_Failure;
226 0 : goto end_and_cleanup;
227 : }
228 :
229 : /* -------------------------------------------------------------------- */
230 : /* Read source data. */
231 : /* -------------------------------------------------------------------- */
232 : GDALRasterIO( hRed, GF_Read, 0, iScanline, nXSize, 1,
233 50 : pabyRed, nXSize, 1, GDT_Byte, 0, 0 );
234 : GDALRasterIO( hGreen, GF_Read, 0, iScanline, nXSize, 1,
235 50 : pabyGreen, nXSize, 1, GDT_Byte, 0, 0 );
236 : GDALRasterIO( hBlue, GF_Read, 0, iScanline, nXSize, 1,
237 50 : pabyBlue, nXSize, 1, GDT_Byte, 0, 0 );
238 :
239 : /* -------------------------------------------------------------------- */
240 : /* Apply the error from the previous line to this one. */
241 : /* -------------------------------------------------------------------- */
242 2550 : for( i = 0; i < nXSize; i++ )
243 : {
244 2500 : pabyRed[i] = (GByte)
245 2500 : MAX(0,MIN(255,(pabyRed[i] + panError[i*3+0+3])));
246 2500 : pabyGreen[i] = (GByte)
247 2500 : MAX(0,MIN(255,(pabyGreen[i] + panError[i*3+1+3])));
248 2500 : pabyBlue[i] = (GByte)
249 2500 : MAX(0,MIN(255,(pabyBlue[i] + panError[i*3+2+3])));
250 : }
251 :
252 50 : memset( panError, 0, sizeof(int) * (nXSize+2) * 3 );
253 :
254 : /* -------------------------------------------------------------------- */
255 : /* Figure out the nearest color to the RGB value. */
256 : /* -------------------------------------------------------------------- */
257 50 : nLastRedError = 0;
258 50 : nLastGreenError = 0;
259 50 : nLastBlueError = 0;
260 :
261 2550 : for( i = 0; i < nXSize; i++ )
262 : {
263 : int iIndex, nError, nSixth, iRed, iGreen, iBlue;
264 : int nRedValue, nGreenValue, nBlueValue;
265 :
266 2500 : nRedValue = MAX(0,MIN(255, pabyRed[i] + nLastRedError));
267 2500 : nGreenValue = MAX(0,MIN(255, pabyGreen[i] + nLastGreenError));
268 2500 : nBlueValue = MAX(0,MIN(255, pabyBlue[i] + nLastBlueError));
269 :
270 2500 : iRed = nRedValue * C_LEVELS / 256;
271 2500 : iGreen = nGreenValue * C_LEVELS / 256;
272 2500 : iBlue = nBlueValue * C_LEVELS / 256;
273 :
274 : iIndex = pabyColorMap[iRed + iGreen * C_LEVELS
275 2500 : + iBlue * C_LEVELS * C_LEVELS];
276 :
277 2500 : pabyIndex[i] = (GByte) iIndex;
278 :
279 : /* -------------------------------------------------------------------- */
280 : /* Compute Red error, and carry it on to the next error line. */
281 : /* -------------------------------------------------------------------- */
282 2500 : nError = nRedValue - anPCT[iIndex ];
283 2500 : nSixth = nError / 6;
284 :
285 2500 : panError[i*3 ] += nSixth;
286 2500 : panError[i*3+6 ] = nSixth;
287 2500 : panError[i*3+3 ] += nError - 5 * nSixth;
288 :
289 2500 : nLastRedError = 2 * nSixth;
290 :
291 : /* -------------------------------------------------------------------- */
292 : /* Compute Green error, and carry it on to the next error line. */
293 : /* -------------------------------------------------------------------- */
294 2500 : nError = nGreenValue - anPCT[iIndex+256];
295 2500 : nSixth = nError / 6;
296 :
297 2500 : panError[i*3 +1] += nSixth;
298 2500 : panError[i*3+6+1] = nSixth;
299 2500 : panError[i*3+3+1] += nError - 5 * nSixth;
300 :
301 2500 : nLastGreenError = 2 * nSixth;
302 :
303 : /* -------------------------------------------------------------------- */
304 : /* Compute Blue error, and carry it on to the next error line. */
305 : /* -------------------------------------------------------------------- */
306 2500 : nError = nBlueValue - anPCT[iIndex+512];
307 2500 : nSixth = nError / 6;
308 :
309 2500 : panError[i*3 +2] += nSixth;
310 2500 : panError[i*3+6+2] = nSixth;
311 2500 : panError[i*3+3+2] += nError - 5 * nSixth;
312 :
313 2500 : nLastBlueError = 2 * nSixth;
314 : }
315 :
316 : /* -------------------------------------------------------------------- */
317 : /* Write results. */
318 : /* -------------------------------------------------------------------- */
319 : GDALRasterIO( hTarget, GF_Write, 0, iScanline, nXSize, 1,
320 50 : pabyIndex, nXSize, 1, GDT_Byte, 0, 0 );
321 : }
322 :
323 1 : pfnProgress( 1.0, NULL, pProgressArg );
324 :
325 : /* -------------------------------------------------------------------- */
326 : /* Cleanup */
327 : /* -------------------------------------------------------------------- */
328 : end_and_cleanup:
329 1 : CPLFree( pabyRed );
330 1 : CPLFree( pabyGreen );
331 1 : CPLFree( pabyBlue );
332 1 : CPLFree( pabyIndex );
333 1 : CPLFree( panError );
334 1 : CPLFree( pabyColorMap );
335 :
336 1 : return err;
337 : }
338 :
339 : /************************************************************************/
340 : /* FindNearestColor() */
341 : /* */
342 : /* Finear near PCT color for any RGB color. */
343 : /************************************************************************/
344 :
345 1 : static void FindNearestColor( int nColors, int *panPCT, GByte *pabyColorMap )
346 :
347 : {
348 : int iBlue, iGreen, iRed;
349 : int iColor;
350 :
351 : /* -------------------------------------------------------------------- */
352 : /* Loop over all the cells in the high density cube. */
353 : /* -------------------------------------------------------------------- */
354 33 : for( iBlue = 0; iBlue < C_LEVELS; iBlue++ )
355 : {
356 1056 : for( iGreen = 0; iGreen < C_LEVELS; iGreen++ )
357 : {
358 33792 : for( iRed = 0; iRed < C_LEVELS; iRed++ )
359 : {
360 : int nRedValue, nGreenValue, nBlueValue;
361 32768 : int nBestDist = 768, nBestIndex = 0;
362 :
363 32768 : nRedValue = (iRed * 255) / (C_LEVELS-1);
364 32768 : nGreenValue = (iGreen * 255) / (C_LEVELS-1);
365 32768 : nBlueValue = (iBlue * 255) / (C_LEVELS-1);
366 :
367 294912 : for( iColor = 0; iColor < nColors; iColor++ )
368 : {
369 : int nThisDist;
370 :
371 262144 : nThisDist = ABS(nRedValue - panPCT[iColor ])
372 262144 : + ABS(nGreenValue - panPCT[iColor+256])
373 524288 : + ABS(nBlueValue - panPCT[iColor+512]);
374 :
375 262144 : if( nThisDist < nBestDist )
376 : {
377 102835 : nBestIndex = iColor;
378 102835 : nBestDist = nThisDist;
379 : }
380 : }
381 :
382 : pabyColorMap[iRed + iGreen*C_LEVELS
383 32768 : + iBlue*C_LEVELS*C_LEVELS] = (GByte)nBestIndex;
384 : }
385 : }
386 : }
387 1 : }
388 :
|
