VMTK
importDicom.cpp
1 /*
2  * importDicom.cpp: Ler um arquivo no formato DICOM, extrair a
3  * imagem 3D e reorientá-la para que fique positivamente
4  * orientada no referencial LPS
5  *
6  * Copyright (C) 2013 Wu Shin-Ting, FEEC, Unicamp
7  *
8  * This program is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program. If not, see <http://www.gnu.org/licenses/>.
20  *
21  */
22 #include <sys/types.h>
23 #include <sys/stat.h>
24 #include <cstdio>
25 #include <cstring>
26 #include <string>
27 #include <iostream>
28 #include <stdio.h> /* for printf */
29 #include <stdint.h>
30 #include <stdlib.h> /* for exit */
31 #include <string.h>
32 #include <limits.h>
33 #include <math.h>
34 #ifndef WIN32
35 #include <getopt.h>
36 #endif
37 
38 /*============================================================
39  Incluir definicoes da bilbioteca gdcm
40 
41 http://www.creatis.insa-lyon.fr/software/public/Gdcm/
42  ===========================================================*/
43 #include "gdcmReader.h"
44 #include "gdcmImageReader.h"
45 #include "gdcmMediaStorage.h"
46 #include "gdcmFile.h"
47 #include "gdcmDataSet.h"
48 #include "gdcmUIDs.h"
49 #include "gdcmGlobal.h"
50 #include "gdcmModules.h"
51 #include "gdcmDefs.h"
52 #include "gdcmOrientation.h"
53 #include "gdcmVersion.h"
54 #include "gdcmMD5.h"
55 #include "gdcmSystem.h"
56 #include "gdcmDirectory.h"
57 #include "gdcmSorter.h"
58 #include "gdcmScanner.h"
59 #include "gdcmDataSet.h"
60 #include "gdcmAttribute.h"
61 
62 // Estrutura de dados para armazenar a informacao do volume de dados
63 #include "importDicom.h"
64 
65 //http://nipy.sourceforge.net/nibabel/dicom/dicom_orientation.html
66 //Reorienta o volume de dados para que seja LPS
67 
68 bool Import::LPS_ReorientImageVolume(ImgFormat *volData)
69 {
70  unsigned int ix, iy, iz, idx, idy, idz;
71  unsigned int dimx, dimy, dimz;
72  float sx, sy, sz;
73  unsigned char flag;
74  unsigned char code;
75 
76  unsigned char *tmp = new unsigned char[2*volData->dims[0]*volData->dims[1]];
77 
78  std::cout << "Reorient the samples in LPS-reference ..." << std::endl;
79 
80  code = 0;
81 
82  // Reflection
83  // X-axis
84  flag = 0;
85  if (fabs(volData->dircos[0][0]) > fabs(volData->dircos[0][1]) &&
86  fabs(volData->dircos[0][0]) > fabs(volData->dircos[0][2])) {
87  if (volData->dircos[0][0] < 0) { // RL (21/05/2013 - Ting)
88  flag = 1;
89  volData->dircos[0][0] *= (-1); // LR
90  }
91  code += 0;
92  } else if (fabs(volData->dircos[0][1]) > fabs(volData->dircos[0][0])
93  && fabs(volData->dircos[0][1]) > fabs(volData->dircos[0][2])) {
94  if (volData->dircos[0][1] < 0) { // PA
95  flag = 1;
96  volData->dircos[0][1] *= (-1); // AP
97  }
98  code += 16;
99  } else {
100  if (volData->dircos[0][2] < 0) { // HF
101  flag = 1;
102  volData->dircos[0][2] *= (-1); // FH
103  }
104  code += 32;
105  }
106  if (flag) {
107  switch (volData->nbitsalloc) {
108  case 8:
109  for (idz = 0; idz < volData->dims[2]; idz++) {
110  for (idy = 0; idy < volData->dims[1]; idy++) {
111  for (ix=0,idx = volData->dims[0]-1; ix < idx; idx--,ix++) {
113  tmp[0]=volData->buffer[
114  idz * volData->dims[0] * volData->dims[1] +
115  idy * volData->dims[0] + idx];
116  volData->buffer[idz * volData->dims[0] * volData->dims[1] +
117  idy * volData->dims[0] + idx] =
118  volData->buffer[idz * volData->dims[0] * volData->dims[1] +
119  idy * volData->dims[0] + ix];
120  volData->buffer[idz * volData->dims[0] * volData->dims[1] +
121  idy * volData->dims[0] + ix] = tmp[0];
122  }
123  }
124  }
125  break;
126  case 16:
127  for (idz = 0; idz < volData->dims[2]; idz++) {
128  for (idy = 0; idy < volData->dims[1]; idy++) {
129  for (ix=0,idx = volData->dims[0]-1; ix < idx; idx=idx--,ix=ix++) {
131  memcpy (tmp,
132  &(volData->buffer[2*(idz * volData->dims[0] *
133  volData->dims[1] + idy * volData->dims[0] + idx)]),
134  sizeof(unsigned char)*2);
135  memcpy (&(volData->buffer[
136  2*(idz * volData->dims[0] * volData->dims[1] +
137  idy * volData->dims[0] + idx)]),
138  &(volData->buffer[
139  2*(idz * volData->dims[0] * volData->dims[1] +
140  idy * volData->dims[0] + ix)]),
141  sizeof(unsigned char)*2);
142  memcpy(&(volData->buffer[
143  2*(idz * volData->dims[0] * volData->dims[1] +
144  idy * volData->dims[1] + ix)]), tmp,
145  sizeof(unsigned char)*2);
146  }
147  }
148  }
149  break;
150  }
151  }
152 
154  flag = 0;
155  if (fabs(volData->dircos[1][0]) > fabs(volData->dircos[1][1])
156  &&
157  fabs(volData->dircos[1][0]) > fabs(volData->dircos[1][2])) {
158  if (volData->dircos[1][0] < 0) { // RL (21/05/2013 - Ting)
159  flag = 1;
160  volData->dircos[1][0] *= (-1); // LR
161  }
162  code += 0;
163  } else if (fabs(volData->dircos[1][1]) > fabs(volData->dircos[1][0])
164  &&
165  fabs(volData->dircos[1][1]) > fabs(volData->dircos[1][2])) {
166  if (volData->dircos[1][1] < 0) { // PA
167  flag = 1;
168  volData->dircos[1][1] *= (-1); // AP
169  }
170  code += 4;
171  } else {
172  if (volData->dircos[1][2] < 0) { // HF
173  flag = 1;
174  volData->dircos[1][2] *= (-1); // FH
175  }
176  code += 8;
177  }
178  if (flag) {
179  switch (volData->nbitsalloc) {
180  case 8:
181  for (idz = 0; idz < volData->dims[2]; idz++) {
182  for (iy=0,idy = volData->dims[1]-1; iy < idy; idy--,iy++) {
183  memcpy (tmp, &(volData->buffer[idz * volData->dims[0] *
184  volData->dims[1] + idy * volData->dims[0]]),
185  sizeof(unsigned char)*volData->dims[0]);
186  memcpy(&(volData->buffer[idz * volData->dims[0] *
187  volData->dims[1] + idy * volData->dims[0]]),
188  &(volData->buffer[idz * volData->dims[0] *
189  volData->dims[1] + iy * volData->dims[0]]),
190  sizeof(unsigned char)*volData->dims[0]);
191  memcpy(&(volData->buffer[idz * volData->dims[0] *
192  volData->dims[1] + iy * volData->dims[0]]), tmp,
193  sizeof(unsigned char)*volData->dims[0]);
194  }
195  }
196  break;
197  case 16:
198  for (idz = 0; idz < volData->dims[2]; idz++) {
199  for (iy=0,idy = volData->dims[1]-1; iy < idy; idy--,iy++) {
200  memcpy (tmp, &(volData->buffer[2*(idz*volData->dims[0] *
201  volData->dims[1] + idy*volData->dims[0])]),
202  sizeof(unsigned char)*2*volData->dims[0]);
203  memcpy(&(volData->buffer[2*(idz * volData->dims[0] *
204  volData->dims[1] + idy * volData->dims[0])]),
205  &(volData->buffer[2*(idz * volData->dims[0] *
206  volData->dims[1] + iy * volData->dims[0])]),
207  sizeof(unsigned char)*2*volData->dims[0]);
208  memcpy(&(volData->buffer[2*(idz * volData->dims[0] *
209  volData->dims[1] + iy * volData->dims[0])]), tmp,
210  sizeof(unsigned char)*2*volData->dims[0]);
211  }
212  }
213  break;
214  }
215  }
216 
218  flag = 0;
219  // Recompute the normal direction vector
220  if (fabs(volData->dircos[2][0]) > fabs(volData->dircos[2][1])
221  &&
222  fabs(volData->dircos[2][0]) > fabs(volData->dircos[2][2])) {
223  if (volData->dircos[2][0] < 0) { // RL (21/05/2013 - Ting)
224  flag = 1;
225  volData->dircos[2][0] *= (-1); // LR
226  }
227  code += 0;
228  } else if (fabs(volData->dircos[2][1]) > fabs(volData->dircos[2][0])
229  &&
230  fabs(volData->dircos[2][1]) > fabs(volData->dircos[2][2])) {
231  if (volData->dircos[2][1] < 0) { // PA
232  flag = 1;
233  volData->dircos[2][1] *= (-1); // AP
234  }
235  code += 1;
236  } else {
237  if (volData->dircos[2][2] < 0) { // HF
238  flag = 1;
239  volData->dircos[2][2] *= (-1); // FH
240  }
241  code += 2;
242  }
243  if (flag) {
244  switch (volData->nbitsalloc) {
245  case 8:
246  for (iz=0,idz = volData->dims[2]-1; iz < idz; idz-=1,iz+=1) {
247  memcpy (tmp, &(volData->buffer[idz * volData->dims[0] * volData->dims[1]]),
248  sizeof(unsigned char)*(volData->dims[0] * volData->dims[1]));
249  memcpy (&(volData->buffer[idz * volData->dims[0] * volData->dims[1]]),
250  &(volData->buffer[iz * volData->dims[0] * volData->dims[1]]),
251  sizeof(unsigned char)*(volData->dims[0] * volData->dims[1]));
252  memcpy(&(volData->buffer[iz * volData->dims[0] * volData->dims[1]]), tmp,
253  sizeof(unsigned char)*(volData->dims[0] * volData->dims[1]));
254  }
255  break;
256  case 16:
257  for (iz=0,idz = volData->dims[2]-1; iz < idz; idz-=1,iz+=1) {
258  memcpy (tmp, &(volData->buffer[idz * 2 * volData->dims[0] * volData->dims[1]]),
259  sizeof(unsigned char)*2*(volData->dims[0] * volData->dims[1]));
260  memcpy (&(volData->buffer[idz * 2 * volData->dims[0] * volData->dims[1]]),
261  &(volData->buffer[iz * 2* volData->dims[0] * volData->dims[1]]),
262  sizeof(unsigned char)*2*(volData->dims[0] * volData->dims[1]));
263  memcpy (&(volData->buffer[iz * 2 * volData->dims[0] * volData->dims[1]]), tmp,
264  sizeof(unsigned char)*2*(volData->dims[0] * volData->dims[1]));
265  }
266  break;
267  }
268  }
269 
270  delete[] tmp;
271 
272  std::cout << "code = " << (int)code << std::endl;
273 
274  if (code == 6) return true;
275 
276  // Rotation to LPS
277  {
278  double dir[3][3];
279 
280  switch (volData->nbitsalloc) {
281  case 8:
282  tmp = new unsigned char[volData->dims[0]*volData->dims[1]*volData->dims[2]];
283  memcpy (tmp,volData->buffer,
284  sizeof(unsigned char)*(volData->dims[0] * volData->dims[1] * volData->dims[2]));
285  memset(volData->buffer, 0x00, volData->dims[0]*volData->dims[1]*volData->dims[2]);
286 
287  if (code == 24) { // PSL
288  std::cout << "PSL" << std::endl;
289  for (ix = 0, idy = 0; idy < volData->dims[1]; idy++) {
290  for (idx = 0; idx < volData->dims[0]; idx++) {
291  for (idz=0; idz < volData->dims[2]; idz++) {
292  volData->buffer[ix++] =
293  tmp[idz * volData->dims[0] * volData->dims[1] +
294  idy * volData->dims[0] + idx];
295  }
296  }
297  }
298  sx = volData->space[2]; sy = volData->space[0]; sz = volData->space[1];
299  dimx = volData->dims[2]; dimy = volData->dims[0]; dimz = volData->dims[1];
300  dir[0][0] = volData->dircos[2][0];
301  dir[0][1] = volData->dircos[2][1];
302  dir[0][2] = volData->dircos[2][2];
303  dir[1][0] = volData->dircos[0][0];
304  dir[1][1] = volData->dircos[0][1];
305  dir[1][2] = volData->dircos[0][2];
306  dir[2][0] = volData->dircos[1][0];
307  dir[2][1] = volData->dircos[1][1];
308  dir[2][2] = volData->dircos[1][2];
309  } else if (code == 33) { // SLP
310  std::cout << "SLP" << std::endl;
311  for (ix = 0, idx = 0; idx < volData->dims[0]; idx++) {
312  for (idz = 0; idz < volData->dims[2]; idz++) {
313  for (idy=0; idy < volData->dims[1]; idy++) {
314  volData->buffer[ix++] =
315  tmp[idz * volData->dims[0] * volData->dims[1] +
316  idy * volData->dims[0] + idx];
317  }
318  }
319  }
320  sx = volData->space[1]; sy = volData->space[2]; sz = volData->space[0];
321  dimx = volData->dims[1]; dimy = volData->dims[2]; dimz = volData->dims[0];
322  dir[0][0] = volData->dircos[1][0];
323  dir[0][1] = volData->dircos[1][1];
324  dir[0][2] = volData->dircos[1][2];
325  dir[1][0] = volData->dircos[2][0];
326  dir[1][1] = volData->dircos[2][1];
327  dir[1][2] = volData->dircos[2][2];
328  dir[2][0] = volData->dircos[0][0];
329  dir[2][1] = volData->dircos[0][1];
330  dir[2][2] = volData->dircos[0][2];
331  }
332  break;
333  case 16:
334  tmp = new unsigned char[2*volData->dims[0]*volData->dims[1]*volData->dims[2]];
335  memcpy (tmp,volData->buffer,
336  sizeof(unsigned char)*(2*volData->dims[0]*volData->dims[1]*volData->dims[2]));
337  memset(volData->buffer, 0x00, 2*volData->dims[0]*volData->dims[1]*volData->dims[2]);
338  if (code == 24) { // PSL
339  std::cout << "PSL" << std::endl;
340  for (ix=0, idy = 0; idy < volData->dims[1]; idy++) { //PLS
341  for (idx = 0; idx < volData->dims[0]; idx++) {
342  for (idz = 0; idz < volData->dims[2]; idz++) {
343  memcpy(&(volData->buffer[2*ix]),&(tmp[2*(idz*volData->dims[0]*volData->dims[1]+idy*volData->dims[0]+idx)]), sizeof(unsigned short));
344  ix++;
345  }
346  }
347  }
348  sx = volData->space[2]; sy = volData->space[0]; sz = volData->space[1];
349  dimx = volData->dims[2]; dimy = volData->dims[0]; dimz = volData->dims[1];
350  dir[0][0] = volData->dircos[2][0];
351  dir[0][1] = volData->dircos[2][1];
352  dir[0][2] = volData->dircos[2][2];
353  dir[1][0] = volData->dircos[0][0];
354  dir[1][1] = volData->dircos[0][1];
355  dir[1][2] = volData->dircos[0][2];
356  dir[2][0] = volData->dircos[1][0];
357  dir[2][1] = volData->dircos[1][1];
358  dir[2][2] = volData->dircos[1][2];
359  } else if (code == 33) { // SLP
360  std::cout << "SLP" << std::endl;
361  for (ix = 0, idx = 0; idx < volData->dims[0]; idx++) {
362  for (idz = 0; idz < volData->dims[2]; idz++) {
363  for (idy=0; idy < volData->dims[1]; idy++) {
364  memcpy(&(volData->buffer[2*ix]),&(tmp[2*(idz*volData->dims[0]*volData->dims[1]+idy*volData->dims[0]+idx)]), sizeof(unsigned short));
365  ix++;
366  }
367  }
368  }
369  sx = volData->space[1]; sy = volData->space[2]; sz = volData->space[0];
370  dimx = volData->dims[1]; dimy = volData->dims[2]; dimz = volData->dims[0];
371  dir[0][0] = volData->dircos[1][0];
372  dir[0][1] = volData->dircos[1][1];
373  dir[0][2] = volData->dircos[1][2];
374  dir[1][0] = volData->dircos[2][0];
375  dir[1][1] = volData->dircos[2][1];
376  dir[1][2] = volData->dircos[2][2];
377  dir[2][0] = volData->dircos[0][0];
378  dir[2][1] = volData->dircos[0][1];
379  dir[2][2] = volData->dircos[0][2];
380  }
381  break;
382  }
383 
384  // Update the dimensions, spacings and directions
385  volData->space[0]=sx; volData->space[1]=sy; volData->space[2]=sz;
386  volData->dims[0]=dimx; volData->dims[1]=dimy; volData->dims[2]=dimz;
387  volData->dircos[0][0] = dir[0][0]; volData->dircos[0][1] = dir[0][1]; volData->dircos[0][2] = dir[0][2];
388  volData->dircos[1][0] = dir[1][0]; volData->dircos[1][1] = dir[1][1]; volData->dircos[1][2] = dir[1][2];
389  volData->dircos[2][0] = dir[2][0]; volData->dircos[2][1] = dir[2][1]; volData->dircos[2][2] = dir[2][2];
390 
391  delete [] tmp;
392  }
393 
394  return true;
395 }
396 
397 bool Import::ReadPixelData(const gdcm::Image &image, const gdcm::DataSet &ds,
398  const gdcm::PixelFormat::ScalarType & stype,
399  unsigned short nbitsalloc, unsigned short nbitsstored,
400  float slope, float intercept,
401  unsigned int *len, char **buf,
402  int *umin, int *umax)
403 {
404  // Get Pixeldata
405  gdcm::Tag rawTag(0x7fe0, 0x0010); // Default to Pixel Data
406  const gdcm::DataElement& pixel_data = ds.GetDataElement(rawTag);
407  const gdcm::ByteValue *bv = pixel_data.GetByteValue();
408  *len = bv->GetLength();
409 
410  //std::cout << "Image length: " << *len << std::endl;
411  int j, step, nShiftBits;
412  unsigned int length;
413  unsigned char mask, v1;
414  unsigned short ivalue;
415  int value;
416  char *tmpbuf8;
417  unsigned short *tmpbuf16;
418 
419  switch (stype) {
420  case 0: // UINT8
421  case 1: // INT8
422  tmpbuf8 = new char[*len];
423  bv->GetBuffer(tmpbuf8, *len);
424  for(unsigned int di = 0; di < *len; di += 1)
425  {
426  if (stype == 0)
427  value = (int)((reinterpret_cast<unsigned char*>(tmpbuf8))[di]);
428  else if (stype == 1)
429  value = (int)((reinterpret_cast<char*>(tmpbuf8))[di]);
430 
431  if(*umax < value) *umax = value;
432  if(*umin > value) *umin = value;
433  }
434  *buf = tmpbuf8;
435  break;
436  case 2: // UINT12
437  break;
438  case 3: // INT12
439  break;
440  case 4: // UINT16
441  case 5: // INT16
442  length = (*len)/2;
443  nShiftBits = nbitsalloc - nbitsstored;
444  step = nbitsalloc / (sizeof(char) * 8);
445  mask = (255 >> nShiftBits);
446 
447  tmpbuf16 = new unsigned short[length];
448  tmpbuf8 = new char[*len];
449  bv->GetBuffer(reinterpret_cast<char*>(tmpbuf8), *len);
450 
451  for(unsigned int di = 0; di < *len; di += step)
452  {
453  // Little endian
455  ivalue = 0;
456  v1 = mask & ((unsigned char) tmpbuf8[di+step-1]);
458  for (j = step - 1; j > 0; j--) {
459  ivalue += (v1 * pow(256, j));
460  v1 = (unsigned char) tmpbuf8[di + j - 1];
461  }
462  ivalue += v1;
463 
464  tmpbuf16[di/2] = ivalue;
465 
466  if (stype == 4)
467  value = ivalue;
468  else if (stype == 5)
469  value = (int)(*(reinterpret_cast<short *>(&ivalue)));
470 
471  if(*umax < value) *umax = value;
472  if(*umin > value) *umin = value;
473  }
474  delete[] tmpbuf8;
475  *buf = reinterpret_cast<char*>(tmpbuf16);
476  break;
477  case 6: // UINT32
478  break;
479  case 7: // INT32
480  break;
481  case 8: // FLOAT16
482  break;
483  case 9: // FLOAT32
484  break;
485  case 10: // FLOAT64
486  break;
487  case 11: // UNKONOWN
488  break;
489  }
490  return 1;
491 }
492 
493 int Import::ProcessOneFile( std::string const & filename,
494  gdcm::Defs const & defs ,
495  std::string *series_desc_str,
496  std::string *patient_name_str,
497  std::string *patient_code_str,
498  std::string *scalar,
499  gdcm::PixelFormat::ScalarType *stype,
500  ImgFormat *volData)
501 {
502  gdcm::ImageReader reader;
503  reader.SetFileName( filename.c_str() );
504  if( !reader.Read() )
505  {
506  std::cerr << "Could not read image from: "
507  << filename << std::endl;
508  return 0;
509  }
510  const gdcm::File &file = reader.GetFile();
511  const gdcm::DataSet &ds = file.GetDataSet();
512  const gdcm::Image &image = reader.GetImage();
513 
514  // Get Series Description
515  gdcm::Tag att(0x0008,0x103e);
516  if( ds.FindDataElement(att) ) {
517  const gdcm::DataElement& series_desc = ds.GetDataElement( att );
518  //const gdcm::ByteValue * bv = series_desc.GetByteValue();
519  (*series_desc_str) += (series_desc.GetByteValue()->GetPointer());
520  (*series_desc_str).resize(series_desc.GetVL());
521  std::cout << "Serie Description: " << (*series_desc_str)
522  << std::endl;
523  }
524  // Get Patient Name
525  if( ds.FindDataElement(gdcm::Tag(0x0010, 0x0010)) ) {
526  const gdcm::DataElement& patient_name =
527  ds.GetDataElement( gdcm::Tag(0x0010, 0x0010) );
528  (*patient_name_str) += (patient_name.GetByteValue()->GetPointer());
529  (*patient_name_str).resize(patient_name.GetVL());
530  std::cout << "Patient Name: " << (*patient_name_str) << std::endl;
531  }
532  // Get Patient ID
533  if( ds.FindDataElement(gdcm::Tag(0x0010, 0x0020)) ) {
534  const gdcm::DataElement& patient_code =
535  ds.GetDataElement(gdcm::Tag(0x0010, 0x0020) );
536  (*patient_code_str) += (patient_code.GetByteValue()->GetPointer());
537  (*patient_code_str).resize(patient_code.GetVL());
538  std::cout << "Patient Code: " << (*patient_code_str)
539  << std::endl;
540  }
541 
542  // Read data set information for each file
543  // Get Slice Location (<0x0020,0x1041>)
544  if(ds.FindDataElement(gdcm::Tag(0x0020, 0x1041)) ) {
545  const gdcm::DataElement& location_ptr=
546  ds.GetDataElement( gdcm::Tag(0x0020, 0x1041) );
547  std::string location_str
548  (location_ptr.GetByteValue()->GetPointer());
549  location_str.resize(location_ptr.GetVL());
550  const float location = atof(location_str.c_str());
551  std::cout << "Location : " << location << std::endl ;
552  }
553  // Get Dimension Sizes
554  const unsigned int ndim = image.GetNumberOfDimensions();
555  std::cout << "Number Of Dimensions: " << ndim << std::endl;
556 
557  const unsigned int *dims=image.GetDimensions();
558  volData->dims[0] = dims[1];
559  volData->dims[1] = dims[0];
560  if (ndim == 3)
561  volData->dims[2] = dims[2];
562  else if (ndim == 2)
563  volData->dims[2] = 0;
564  std::cout << "Dimensions: (" << volData->dims[0] << ","
565  << volData->dims[1] << ","
566  << volData->dims[2] << ")" << std::endl ;
567  // Get Origin (<0x0020,0x0032>)
568  const double *origin = image.GetOrigin();
569 
570  volData->origin[0] = origin[0];
571  volData->origin[1] = origin[1];
572  volData->origin[2] = origin[2];
573  std::cout << "Origin: (";
574  int i;
575  for(i = 0; i < 2; ++i)
576  {
577  std::cout << origin[i] << ",";
578  }
579  std::cout << origin[i] << ")" << std::endl;
580  // Get Pixel Spacing (<0x0028, 0x0030>)
581  volData->space[2] = 0.0;
582  const double *space = image.GetSpacing();
583  volData->space[0] = space[0];
584  volData->space[1] = space[1];
585  volData->space[2] = space[2]; // (21/05/2013 - Ting)
586  // Get Slice Thickness
587  if (fabs(volData->space[2]) < 1.e-7) { // Bug fixed - 21/05/2013
588  const gdcm::DataElement& slice_thickness =
589  ds.GetDataElement(gdcm::Tag(0x0018, 0x0050));
590  const gdcm::DataElement& spacing_between_slices =
591  ds.GetDataElement(gdcm::Tag(0x0018, 0x0088));
592  if(!slice_thickness.IsEmpty()) {
593  std::string slice_thickness_str
594  (slice_thickness.GetByteValue()->GetPointer());
595  slice_thickness_str.resize(slice_thickness.GetVL());
596  volData->space[2] = atof(slice_thickness_str.c_str());
597  } else if (spacing_between_slices.IsEmpty()) {
598  std::string spacing_between_slices_str
599  (spacing_between_slices.GetByteValue()->GetPointer());
600  spacing_between_slices_str.resize(spacing_between_slices.GetVL());
601  volData->space[2] = atof(spacing_between_slices_str.c_str());
602  }
603  }
604  std::cout << "Spacing: (" << volData->space[0] << ","
605  << volData->space[1] << ","
606  << volData->space[2] << ")" << std::endl;
607  // Get Pixel Format
608  const gdcm::PixelFormat &ptype=image.GetPixelFormat();
609  // Get Samples Per Pixel (<0x0028,0x0002>)
610  volData->samples = ptype.GetSamplesPerPixel();
611  std::cout << "Samples per pixel: " << volData->samples << std::endl;
612  // Get Bits Allocated (<0x0028,0x0100>)
613  volData->nbitsalloc = ptype.GetBitsAllocated();
614  std::cout << "Bits Allocated: " << volData->nbitsalloc << std::endl;
615  // Get Bits Stored (<0x0028,0x0101>)
616  volData->nbitsstored = ptype.GetBitsStored();
617  std::cout << "Bits Stored: " << volData->nbitsstored << std::endl;
618  // Get High Bit (<0x0028,0x0102>)
619  volData->nhighbit = ptype.GetHighBit();
620  std::cout << "High Bit: " << volData->nhighbit << std::endl;
621  // Get Pixel Representation (<0x0028,0x0103>)
622  const unsigned short pixelRep = ptype.GetPixelRepresentation();
623  std::cout << "Pixel Representation: " << pixelRep << std::endl;
624  // Get Scalar Type (???)
625  *stype = ptype.GetScalarType();
626  const char *stypestring = ptype.GetScalarTypeAsString();
627  *scalar += stypestring;
628  std::cout << "Scalar Type: " << *stype << ", " << *scalar << std::endl;
629  // Get Orientation (<0x0020,0x0035> ou <0x0020,0x0037>)
630  const double *dircos = image.GetDirectionCosines();
631  volData->dircos[0][0] = dircos[0];
632  volData->dircos[0][1] = dircos[1];
633  volData->dircos[0][2] = dircos[2];
634  volData->dircos[1][0] = dircos[3];
635  volData->dircos[1][1] = dircos[4];
636  volData->dircos[1][2] = dircos[5];
637  volData->dircos[2][0] = dircos[1]*dircos[5]-dircos[2]*dircos[4];
638  volData->dircos[2][1] = dircos[2]*dircos[3]-dircos[0]*dircos[5];
639  volData->dircos[2][2] = dircos[0]*dircos[4]-dircos[1]*dircos[3];
640  std::cout << "X: (" << volData->dircos[0][0] << "," << volData->dircos[0][1] << "," << volData->dircos[0][2] << ")" << std::endl;
641  std::cout << "Y: (" << volData->dircos[1][0] << "," << volData->dircos[1][1] << "," << volData->dircos[1][2] << ")" << std::endl;
642  std::cout << "Z: (" << volData->dircos[2][0] << "," << volData->dircos[2][1] << "," << volData->dircos[2][2] << ")" << std::endl;
643 
644  gdcm::Orientation::OrientationType type =
645  gdcm::Orientation::GetType(dircos);
646  std::cout << "Type : " << type << std::endl;
647  const char *label = gdcm::Orientation::GetLabel( type );
648  std::cout << "Orientation Label: " << label << std::endl;
649  // Get Photometric Interpretation (<0x0028,0x0004>)
650  const gdcm::PhotometricInterpretation &pi =
651  image.GetPhotometricInterpretation();
652  std::cout << "PhotometricInterpretation: " << pi << std::endl;
653 
654  // Get Rescale Slope (<0x0028,0x1053>)
655  float slope, intercept;
656 
657  if(ds.FindDataElement(gdcm::Tag(0x0028, 0x1053)) ) {
658  const gdcm::DataElement& slope_ptr=
659  ds.GetDataElement( gdcm::Tag(0x0028, 0x1053) );
660  std::string slope_str (slope_ptr.GetByteValue()->GetPointer());
661  slope_str.resize(slope_ptr.GetVL());
662  volData->slope = atof(slope_str.c_str());
663  } else
664  volData->slope = image.GetSlope();
665 
666  // Get Rescale Intercept (<0x0028,0x1052>)
667  if(ds.FindDataElement(gdcm::Tag(0x0028, 0x1052)) ) {
668  const gdcm::DataElement& intercept_ptr=
669  ds.GetDataElement( gdcm::Tag(0x0028, 0x1052) );
670  std::string intercept_str
671  (intercept_ptr.GetByteValue()->GetPointer());
672  intercept_str.resize(intercept_ptr.GetVL());
673  volData->intercept = atof(intercept_str.c_str());
674  } else
675  volData->intercept = image.GetIntercept();
676 
677  ReadPixelData(image, ds, *stype, volData->nbitsalloc,
678  volData->nbitsstored, volData->slope,
679  volData->intercept, &volData->length, &volData->buffer,
680  &volData->umin,
681  &volData->umax);
682 
683  return 1;
684 }
685 
686 bool Import::ValidateMediaStorageIsImage (std::string const & filename,
687  gdcm::Defs const & defs )
688 {
689  gdcm::Reader reader;
690  reader.SetFileName( filename.c_str() );
691  if( !reader.Read() )
692  {
693  std::cerr << "Failed to read: " << filename << std::endl;
694  return 0;
695  }
696  const gdcm::File &file = reader.GetFile();
697  const gdcm::DataSet &ds = file.GetDataSet();
698 
699  gdcm::MediaStorage ms;
700  ms.SetFromFile(file);
701  /*
702  * Until gdcm::MediaStorage is fixed only *compile* time constant
703  * will be handled
704  * see -> http://chuckhahm.com/Ischem/Zurich/XX_0134
705  * which make gdcm::UIDs useless :(
706  */
707  if( ms.IsUndefined() )
708  {
709  std::cerr << "Unknown MediaStorage" << std::endl;
710  return 0;
711  }
712 
713  gdcm::UIDs uid;
714  uid.SetFromUID( ms.GetString() );
715  // std::cout << "MediaStorage is " << ms << " [" << uid.GetName()
716  // << "]" << std::endl;
717  const gdcm::TransferSyntax &ts =
718  file.GetHeader().GetDataSetTransferSyntax();
719  uid.SetFromUID( ts.GetString() );
720  // std::cout << "TransferSyntax is " << ts << " [" << uid.GetName()
721  // << "]" << std::endl;
722 
723  return gdcm::MediaStorage::IsImage( ms );
724 }
725 
726 int Import::ImportFile(std::string &filename, const gdcm::Defs &defs,
727  ImgFormat *volData) {
728  std::cout << "A file is processed ..." << std::endl;
729  unsigned int len;
730  std::string series_desc_str=""; // empty string
731  std::string patient_name_str=""; // empty string
732  std::string patient_code_str=""; // empty string
733  std::string type = "";
734 
735  series_desc_str.clear();
736  patient_name_str.clear();
737  patient_code_str.clear();
738  type.clear();
739  gdcm::PixelFormat::ScalarType stype;
740 
741  if (ValidateMediaStorageIsImage (filename, defs)) {
742  gdcm::Reader reader;
743  reader.SetFileName( filename.c_str() );
744  if( reader.Read() ) {
745  const gdcm::File &file = reader.GetFile();
746  const gdcm::DataSet &ds = file.GetDataSet();
747  // Read Series Instance UID
748  gdcm::Tag att(0x0020,0x000e); // Series Instance UID
749  const gdcm::DataElement& series_value= ds.GetDataElement( att );
750  if( ds.FindDataElement(att) ) {
751  const gdcm::ByteValue *bv = series_value.GetByteValue();
752  len = 0;
753  // Get Series Description
754  if( ds.FindDataElement(gdcm::Tag(0x0008,0x103e)) ) {
755  const gdcm::DataElement& series_desc =
756  ds.GetDataElement(gdcm::Tag(0x0008,0x103e));
757  (series_desc_str) +=
758  (series_desc.GetByteValue()->GetPointer());
759  (series_desc_str).resize(series_desc.GetVL());
760  std::cout << "Serie Description: " <<
761  (series_desc_str) << std::endl;
762  }
763 
764  series_desc_str.clear();
765 
766  if (ProcessOneFile(filename, defs, &series_desc_str,
767  &patient_name_str, &patient_code_str,
768  &type, &stype, volData)) {
769  return 1;
770  }
771  }
772  }
773  }
774  return 0;
775 }
776 
777 /*==========================================================
778  Public Methods
779  ===========================================================*/
780 int Import::DICOMImage (std::string &filename, ImgFormat *volData)
781 {
782  std::string xmlpath;
783 
784  if( !gdcm::System::FileExists(filename.c_str()) )
785  {
786  std::cerr << filename << " does not exist on the system"
787  << std::endl;
788  return 0;
789  }
790 
791  // Return the sinleton instance
792  gdcm::Global& g = gdcm::Global::GetInstance();
793 
794  // Locate the XML dict
795 
796  const char *xmlpathenv = getenv("GDCM_RESOURCES_PATH");
797  if( xmlpathenv )
798  {
799  // Make sure to look for XML dict in user explicitly
800  // specified dir first:
801  xmlpath = xmlpathenv;
802  // Prepend path at the begining of the path list
803  if( !g.Prepend( xmlpath.c_str() ) )
804  {
805  std::cerr << "specified Resources Path is not valid: "
806  << xmlpath << std::endl;
807  return 0;
808  }
809  }
810 
811 
812  // Retrieve the default/internal (Part 3)
813  const gdcm::Defs &defs = g.GetDefs();
814 
815  int res = 0;
816 
817  if(gdcm::System::FileIsDirectory(filename.c_str()) ) {
818  std::cout << "Reading file directory is not supported ..." << std::endl;
819  } else {
820  res = ImportFile (filename, defs, volData);
821  LPS_ReorientImageVolume(volData);
822  }
823 
824  return res;
825 }
826 
int DICOMImage(std::string &filename, ImgFormat *imgformat)
Destructor.