1 /*
2  *
3  *  Copyright (C) 2002-2016, OFFIS e.V.
4  *  All rights reserved.  See COPYRIGHT file for details.
5  *
6  *  This software and supporting documentation were developed by
7  *
8  *    OFFIS e.V.
9  *    R&D Division Health
10  *    Escherweg 2
11  *    D-26121 Oldenburg, Germany
12  *
13  *
14  *  Module:  dcmimage
15  *
16  *  Author:  Marco Eichelberg
17  *
18  *  Purpose: class DcmQuantColorTable
19  *
20  */
21 
22 #include "dcmtk/config/osconfig.h"
23 
24 #include "dcmtk/dcmimage/diqtctab.h"
25 #include "dcmtk/dcmimage/diqtpbox.h"
26 #include "dcmtk/dcmdata/dcerror.h"   /* for EC_IllegalCall */
27 #include "dcmtk/dcmdata/dcelem.h"    /* for DcmElement */
28 #include "dcmtk/dcmdata/dcitem.h"    /* for DcmItem */
29 #include "dcmtk/dcmdata/dcvrus.h"    /* for DcmUnsignedShort */
30 #include "dcmtk/dcmdata/dcvrobow.h"  /* for DcmOtherByteOtherWord */
31 #include "dcmtk/dcmdata/dcswap.h"    /* for swapIfNecessary() */
32 #include "dcmtk/dcmdata/dcdeftag.h"  /* for tag constants */
33 #include "dcmtk/dcmdata/dcuid.h"     /* for OFFIS_DCMTK_VERSION */
34 
35 #define INCLUDE_CSTDIO
36 #include "dcmtk/ofstd/ofstdinc.h"
37 
38 /* ------------------------------------------------------------ */
39 
40 // static comparison functions for qsort
41 
42 BEGIN_EXTERN_C
redcompare(const void * x1,const void * x2)43 static int redcompare(const void *x1, const void *x2)
44 {
45   return OFstatic_cast(int, (*OFstatic_cast(const DcmQuantHistogramItemPointer *, x1))->getRed())
46        - OFstatic_cast(int, (*OFstatic_cast(const DcmQuantHistogramItemPointer *, x2))->getRed());
47 }
48 
greencompare(const void * x1,const void * x2)49 static int greencompare(const void *x1, const void *x2)
50 {
51   return OFstatic_cast(int, (*OFstatic_cast(const DcmQuantHistogramItemPointer *, x1))->getGreen())
52        - OFstatic_cast(int, (*OFstatic_cast(const DcmQuantHistogramItemPointer *, x2))->getGreen());
53 }
54 
bluecompare(const void * x1,const void * x2)55 static int bluecompare(const void *x1, const void *x2)
56 {
57   return OFstatic_cast(int, (*OFstatic_cast(const DcmQuantHistogramItemPointer *, x1))->getBlue())
58        - OFstatic_cast(int, (*OFstatic_cast(const DcmQuantHistogramItemPointer *, x2))->getBlue());
59 }
60 END_EXTERN_C
61 
62 /* ------------------------------------------------------------ */
63 
DcmQuantColorTable()64 DcmQuantColorTable::DcmQuantColorTable()
65 : array(NULL)
66 , numColors(0)
67 , maxval(0)
68 {
69 }
70 
71 
~DcmQuantColorTable()72 DcmQuantColorTable::~DcmQuantColorTable()
73 {
74   clear();
75 }
76 
77 
clear()78 void DcmQuantColorTable::clear()
79 {
80   if (array)
81   {
82     for (unsigned long i=0; i < numColors; i++) delete array[i];
83     delete[] array;
84     array = NULL;
85   }
86   numColors = 0;
87   maxval = 0;
88 }
89 
90 
computeHistogram(DicomImage & image,unsigned long maxcolors)91 OFCondition DcmQuantColorTable::computeHistogram(
92   DicomImage& image,
93   unsigned long maxcolors)
94 {
95   // reset object to initial state
96   clear();
97 
98   // compute initial maxval
99   maxval = OFstatic_cast(DcmQuantComponent, -1);
100   DcmQuantColorHashTable *htable = NULL;
101 
102   // attempt to make a histogram of the colors, unclustered.
103   // If at first we don't succeed, lower maxval to increase color
104   // coherence and try again.  This will eventually terminate.
105   OFBool done = OFFalse;
106   while (! done)
107   {
108     htable = new DcmQuantColorHashTable();
109     numColors = htable->addToHashTable(image, maxval, maxcolors);
110     if (numColors > 0) done = OFTrue;
111     else
112     {
113       delete htable;
114       maxval = maxval/2;
115     }
116   }
117 
118   numColors = htable->createHistogram(array);
119   delete htable;
120   return EC_Normal;
121 }
122 
123 
medianCut(DcmQuantColorTable & histogram,unsigned long sum,unsigned long theMaxval,unsigned long numberOfColors,DcmLargestDimensionType largeType,DcmRepresentativeColorType repType)124 OFCondition DcmQuantColorTable::medianCut(
125   DcmQuantColorTable& histogram,
126   unsigned long sum,
127   unsigned long theMaxval,
128   unsigned long numberOfColors,
129   DcmLargestDimensionType largeType,
130   DcmRepresentativeColorType repType)
131 {
132   // reset object to initial state
133   clear();
134 
135   // update maxval
136   maxval = theMaxval;
137 
138   // initialize color table
139   array = new DcmQuantHistogramItemPointer[numberOfColors];
140   for (unsigned int xx=0; xx < numberOfColors; xx++) array[xx] = new DcmQuantHistogramItem();
141   numColors = numberOfColors;
142 
143   int i;
144   unsigned int bi;
145   DcmQuantPixelBoxArray bv(numberOfColors);
146 
147   // Set up the initial box.
148   bv[0].ind = 0;
149   bv[0].colors = OFstatic_cast(int, histogram.getColors());
150   bv[0].sum = sum;
151   unsigned int boxes = 1;
152 
153   // Main loop: split boxes until we have enough.
154   while ( boxes < numberOfColors )
155   {
156       int indx, clrs;
157       unsigned long sm;
158       int minr, maxr, ming, maxg, minb, maxb, v;
159       unsigned long halfsum, lowersum;
160 
161       // Find the first splittable box.
162       for ( bi = 0; bi < boxes; ++bi )
163           if ( bv[bi].colors >= 2 )
164           break;
165       if ( bi == boxes )
166           break;  /* ran out of colors! */
167       indx = bv[bi].ind;
168       clrs = bv[bi].colors;
169       sm = bv[bi].sum;
170 
171       // Go through the box finding the minimum and maximum of each
172       // component - the boundaries of the box.
173       minr = maxr = histogram.array[indx]->getRed();
174       ming = maxg = histogram.array[indx]->getGreen();
175       minb = maxb = histogram.array[indx]->getBlue();
176       for ( i = 1; i < clrs; ++i )
177       {
178           v = histogram.array[indx+i]->getRed();
179           if ( v < minr ) minr = v;
180           if ( v > maxr ) maxr = v;
181           v = histogram.array[indx+i]->getGreen();
182           if ( v < ming ) ming = v;
183           if ( v > maxg ) maxg = v;
184           v = histogram.array[indx+i]->getBlue();
185           if ( v < minb ) minb = v;
186           if ( v > maxb ) maxb = v;
187       }
188 
189       /*
190       ** Find the largest dimension, and sort by that component.  I have
191       ** included two methods for determining the "largest" dimension;
192       ** first by simply comparing the range in RGB space, and second
193       ** by transforming into luminosities before the comparison.
194       */
195       if (largeType == DcmLargestDimensionType_default)
196       {
197           if ( maxr - minr >= maxg - ming && maxr - minr >= maxb - minb )
198               qsort(OFreinterpret_cast(char*, &(histogram.array[indx])), clrs, sizeof(DcmQuantHistogramItemPointer), redcompare);
199           else if ( maxg - ming >= maxb - minb )
200               qsort(OFreinterpret_cast(char*, &(histogram.array[indx])), clrs, sizeof(DcmQuantHistogramItemPointer), greencompare);
201           else
202               qsort(OFreinterpret_cast(char*, &(histogram.array[indx])), clrs, sizeof(DcmQuantHistogramItemPointer), bluecompare);
203        }
204        else // DcmLargestDimensionType_luminance
205        {
206           double rl, gl, bl;
207           DcmQuantPixel p;
208 
209           p.assign(maxr - minr, 0, 0);
210           rl = p.luminance();
211 
212           p.assign(0, maxg - ming, 0);
213           gl = p.luminance();
214 
215           p.assign(0, 0, maxb - minb);
216           bl = p.luminance();
217 
218           if ( rl >= gl && rl >= bl )
219               qsort(OFreinterpret_cast(char*, &(histogram.array[indx])), clrs, sizeof(DcmQuantHistogramItemPointer), redcompare);
220           else if ( gl >= bl )
221               qsort(OFreinterpret_cast(char*, &(histogram.array[indx])), clrs, sizeof(DcmQuantHistogramItemPointer), greencompare);
222           else
223               qsort(OFreinterpret_cast(char*, &(histogram.array[indx])), clrs, sizeof(DcmQuantHistogramItemPointer), bluecompare);
224       }
225 
226       /*
227       ** Now find the median based on the counts, so that about half the
228       ** pixels (not colors, pixels) are in each subdivision.
229       */
230       lowersum = histogram.array[indx]->getValue();
231       halfsum = sm / 2;
232       for ( i = 1; i < clrs - 1; ++i )
233       {
234           if ( lowersum >= halfsum )
235             break;
236           lowersum += histogram.array[indx+i]->getValue();
237       }
238 
239       // Split the box, and sort to bring the biggest boxes to the top.
240       bv[bi].colors = i;
241       bv[bi].sum = lowersum;
242       bv[boxes].ind = indx + i;
243       bv[boxes].colors = clrs - i;
244       bv[boxes].sum = sm - lowersum;
245       ++boxes;
246       bv.sort(boxes);
247   }
248 
249   /*
250   ** Ok, we've got enough boxes.  Now choose a representative color for
251   ** each box.  There are a number of possible ways to make this choice.
252   ** One would be to choose the center of the box; this ignores any structure
253   ** within the boxes.  Another method would be to average all the colors in
254   ** the box - this is the method specified in Heckbert's paper.  A third
255   ** method is to average all the pixels in the box.
256   */
257   if (repType == DcmRepresentativeColorType_centerOfBox)
258   {
259       for ( bi = 0; bi < boxes; ++bi )
260       {
261           int indx = bv[bi].ind;
262           int clrs = bv[bi].colors;
263           int minr, maxr, ming, maxg, minb, maxb, v;
264 
265           minr = maxr = histogram.array[indx]->getRed();
266           ming = maxg = histogram.array[indx]->getGreen();
267           minb = maxb = histogram.array[indx]->getBlue();
268           for ( i = 1; i < clrs; ++i )
269           {
270               v = histogram.array[indx+i]->getRed();
271               minr = (minr < v ? minr : v);
272               maxr = (minr > v ? minr : v);
273               v = histogram.array[indx+i]->getGreen();
274               ming = (ming < v ? ming : v);
275               maxg = (ming > v ? ming : v);
276               v = histogram.array[indx+i]->getBlue();
277               minb = (minb < v ? minb : v);
278               maxb = (minb > v ? minb : v);
279           }
280           array[bi]->assign(( minr + maxr ) / 2, ( ming + maxg ) / 2, ( minb + maxb ) / 2);
281       }
282   }
283   else if (repType == DcmRepresentativeColorType_default)
284   {
285       for ( bi = 0; bi < boxes; ++bi )
286       {
287           int indx = bv[bi].ind;
288           int clrs = bv[bi].colors;
289           long r = 0, g = 0, b = 0;
290 
291           for ( i = 0; i < clrs; ++i )
292           {
293               r += histogram.array[indx+i]->getRed();
294               g += histogram.array[indx+i]->getGreen();
295               b += histogram.array[indx+i]->getBlue();
296           }
297           r = r / clrs;
298           g = g / clrs;
299           b = b / clrs;
300           array[bi]->assign(OFstatic_cast(DcmQuantComponent, r), OFstatic_cast(DcmQuantComponent, g), OFstatic_cast(DcmQuantComponent, b));
301       }
302   }
303   else // DcmRepresentativeColorType_averagePixels
304   {
305       for ( bi = 0; bi < boxes; ++bi )
306       {
307           int indx = bv[bi].ind;
308           int clrs = bv[bi].colors;
309           unsigned long r = 0, g = 0, b = 0, sumVal = 0;
310 
311           for ( i = 0; i < clrs; ++i )
312           {
313               r += histogram.array[indx+i]->getRed()   * histogram.array[indx+i]->getValue();
314               g += histogram.array[indx+i]->getGreen() * histogram.array[indx+i]->getValue();
315               b += histogram.array[indx+i]->getBlue()  * histogram.array[indx+i]->getValue();
316               sumVal += histogram.array[indx+i]->getValue();
317           }
318           r = r / sumVal;
319           if ( r > maxval ) r = maxval;   /* avoid math errors */
320           g = g / sumVal;
321           if ( g > maxval ) g = maxval;
322           b = b / sumVal;
323           if ( b > maxval ) b = maxval;
324           array[bi]->assign(OFstatic_cast(DcmQuantComponent, r), OFstatic_cast(DcmQuantComponent, g), OFstatic_cast(DcmQuantComponent, b));
325       }
326   }
327 
328   // All done, now compute clusters
329   computeClusters();
330   return EC_Normal;
331 }
332 
333 
computeClusters()334 void DcmQuantColorTable::computeClusters()
335 {
336   unsigned long i;
337   unsigned long j;
338   unsigned long k=0;
339   int cluster;
340   int newdist;
341   int r1, g1, b1;
342   int r2, g2, b2;
343 
344   // initialize clusters
345   for (i = 0; i < numColors; ++i) array[i]->setValue(2000000000);
346 
347   for (i = 0; i < numColors-1; ++i)
348   {
349     cluster = array[i]->getValue();
350     r1 = OFstatic_cast(int, array[i]->getRed());
351     g1 = OFstatic_cast(int, array[i]->getGreen());
352     b1 = OFstatic_cast(int, array[i]->getBlue());
353 
354     for (j = i+1; j < numColors; ++j)
355     {
356       // compute euclidean distance between i and j
357       r2 = r1 - OFstatic_cast(int, array[j]->getRed());
358       g2 = g1 - OFstatic_cast(int, array[j]->getGreen());
359       b2 = b1 - OFstatic_cast(int, array[j]->getBlue());
360       newdist = (r2*r2 + g2*g2 + b2*b2)/2;
361       if (newdist < cluster)
362       {
363         cluster = newdist;
364         k = j;
365       }
366     }
367     array[i]->setValue(cluster);
368     array[k]->setValue(cluster);
369   }
370 }
371 
372 
write(DcmItem & target,OFBool writeAsOW,OFBool write16BitEntries)373 OFCondition DcmQuantColorTable::write(
374   DcmItem& target,
375   OFBool writeAsOW,
376   OFBool write16BitEntries)
377 {
378   if (numColors == 0) return EC_IllegalCall;
379 
380   // if we're using 16 bit per sample anyway, force 16 bit palette entries
381   if (sizeof(DcmQuantComponent) > 1) write16BitEntries = OFTrue;
382 
383   OFCondition result = EC_Normal;
384   if (array)
385   {
386     // create palette color LUT descriptor
387     Uint16 descriptor[3];
388     descriptor[0] = (numColors > 65535) ? 0 : OFstatic_cast(Uint16, numColors); // number of entries
389     descriptor[1] = 0; // first pixel value mapped
390     descriptor[2] = write16BitEntries ? 16 : 8; // bits per entry, must be 8 or 16.
391 
392     // if we're writing a 16-bit LUT with 64K entries, we must write as OW because
393     // otherwise the US length field will overflow when writing in explicit VR!
394     if ((descriptor[0] == 0) && write16BitEntries) writeAsOW = OFTrue;
395 
396     DcmElement *elem;
397     if (result.good())
398     {
399       elem = new DcmUnsignedShort(DCM_RedPaletteColorLookupTableDescriptor);
400       if (elem)
401       {
402         result = elem->putUint16Array(descriptor, 3);
403         if (result.good()) result = target.insert(elem, OFTrue);
404         if (result.bad()) delete elem;
405       } else result = EC_MemoryExhausted;
406     }
407 
408     if (result.good())
409     {
410       elem = new DcmUnsignedShort(DCM_GreenPaletteColorLookupTableDescriptor);
411       if (elem)
412       {
413         result = elem->putUint16Array(descriptor, 3);
414         if (result.good()) result = target.insert(elem, OFTrue);
415         if (result.bad()) delete elem;
416       } else result = EC_MemoryExhausted;
417     }
418 
419     if (result.good())
420     {
421       elem = new DcmUnsignedShort(DCM_BluePaletteColorLookupTableDescriptor);
422       if (elem)
423       {
424         result = elem->putUint16Array(descriptor, 3);
425         if (result.good()) result = target.insert(elem, OFTrue);
426         if (result.bad()) delete elem;
427       } else result = EC_MemoryExhausted;
428     }
429 
430     // now create the LUT content
431     if (result.good())
432     {
433       Uint16* rLUT = NULL;
434       Uint16* gLUT = NULL;
435       Uint16* bLUT = NULL;
436       unsigned long numWords = 0;
437       double factor = 1.0;
438       if (write16BitEntries)
439       {
440         numWords = numColors;
441         rLUT = new Uint16[numWords];
442         gLUT = new Uint16[numWords];
443         bLUT = new Uint16[numWords];
444         factor = 65535.0 / maxval;
445         if (rLUT && gLUT && bLUT)
446         {
447           for (unsigned long i=0; i<numColors; i++)
448           {
449             rLUT[i] = OFstatic_cast(Uint16, OFstatic_cast(double, array[i]->getRed()) * factor);
450             gLUT[i] = OFstatic_cast(Uint16, OFstatic_cast(double, array[i]->getGreen()) * factor);
451             bLUT[i] = OFstatic_cast(Uint16, OFstatic_cast(double, array[i]->getBlue()) * factor);
452 
453             // if the source data is only 8 bits per entry, replicate high and low bytes
454             if (sizeof(DcmQuantComponent) == 1)
455             {
456               rLUT[i] |= rLUT[i] << 8;
457               gLUT[i] |= gLUT[i] << 8;
458               bLUT[i] |= bLUT[i] << 8;
459             }
460           }
461         } else result = EC_MemoryExhausted;
462       }
463       else
464       {
465         // number of Uint16 words needed to store numColors Uint8 values plus padding
466         numWords = (numColors+1)/2;
467         rLUT = new Uint16[numWords];
468         gLUT = new Uint16[numWords];
469         bLUT = new Uint16[numWords];
470         rLUT[numWords-1] = 0;
471         gLUT[numWords-1] = 0;
472         bLUT[numWords-1] = 0;
473         factor = 255.0 / maxval;
474         if (rLUT && gLUT && bLUT)
475         {
476           Uint8 *rLUT8 = OFreinterpret_cast(Uint8 *, rLUT);
477           Uint8 *gLUT8 = OFreinterpret_cast(Uint8 *, gLUT);
478           Uint8 *bLUT8 = OFreinterpret_cast(Uint8 *, bLUT);
479           for (unsigned long i=0; i<numColors; i++)
480           {
481             rLUT8[i] = OFstatic_cast(Uint8, OFstatic_cast(double, array[i]->getRed()) * factor);
482             gLUT8[i] = OFstatic_cast(Uint8, OFstatic_cast(double, array[i]->getGreen()) * factor);
483             bLUT8[i] = OFstatic_cast(Uint8, OFstatic_cast(double, array[i]->getBlue()) * factor);
484           }
485           // we have written the byte array in little endian order, now swap to US/OW if neccessary
486           swapIfNecessary(gLocalByteOrder, EBO_LittleEndian, rLUT, numWords*sizeof(Uint16), sizeof(Uint16));
487           swapIfNecessary(gLocalByteOrder, EBO_LittleEndian, gLUT, numWords*sizeof(Uint16), sizeof(Uint16));
488           swapIfNecessary(gLocalByteOrder, EBO_LittleEndian, bLUT, numWords*sizeof(Uint16), sizeof(Uint16));
489         } else result = EC_MemoryExhausted;
490       }
491 
492       // the LUT data is prepared, now create the corresponding DICOM elements
493       if (result.good())
494       {
495         if (writeAsOW) elem = new DcmOtherByteOtherWord(DcmTag(DCM_RedPaletteColorLookupTableData, EVR_OW));
496         else elem = new DcmUnsignedShort(DcmTag(DCM_RedPaletteColorLookupTableData, EVR_US));
497         if (elem)
498         {
499           result = elem->putUint16Array(rLUT, numWords);
500           if (result.good()) result = target.insert(elem, OFTrue);
501           if (result.bad()) delete elem;
502         } else result = EC_MemoryExhausted;
503       }
504 
505       if (result.good())
506       {
507         if (writeAsOW) elem = new DcmOtherByteOtherWord(DcmTag(DCM_GreenPaletteColorLookupTableData, EVR_OW));
508         else elem = new DcmUnsignedShort(DcmTag(DCM_GreenPaletteColorLookupTableData, EVR_US));
509         if (elem)
510         {
511           result = elem->putUint16Array(gLUT, numWords);
512           if (result.good()) result = target.insert(elem, OFTrue);
513           if (result.bad()) delete elem;
514         } else result = EC_MemoryExhausted;
515       }
516 
517       if (result.good())
518       {
519         if (writeAsOW) elem = new DcmOtherByteOtherWord(DcmTag(DCM_BluePaletteColorLookupTableData, EVR_OW));
520         else elem = new DcmUnsignedShort(DcmTag(DCM_BluePaletteColorLookupTableData, EVR_US));
521         if (elem)
522         {
523           result = elem->putUint16Array(bLUT, numWords);
524           if (result.good()) result = target.insert(elem, OFTrue);
525           if (result.bad()) delete elem;
526         } else result = EC_MemoryExhausted;
527       }
528       delete[] rLUT;
529       delete[] gLUT;
530       delete[] bLUT;
531     }
532   } else result = EC_IllegalCall;
533 
534   return result;
535 }
536 
setDescriptionString(OFString & str) const537 void DcmQuantColorTable::setDescriptionString(OFString& str) const
538 {
539   char buf[100];
540 
541   sprintf(buf, "Converted to PALETTE COLOR %lu/0/%u with DCMTK %s",
542     (numColors > 65535) ? 0 : numColors, (sizeof(DcmQuantComponent) == 1) ? 8 : 16,
543     OFFIS_DCMTK_VERSION);
544 
545   str = buf;
546 }
547