src/gdcmSegmentedPalette.h

   1 /*=========================================================================
   2
   3   Program:   gdcm
   4   Module:    $RCSfile: gdcmSegmentedPalette.h,v $
   5   Language:  C++
   6   Date:      $Date: 2007/11/09 17:53:19 $
   7   Version:   $Revision: 1.19 $
   8
   9   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
  10   l'Image). All rights reserved. See Doc/License.txt or
  11   http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
  12
  13      This software is distributed WITHOUT ANY WARRANTY; without even
  14      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15      PURPOSE.  See the above copyright notices for more information.
  16
  17 =========================================================================*/
  18
  19 #ifndef _GDCMSEGMENTEDPALETTE_H_
  20 #define _GDCMSEGMENTEDPALETTE_H_
  21
  22 #include "gdcmFile.h"
  23 #include "gdcmTagKey.h"
  24 #include "gdcmDataEntry.h"
  25
  26 // Notepad from Satomi Takeo:
  27 // http://blog.goo.ne.jp/satomi_takeo/e/3643e5249b2a9650f9e10ef1c830e8b8
  28 // I bet the code was compiled on VS6. Make it compile on other platform:
  29 // * typedef are not inherited
  30 // * need to explicitely add typename keyword
  31 // * Uint8 / Uint16 are neither C nor C++
  32 // * replace all dcmtk code with equivalent gdcm code
  33 // * Add extra parameter to ReadPaletteInto to save extracted info
  34
  35 // MM: Over 50000 DICOM files I only found two that had a Segmented Palette LUT
  36 // And both were coming from a ALOKA SSD-4000 station
  37 //
  38 // JPRx :
  39 // It compiles fine under gcc 4.1.2, msvc 7 and DarwinG5-g++
  40 // It doesn't compile on msvc 6, borland and SunOS
  41 // Any fix welcome !
  42
  43 #include <assert.h>
  44 #include <algorithm>
  45 #include <deque>
  46 #include <map>
  47 #include <vector>
  48 #include <iterator>
  49
  50 // Hack for VS6
  51 #if defined(_MSC_VER) && (_MSC_VER < 1310)
  52 #define GDCM_TYPENAME
  53 #else
  54 #define GDCM_TYPENAME typename
  55 #endif
  56
  57 // Hack for Borland
  58 #if (defined(_MSC_VER) && (_MSC_VER < 1310)) || defined(__BORLANDC__)
  59 #define GDCM_TYPENAME2
  60 #else
  61 #define GDCM_TYPENAME2 typename
  62 #endif
  63
  64
  65 namespace GDCM_NAME_SPACE
  66 {
  67 // Long stody short: Sun compiler only provide the second interface of std::distance, since the implemenation
  68 // is so trivial, I'd rather redo it myself.
  69 // Ref:
  70 // http://www.sgi.com/tech/stl/distance.html#2
  71 // The second version of distance was the one defined in the original STL, and the first version is the one defined in the draft C++ standard; the definition was changed because the older interface was clumsy and error-prone. The older interface required the use of a temporary variable, and it has semantics that are somewhat nonintuitive: it increments n by the distance from first to last, rather than storing that distance in n
  72   template<typename InputIterator>
  73     inline int32_t
  74     mydistance(InputIterator first, InputIterator last)
  75       {
  76       int32_t n = 0;
  77       while (first != last)
  78         {
  79         ++first;
  80         ++n;
  81         }
  82       return n;
  83       }
  84
  85     // abstract class for segment.
  86     template <typename EntryType>
  87     class Segment {
  88     public:
  89         typedef std::map<const EntryType* const, const Segment*> SegmentMap;
  90         virtual bool Expand(const SegmentMap& instances,
  91             std::vector<EntryType>& expanded) const = 0;
  92         const EntryType* First() const { return _first; }
  93         const EntryType* Last()  const { return _last;  }
  94         struct ToMap {
  95             typename SegmentMap::value_type
  96                 operator()(const Segment* segment) const
  97             { return std::make_pair(segment->First(), segment); }
  98         };
  99     protected:
 100         Segment(const EntryType* first, const EntryType* last) {
 101             _first = first; _last = last;
 102         }
 103         const EntryType* _first;
 104         const EntryType* _last;
 105     };
 106
 107     // discrete segment (opcode = 0)
 108     template <typename EntryType>
 109     class DiscreteSegment : public Segment<EntryType> {
 110     public:
 111 #if !defined(__sun__)
 112         typedef typename Segment<EntryType>::SegmentMap SegmentMap;
 113 #endif
 114         DiscreteSegment(const EntryType* first)
 115             : Segment<EntryType>(first, first+2+*(first+1)) {}
 116         virtual bool Expand(const SegmentMap&,
 117             std::vector<EntryType>& expanded) const
 118         {
 119             std::copy(this->_first + 2, this->_last, std::back_inserter(expanded));
 120             return true;
 121         }
 122     };
 123
 124     // linear segment (opcode = 1)
 125     template <typename EntryType>
 126     class LinearSegment : public Segment<EntryType> {
 127     public:
 128 #if !defined(__sun__)
 129         typedef typename Segment<EntryType>::SegmentMap SegmentMap;
 130 #endif
 131         LinearSegment(const EntryType* first)
 132             : Segment<EntryType>(first, first+3) {}
 133         virtual bool Expand(const SegmentMap&,
 134             std::vector<EntryType>& expanded) const
 135         {
 136             if ( expanded.empty() ) {
 137                 // linear segment can't be the first segment.
 138                 return false;
 139             }
 140             EntryType length = *(this->_first + 1);
 141             EntryType y0 = expanded.back();
 142             EntryType y1 = *(this->_first + 2);
 143             double y01 = y1 - y0;
 144             for ( EntryType i = 0; i <length; ++i ) {
 145                 double value_float
 146                     = static_cast<double>(y0)
 147                     + (static_cast<double>(i)/static_cast<double>(length)) * y01;
 148                 EntryType value_int = static_cast<EntryType>(value_float + 0.5);
 149                 expanded.push_back(value_int);
 150             }
 151             return true;
 152         }
 153     };
 154
 155     // indirect segment (opcode = 2)
 156     template <typename EntryType>
 157     class IndirectSegment : public Segment<EntryType> {
 158     public:
 159 #if !defined(__sun__)
 160         typedef typename Segment<EntryType>::SegmentMap SegmentMap;
 161 #endif
 162         IndirectSegment(const EntryType* first)
 163             : Segment<EntryType>(first, first+2+4/sizeof(EntryType)) {}
 164         virtual bool Expand(const SegmentMap& instances,
 165             std::vector<EntryType>& expanded) const
 166         {
 167             if ( instances.empty() ) {
 168                 // some other segments are required as references.
 169                 return false;
 170             }
 171             const EntryType* first_segment = instances.begin()->first;
 172             const unsigned short* pOffset
 173                 = reinterpret_cast<const unsigned short*>(this->_first + 2);
 174             unsigned long offsetBytes
 175                 = (*pOffset) | (static_cast<unsigned long>(*(pOffset + 1)) << 16);
 176             const EntryType* copied_part_head
 177                 = first_segment + offsetBytes / sizeof(EntryType);
 178             GDCM_TYPENAME SegmentMap::const_iterator ppHeadSeg = instances.find(copied_part_head);
 179             if ( ppHeadSeg == instances.end() ) {
 180                 // referred segment not found
 181                 return false;
 182             }
 183             EntryType nNumCopies = *(this->_first + 1);
 184             GDCM_TYPENAME SegmentMap::const_iterator ppSeg = ppHeadSeg;
 185             while ( mydistance(ppHeadSeg, ppSeg) < nNumCopies ) {
 186                 // assert( mydistance(ppHeadSeg, ppSeg) == std::distance(ppHeadSeg, ppSeg) );
 187                 assert( ppSeg != instances.end() );
 188                 ppSeg->second->Expand(instances, expanded);
 189                 ++ppSeg;
 190             }
 191             return true;
 192         }
 193     };
 194
 195     template <typename EntryType>
 196     void ExpandPalette(const EntryType* raw_values, uint32_t length,
 197         std::vector<EntryType>& palette)
 198     {
 199         typedef std::deque<Segment<EntryType>*> SegmentList;
 200         SegmentList segments;
 201         const EntryType* raw_seg = raw_values;
 202         while ( (mydistance(raw_values, raw_seg) * sizeof(EntryType)) < length ) {
 203             // assert( mydistance(raw_values, raw_seg) == std::distance(raw_values, raw_seg) );
 204             Segment<EntryType>* segment = NULL;
 205             if ( *raw_seg == 0 ) {
 206                 segment = new DiscreteSegment<EntryType>(raw_seg);
 207             } else if ( *raw_seg == 1 ) {
 208                 segment = new LinearSegment<EntryType>(raw_seg);
 209             } else if ( *raw_seg == 2 ) {
 210                 segment = new IndirectSegment<EntryType>(raw_seg);
 211             }
 212             if ( segment ) {
 213                 segments.push_back(segment);
 214                 raw_seg = segment->Last();
 215             } else {
 216                 // invalid opcode
 217                 break;
 218             }
 219
 220         }
 221         GDCM_TYPENAME Segment<EntryType>::SegmentMap instances;
 222         std::transform(segments.begin(), segments.end(),
 223             std::inserter(instances, instances.end()), GDCM_TYPENAME2 Segment<EntryType>::ToMap());
 224         GDCM_TYPENAME SegmentList::iterator ppSeg = segments.begin();
 225         GDCM_TYPENAME SegmentList::iterator endOfSegments = segments.end();
 226         for ( ; ppSeg != endOfSegments; ++ppSeg ) {
 227             (*ppSeg)->Expand(instances, palette);
 228         }
 229         ppSeg = segments.begin();
 230         for ( ; ppSeg != endOfSegments; ++ppSeg ) {
 231             delete *ppSeg;
 232         }
 233     }
 234
 235     void ReadPaletteInto(GDCM_NAME_SPACE::File* pds, const GDCM_NAME_SPACE::TagKey& descriptor,
 236       const GDCM_NAME_SPACE::TagKey& segment, uint8_t* lut)
 237       {
 238       unsigned int desc_values[3] = {0,0,0};
 239       unsigned long count = 0;
 240       //if ( pds->findAndGetUint16Array(descriptor, desc_values, &count).good() )
 241       std::string desc_values_str = pds->GetEntryString(descriptor.GetGroup(), descriptor.GetElement() );
 242       count = sscanf( desc_values_str.c_str(), "%u\\%u\\%u", desc_values, desc_values+1, desc_values+2 );
 243         {
 244         assert( count == 3 );
 245         unsigned int num_entries = desc_values[0];
 246         if ( num_entries == 0 ) {
 247           num_entries = 0x10000;
 248         }
 249         unsigned int min_pixel_value = desc_values[1];
 250         assert( min_pixel_value == 0 ); // FIXME
 251         unsigned int entry_size = desc_values[2];
 252         assert( entry_size == 8 || entry_size == 16 );
 253         //DcmElement* pe = NULL;
 254         GDCM_NAME_SPACE::DataEntry* pe = NULL;
 255
 256         pe = pds->GetDataEntry(segment.GetGroup(), segment.GetElement() );
 257           {
 258           //if ( pds->findAndGetElement(segment, pe).good() )
 259           unsigned long length = pe->GetLength();
 260           if ( entry_size == 8 )
 261             {
 262             uint8_t* segment_values = NULL;
 263             //if ( pe->getUint8Array(segment_values).good() )
 264             segment_values = (uint8_t*)pe->GetBinArea();
 265               {
 266               std::vector<uint8_t> palette;
 267               palette.reserve(num_entries);
 268               ExpandPalette(segment_values, length, palette);
 269               memcpy(lut, &palette[0], palette.size() );
 270               }
 271             }
 272           else if ( entry_size == 16 )
 273             {
 274             uint16_t* segment_values = NULL;
 275             segment_values = (uint16_t*)pe->GetBinArea();
 276               {
 277               //if ( pe->getUint16Array(segment_values).good() )
 278               std::vector<uint16_t> palette;
 279               palette.reserve(num_entries);
 280               ExpandPalette(segment_values, length, palette);
 281               memcpy(lut, &palette[0], palette.size()*2 );
 282
 283               }
 284             }
 285           }
 286         }
 287       }
 288 } // end namespace gdcm
 289
 290
 291 // do not pollute namespace:
 292 #undef GDCM_TYPENAME
 293
 294 #endif