1 /*
2 * tkimgMap.c
3 */
4
5 #include <string.h>
6 #include <stdlib.h>
7 #include <math.h>
8 #include "tkimg.h"
9
10 #define IMG_MIN(a,b) ((a)<(b)? (a): (b))
11 #define IMG_MAX(a,b) ((a)>(b)? (a): (b))
12 #define IMG_CLAMP(val,minVal,maxVal) (IMG_MAX (minVal, IMG_MIN (val, maxVal)))
13
14 /* The size of the histogram and lookup tables for Automatic Gain Control */
15 #define IMG_HISTOGRAM_SIZE 256
16 #define IMG_HISTOGRAM_SCALE 255.0
17
18 /* This function determines at runtime, whether we are on an Intel system. */
19
tkimg_IsIntel(void)20 int tkimg_IsIntel (void)
21 {
22 unsigned long val = 513;
23 /* On Intel (little-endian) systems this value is equal to "\01\02\00\00".
24 On big-endian systems this value equals "\00\00\02\01" */
25 return memcmp(&val, "\01\02", 2) == 0;
26 }
27
28 /*
29 *----------------------------------------------------------------------
30 *
31 * tkimg_CreateGammaTable --
32 *
33 * Returns the string representation's byte array pointer and length
34 * for an object.
35 *
36 * Results:
37 * Returns a pointer to the string representation of objPtr. If
38 * lengthPtr isn't NULL, the length of the string representation is
39 * stored at *lengthPtr. The byte array referenced by the returned
40 * pointer must not be modified by the caller. Furthermore, the
41 * caller must copy the bytes if they need to retain them since the
42 * object's string rep can change as a result of other operations.
43 * REMARK: This function reacts a little bit different than
44 * Tcl_GetStringFromObj():
45 * - objPtr is allowed to be NULL. In that case the NULL pointer
46 * will be returned, and the length will be reported to be 0;
47 * In the tkimg_ code there is never a distinction between en empty
48 * string and a NULL pointer, while the latter is easier to check
49 * for. That's the reason for this difference.
50 *
51 * Side effects:
52 * May call the object's updateStringProc to update the string
53 * representation from the internal representation.
54 *
55 *----------------------------------------------------------------------
56 */
57
tkimg_CreateGammaTable(float gammaVal,float * gammaTable)58 void tkimg_CreateGammaTable (float gammaVal, float *gammaTable)
59 {
60 int i;
61
62 for (i = 0; i < IMG_GAMMA_TABLE_SIZE - 1; i++) {
63 gammaTable[i] = pow ((float) i / (float) (IMG_GAMMA_TABLE_SIZE - 2),
64 1.0 / gammaVal);
65 }
66 gammaTable[IMG_GAMMA_TABLE_SIZE - 1] = 1.0f;
67 return;
68 }
69
70 /* Given a pixel value in Float format, "valIn", and a gamma-correction
71 * lookup table, "tab", macro "gcorrectFloat" returns the gamma-corrected
72 * pixel value in "valOut".
73 */
74
tkimg_LookupGammaTable(float val,const float * gammaTable)75 float tkimg_LookupGammaTable (float val, const float *gammaTable)
76 {
77 int gc_i;
78 float gc_t;
79 gc_t = (val) * (IMG_GAMMA_TABLE_SIZE - 2);
80 gc_i = gc_t;
81 gc_t -= gc_i;
82 return gammaTable[gc_i] * (1.0f - gc_t) + gammaTable[gc_i+1] * gc_t;
83 }
84
tkimg_UShortToUByte(int n,const unsigned short * shortIn,const float * gammaTable,unsigned char * ubOut)85 void tkimg_UShortToUByte (int n, const unsigned short *shortIn,
86 const float *gammaTable, unsigned char *ubOut)
87 {
88 const unsigned short *src, *stop;
89 unsigned char *ubdest;
90 float ftmp;
91 int itmp;
92
93 src = shortIn;
94 stop = shortIn + n;
95 ubdest = ubOut;
96
97 /* Handle a gamma value of 1.0 (gammaTable == NULL) as a special case.
98 Quite nice speed improvement for the maybe most used case. */
99 if (gammaTable) {
100 while (src < stop) {
101 ftmp = *src / 65535.0f;
102 ftmp = IMG_CLAMP (ftmp, 0.0f, 1.0f);
103 ftmp = tkimg_LookupGammaTable (ftmp, gammaTable);
104 itmp = (int)(ftmp * 255.0f + 0.5f);
105 *ubdest = IMG_CLAMP (itmp, 0, 255);
106 ++ubdest;
107 ++src;
108 }
109 } else {
110 while (src < stop) {
111 itmp = (int)(*src / 256);
112 *ubdest = IMG_CLAMP (itmp, 0, 255);
113 ++ubdest;
114 ++src;
115 }
116 }
117 return;
118 }
119
tkimg_ShortToUByte(int n,const short * shortIn,const float * gammaTable,unsigned char * ubOut)120 void tkimg_ShortToUByte (int n, const short *shortIn,
121 const float *gammaTable, unsigned char *ubOut)
122 {
123 const short *src, *stop;
124 unsigned char *ubdest;
125 float ftmp;
126 int itmp;
127
128 src = shortIn;
129 stop = shortIn + n;
130 ubdest = ubOut;
131
132 /* Handle a gamma value of 1.0 (gammaTable == NULL) as a special case.
133 Quite nice speed improvement for the maybe most used case. */
134 if (gammaTable) {
135 while (src < stop) {
136 /* Map short values from the range [MIN_SHORT .. MAX_SHORT] to
137 the range [0.0 .. 1.0], do gamma correction and then map into
138 the displayable range [0 .. 255]. */
139 ftmp = (*src * 1.0f / 65535.0f + 0.5f);
140 ftmp = tkimg_LookupGammaTable (ftmp, gammaTable);
141 itmp = (int)(ftmp * 255.0f + 0.5f);
142 *ubdest = IMG_CLAMP (itmp, 0, 255);
143 ++ubdest;
144 ++src;
145 }
146 } else {
147 while (src < stop) {
148 /* Map short values from the range [MIN_SHORT .. MAX_SHORT] to
149 the displayable range [0 .. 255]. */
150 itmp = (int)(*src * 255.0f / 65535.0f + 128);
151 *ubdest = IMG_CLAMP (itmp, 0, 255);
152 ++ubdest;
153 ++src;
154 }
155 }
156 return;
157 }
158
159 /* If no gamma correction is needed (i.e. gamma == 1.0), specify NULL for
160 * parameter gammaTable.
161 */
tkimg_FloatToUByte(int n,const float * floatIn,const float * gammaTable,unsigned char * ubOut)162 void tkimg_FloatToUByte (int n, const float *floatIn,
163 const float *gammaTable, unsigned char *ubOut)
164 {
165 const float *src, *stop;
166 unsigned char *ubdest;
167 float ftmp;
168 int itmp;
169
170 src = floatIn;
171 stop = floatIn + n;
172 ubdest = ubOut;
173
174 /* Handle a gamma value of 1.0 (gammaTable == NULL) as a special case.
175 Quite nice speed improvement for the maybe most used case. */
176 if (gammaTable) {
177 while (src < stop) {
178 ftmp = IMG_CLAMP (*src, 0.0f, 1.0f);
179 ftmp = tkimg_LookupGammaTable (ftmp, gammaTable);
180 itmp = (int)(ftmp * 255.0f + 0.5f);
181 *ubdest = IMG_CLAMP (itmp, 0, 255);
182 ++ubdest;
183 ++src;
184 }
185 } else {
186 while (src < stop) {
187 itmp = (int)(*src * 255.0f + 0.5f);
188 *ubdest = IMG_CLAMP (itmp, 0, 255);
189 ++ubdest;
190 ++src;
191 }
192 }
193 return;
194 }
195
tkimg_ReadUByteRow(tkimg_MFile * handle,unsigned char * pixels,int nBytes)196 int tkimg_ReadUByteRow (tkimg_MFile *handle, unsigned char *pixels, int nBytes)
197 {
198 #if DEBUG_READ == 1
199 printf ("Reading %d UBytes\n", nBytes); fflush (stdout);
200 #endif
201 if (nBytes != tkimg_Read2(handle, (char *)pixels, nBytes)) {
202 return 0;
203 }
204 return 1;
205 }
206
tkimg_ReadUShortRow(tkimg_MFile * handle,unsigned short * pixels,int nShorts,char * buf,int swapBytes)207 int tkimg_ReadUShortRow (tkimg_MFile *handle, unsigned short *pixels, int nShorts,
208 char *buf, int swapBytes)
209 {
210 int i;
211 char *bufPtr = buf;
212 unsigned short *mPtr = pixels;
213
214 #if DEBUG_READ == 1
215 printf ("Reading %d UShorts\n", nShorts); fflush (stdout);
216 #endif
217 if (2 * nShorts != tkimg_Read2(handle, buf, 2 * nShorts)) {
218 return 0;
219 }
220
221 if (swapBytes) {
222 for (i=0; i<nShorts; i++) {
223 ((char *)mPtr)[0] = bufPtr[1];
224 ((char *)mPtr)[1] = bufPtr[0];
225 mPtr++;
226 bufPtr += 2;
227 }
228 } else {
229 for (i=0; i<nShorts; i++) {
230 ((char *)mPtr)[0] = bufPtr[0];
231 ((char *)mPtr)[1] = bufPtr[1];
232 mPtr++;
233 bufPtr += 2;
234 }
235 }
236 return 1;
237 }
238
tkimg_ReadShortRow(tkimg_MFile * handle,short * pixels,int nShorts,char * buf,int swapBytes)239 int tkimg_ReadShortRow (tkimg_MFile *handle, short *pixels, int nShorts,
240 char *buf, int swapBytes)
241 {
242 int i;
243 char *bufPtr = buf;
244 short *mPtr = pixels;
245
246 #if DEBUG_READ == 1
247 printf ("Reading %d Shorts\n", nShorts); fflush (stdout);
248 #endif
249 if (2 * nShorts != tkimg_Read2(handle, buf, 2 * nShorts)) {
250 return 0;
251 }
252
253 if (swapBytes) {
254 for (i=0; i<nShorts; i++) {
255 ((char *)mPtr)[0] = bufPtr[1];
256 ((char *)mPtr)[1] = bufPtr[0];
257 mPtr++;
258 bufPtr += 2;
259 }
260 } else {
261 for (i=0; i<nShorts; i++) {
262 ((char *)mPtr)[0] = bufPtr[0];
263 ((char *)mPtr)[1] = bufPtr[1];
264 mPtr++;
265 bufPtr += 2;
266 }
267 }
268 return 1;
269 }
270
tkimg_ReadFloatRow(tkimg_MFile * handle,float * pixels,int nFloats,char * buf,int swapBytes)271 int tkimg_ReadFloatRow (tkimg_MFile *handle, float *pixels, int nFloats,
272 char *buf, int swapBytes)
273 {
274 int i;
275 char *bufPtr = buf;
276 float *mPtr = pixels;
277
278 #if DEBUG_READ == 1
279 printf ("Reading %d floats\n", nFloats); fflush (stdout);
280 #endif
281 if (4 * nFloats != tkimg_Read2(handle, buf, 4 * nFloats)) {
282 return 0;
283 }
284
285 if (swapBytes) {
286 for (i=0; i<nFloats; i++) {
287 ((char *)mPtr)[0] = bufPtr[3];
288 ((char *)mPtr)[1] = bufPtr[2];
289 ((char *)mPtr)[2] = bufPtr[1];
290 ((char *)mPtr)[3] = bufPtr[0];
291 mPtr++;
292 bufPtr += 4;
293 }
294 } else {
295 for (i=0; i<nFloats; i++) {
296 ((char *)mPtr)[0] = bufPtr[0];
297 ((char *)mPtr)[1] = bufPtr[1];
298 ((char *)mPtr)[2] = bufPtr[2];
299 ((char *)mPtr)[3] = bufPtr[3];
300 mPtr++;
301 bufPtr += 4;
302 }
303 }
304 return 1;
305 }
306
tkimg_ReadUByteFile(tkimg_MFile * handle,unsigned char * buf,int width,int height,int nchan,int verbose,int findMinMax,float * minVals,float * maxVals)307 int tkimg_ReadUByteFile (tkimg_MFile *handle, unsigned char *buf, int width, int height,
308 int nchan, int verbose, int findMinMax,
309 float *minVals, float *maxVals)
310 {
311 int x, y, c;
312 unsigned char *bufPtr = buf;
313
314 #if DEBUG_READ == 1
315 printf ("readUByteFile: Width=%d Height=%d nchan=%d findMinMax=%s\n",
316 width, height, nchan, findMinMax? "yes": "no");
317 fflush (stdout);
318 #endif
319 for (c=0; c<nchan; c++) {
320 minVals[c] = (float) 1.0E30;
321 maxVals[c] = (float)-1.0E30;
322 }
323 for (y=0; y<height; y++) {
324 if (! tkimg_ReadUByteRow (handle, bufPtr, nchan * width)) {
325 return 0;
326 }
327 if (findMinMax) {
328 for (x=0; x<width; x++) {
329 for (c=0; c<nchan; c++) {
330 if (*bufPtr > maxVals[c]) maxVals[c] = *bufPtr;
331 if (*bufPtr < minVals[c]) minVals[c] = *bufPtr;
332 bufPtr++;
333 }
334 }
335 } else {
336 bufPtr += nchan * width;
337 }
338 }
339 if (verbose && findMinMax) {
340 printf ("\tMinimum pixel values :");
341 for (c=0; c<nchan; c++) {
342 printf (" %d", (unsigned char)minVals[c]);
343 }
344 printf ("\n");
345 printf ("\tMaximum pixel values :");
346 for (c=0; c<nchan; c++) {
347 printf (" %d", (unsigned char)maxVals[c]);
348 }
349 printf ("\n");
350 fflush (stdout);
351 }
352 return 1;
353 }
354
tkimg_ReadUShortFile(tkimg_MFile * handle,unsigned short * buf,int width,int height,int nchan,int swapBytes,int verbose,int findMinMax,float * minVals,float * maxVals)355 int tkimg_ReadUShortFile (tkimg_MFile *handle, unsigned short *buf, int width, int height,
356 int nchan, int swapBytes, int verbose, int findMinMax,
357 float *minVals, float *maxVals)
358 {
359 int x, y, c;
360 char *line;
361 unsigned short *bufPtr = buf;
362
363 #if DEBUG_READ == 1
364 printf ("readUShortFile: Width=%d Height=%d nchan=%d swapBytes=%s findMinMax=%s\n",
365 width, height, nchan, swapBytes? "yes": "no", findMinMax? "yes": "no"); fflush (stdout);
366 #endif
367 for (c=0; c<nchan; c++) {
368 minVals[c] = (float) 1.0E30;
369 maxVals[c] = (float)-1.0E30;
370 }
371 line = ckalloc (sizeof (unsigned short) * nchan * width);
372
373 for (y=0; y<height; y++) {
374 if (! tkimg_ReadUShortRow (handle, bufPtr, nchan * width, line, swapBytes)) {
375 return 0;
376 }
377 if (findMinMax) {
378 for (x=0; x<width; x++) {
379 for (c=0; c<nchan; c++) {
380 if (*bufPtr > maxVals[c]) maxVals[c] = *bufPtr;
381 if (*bufPtr < minVals[c]) minVals[c] = *bufPtr;
382 bufPtr++;
383 }
384 }
385 } else {
386 bufPtr += nchan * width;
387 }
388 }
389 if (verbose && findMinMax) {
390 printf ("\tMinimum pixel values :");
391 for (c=0; c<nchan; c++) {
392 printf (" %d", (unsigned short)minVals[c]);
393 }
394 printf ("\n");
395 printf ("\tMaximum pixel values :");
396 for (c=0; c<nchan; c++) {
397 printf (" %d", (unsigned short)maxVals[c]);
398 }
399 printf ("\n");
400 fflush (stdout);
401 }
402 ckfree (line);
403 return 1;
404 }
405
tkimg_ReadFloatFile(tkimg_MFile * handle,float * buf,int width,int height,int nchan,int swapBytes,int verbose,int findMinMax,float * minVals,float * maxVals,float saturation)406 int tkimg_ReadFloatFile (tkimg_MFile *handle, float *buf, int width, int height,
407 int nchan, int swapBytes, int verbose, int findMinMax,
408 float *minVals, float *maxVals, float saturation)
409 {
410 int x, y, c;
411 float value;
412 char *line;
413 float *bufPtr = buf;
414
415 #if DEBUG_READ == 1
416 printf ("readFloatFile: Width=%d Height=%d nchan=%d swapBytes=%s findMinMax=%s\n",
417 width, height, nchan, swapBytes? "yes": "no", findMinMax? "yes": "no"); fflush (stdout);
418 #endif
419 if (saturation <= 0.0f) {
420 saturation = (float) 1.0E30;
421 }
422 for (c=0; c<nchan; c++) {
423 minVals[c] = (float) 1.0E30;
424 maxVals[c] = (float)-1.0E30;
425 }
426 line = ckalloc (sizeof (float) * nchan * width);
427
428 for (y=0; y<height; y++) {
429 if (! tkimg_ReadFloatRow (handle, bufPtr, nchan * width, line, swapBytes)) {
430 return 0;
431 }
432 if (findMinMax) {
433 for (x=0; x<width; x++) {
434 for (c=0; c<nchan; c++) {
435 value = IMG_MIN (*bufPtr, saturation);
436 if (value > maxVals[c]) maxVals[c] = value;
437 if (value < minVals[c]) minVals[c] = value;
438 bufPtr++;
439 }
440 }
441 } else {
442 bufPtr += nchan * width;
443 }
444 }
445 if (verbose && findMinMax) {
446 printf ("\tMinimum pixel values :");
447 for (c=0; c<nchan; c++) {
448 printf (" %f", minVals[c]);
449 }
450 printf ("\n");
451 printf ("\tMaximum pixel values :");
452 for (c=0; c<nchan; c++) {
453 printf (" %f", maxVals[c]);
454 }
455 printf ("\n");
456 fflush (stdout);
457 }
458 ckfree (line);
459 return 1;
460 }
461
tkimg_RemapUShortValues(unsigned short * buf,int width,int height,int nchan,float * minVals,float * maxVals)462 void tkimg_RemapUShortValues (unsigned short *buf, int width, int height, int nchan,
463 float *minVals, float *maxVals)
464 {
465 int x, y, c;
466 float m[IMG_MAX_CHANNELS], t[IMG_MAX_CHANNELS];
467 unsigned int scaledVal;
468 unsigned short *bufPtr = buf;
469
470 for (c=0; c<nchan; c++) {
471 m[c] = (65535.0 - 0.0) / (maxVals[c] - minVals[c]);
472 t[c] = 0.0 - m[c] * minVals[c];
473 }
474 for (y=0; y<height; y++) {
475 for (x=0; x<width; x++) {
476 for (c=0; c<nchan; c++) {
477 scaledVal = *bufPtr * m[c] + t[c];
478 *bufPtr = IMG_CLAMP (scaledVal, 0, 65535);
479 bufPtr++;
480 }
481 }
482 }
483 return;
484 }
485
getHistogramIndex(float val,float minVal,float maxVal)486 static int getHistogramIndex (float val, float minVal, float maxVal )
487 {
488 float scaledVal;
489 int histoInd;
490
491 scaledVal = IMG_HISTOGRAM_SCALE * (IMG_MAX (val - minVal, 0.0) / (maxVal - minVal));
492 histoInd = (int) IMG_CLAMP (scaledVal, 0.0, IMG_HISTOGRAM_SCALE);
493 return histoInd;
494 }
495
histogramFloat(float * buf,int width,int height,int nchan,float minVals[],float maxVals[],int histogram[],int printAgc)496 static void histogramFloat (float *buf, int width, int height, int nchan,
497 float minVals[], float maxVals[], int histogram[],
498 int printAgc)
499 {
500 int x, y, c;
501 int histoInd;
502 float *bufPtr = buf;
503
504 memset (histogram, 0, IMG_HISTOGRAM_SIZE * sizeof (int));
505 /* Currently supported only for 1 channel images. */
506 c = 0;
507 for (y=0; y<height; y++) {
508 for (x=0; x<width; x++) {
509 histoInd = getHistogramIndex (*bufPtr, minVals[c], maxVals[c]);
510 histogram[histoInd]++;
511 bufPtr++;
512 }
513 }
514 if (printAgc) {
515 int i;
516 int count = 0;
517 printf("agc globalMin %f\n", minVals[0]);
518 printf("agc globalMax %f\n", maxVals[0]);
519 for(i=0; i<IMG_HISTOGRAM_SIZE; i++) {
520 printf ("agc histogram %3d %5d\n", i, histogram[i]);
521 if (histogram[i] != 0) count++;
522 }
523 printf( "agc histostat %d %d\n", count, IMG_HISTOGRAM_SIZE - count);
524 }
525 return;
526 }
527
tkimg_RemapFloatValues(float * buf,int width,int height,int nchan,float * minVals,float * maxVals,float agcCutOffPercent,int printAgc)528 void tkimg_RemapFloatValues (float *buf, int width, int height, int nchan,
529 float *minVals, float *maxVals,
530 float agcCutOffPercent, int printAgc)
531 {
532 int x, y, c;
533 float *bufPtr = buf;
534 float m[IMG_MAX_CHANNELS], t[IMG_MAX_CHANNELS];
535 float scaledVal;
536 float minNewVals[IMG_MAX_CHANNELS], maxNewVals[IMG_MAX_CHANNELS];
537
538 for (c=0; c<nchan; c++) {
539 minNewVals[c] = minVals[c];
540 maxNewVals[c] = maxVals[c];
541 }
542
543 if (agcCutOffPercent > 0.0) {
544 int i;
545 int minLutInd = -1, maxLutInd = -1;
546 int histogram[IMG_HISTOGRAM_SIZE];
547 float lut[IMG_HISTOGRAM_SIZE];
548 float agcCutOff = IMG_CLAMP(agcCutOffPercent * 0.01, 0.0, 1.0);
549 float sum = 0.0;
550
551 /* Calculate histogram. */
552 histogramFloat (buf, width, height, nchan, minVals, maxVals,
553 histogram, printAgc);
554
555 /* Accumulate the histogram and divide by the image size. */
556 for (i=0; i<IMG_HISTOGRAM_SIZE; i++) {
557 sum = sum + histogram[i];
558 lut[i] = sum / (double)(width * height);
559 if( lut[i] >= agcCutOff && minLutInd < 0 ) {
560 minLutInd = i;
561 }
562 if (lut[i] >= (1.0 - agcCutOff) && maxLutInd < 0) {
563 maxLutInd = i;
564 }
565 if (printAgc) {
566 printf ("agc lut %3d %.3f\n", i, lut[i]);
567 }
568 }
569
570 for (c=0; c<nchan; c++) {
571 minNewVals[c] = minLutInd * (maxVals[c] - minVals[c]) /
572 IMG_HISTOGRAM_SCALE + minVals[c];
573 maxNewVals[c] = maxLutInd * (maxVals[c] - minVals[c]) /
574 IMG_HISTOGRAM_SCALE + minVals[c];
575 if (printAgc) {
576 printf ("agc cutOff %f\n", agcCutOff);
577 printf ("agc lutMinInd %d\n", minLutInd);
578 printf ("agc lutMaxInd %d\n", maxLutInd);
579 printf ("agc lutMin %f\n", minNewVals[c]);
580 printf ("agc lutMax %f\n", maxNewVals[c]);
581 }
582 }
583 }
584
585 for (c=0; c<nchan; c++) {
586 m[c] = (1.0 - 0.0) / (maxNewVals[c] - minNewVals[c]);
587 t[c] = 0.0 - m[c] * minNewVals[c];
588 }
589
590 for (y=0; y<height; y++) {
591 for (x=0; x<width; x++) {
592 for (c=0; c<nchan; c++) {
593 scaledVal = *bufPtr * m[c] + t[c];
594 *bufPtr = IMG_CLAMP (scaledVal, 0.0, 1.0);
595 bufPtr++;
596 }
597 }
598 }
599 return;
600 }
601
602
603