r639: no message
[ctsim.git] / libctsim / ctndicom.cpp
1 /*****************************************************************************
2 ** FILE IDENTIFICATION
3 **
4 **      Name:           ctndicomp.cpp
5 **  Purpose:      Interface to CTN Dicom classes
6 **      Programmer:   Kevin Rosenberg
7 **      Date Started: March 2001
8 **
9 **  This is part of the CTSim program
10 **  Copyright (c) 1983-2001 Kevin Rosenberg
11 **
12 **  $Id: ctndicom.cpp,v 1.12 2001/03/13 14:53:44 kevin Exp $
13 **
14 **  This program is free software; you can redistribute it and/or modify
15 **  it under the terms of the GNU General Public License (version 2) as
16 **  published by the Free Software Foundation.
17 **
18 **  This program is distributed in the hope that it will be useful,
19 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
20 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
21 **  GNU General Public License for more details.
22 **
23 **  You should have received a copy of the GNU General Public License
24 **  along with this program; if not, write to the Free Software
25 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
26 ******************************************************************************/
27
28 #ifdef HAVE_CONFIG_H
29 #include "config.h"
30 #endif
31
32 #ifdef HAVE_CTN_DICOM
33
34 #include "ctndicom.h" 
35 #include "imagefile.h"
36 #include "projections.h"
37
38
39 DicomImporter::DicomImporter (const char* const pszFile)
40   : m_pFile(NULL), m_strFilename(pszFile), m_bFail(false), m_iContents(DICOM_CONTENTS_INVALID), 
41   m_pImageFile(NULL), m_pProjections(NULL)
42 {
43   unsigned long lOptions = DCM_ORDERLITTLEENDIAN;
44   DCM_Debug (FALSE);
45   if (DCM_OpenFile (pszFile, lOptions, &m_pFile) != DCM_NORMAL) {
46     m_bFail = true;;
47     m_strFailMessage = "Can't open file ";
48     m_strFailMessage += m_strFilename;
49     return;
50   }
51
52   unsigned short iNRows, iNCols, iBitsAllocated, iBitsStored, iHighBit, iPixRep;
53   DCM_ELEMENT aElemRequired[] = {
54     {DCM_IMGROWS, DCM_US, "", 1, sizeof(iNRows), reinterpret_cast<char*>(&iNRows)},
55     {DCM_IMGCOLUMNS, DCM_US, "", 1, sizeof(iNCols), reinterpret_cast<char*>(&iNCols)},
56     {DCM_IMGBITSALLOCATED, DCM_US, "", 1, sizeof(iBitsAllocated), reinterpret_cast<char*>(&iBitsAllocated)},
57     {DCM_IMGBITSSTORED, DCM_US, "", 1, sizeof(iBitsStored), reinterpret_cast<char*>(&iBitsStored)},
58     {DCM_IMGHIGHBIT, DCM_US, "", 1, sizeof(iHighBit), reinterpret_cast<char*>(&iHighBit)},
59     {DCM_IMGPIXELREPRESENTATION, DCM_US, "", 1, sizeof(iPixRep), reinterpret_cast<char*>(&iPixRep)}
60   };
61   int nElemRequired = sizeof (aElemRequired) / sizeof(DCM_ELEMENT);
62
63   if (DCM_ParseObject (&m_pFile, aElemRequired, nElemRequired, NULL, 0, NULL) == DCM_NORMAL) {
64     loadImage (iNRows, iNCols, iBitsAllocated, iBitsStored, iHighBit, iPixRep);
65     return;
66   }
67   unsigned long lRtnLength;
68   DCM_TAG somatomTag = DCM_MAKETAG(TAG_GROUP_SOMATOM, TAG_MEMBER_SOMATOM_DATA);
69   if (DCM_GetElementSize (&m_pFile, somatomTag, &lRtnLength) == DCM_NORMAL)
70     loadProjections();
71 }
72
73 void
74 DicomImporter::loadImage(unsigned short iNRows, unsigned short iNCols, unsigned short iBitsAllocated, 
75                          unsigned short iBitsStored, unsigned short iHighBit, unsigned short iPixRep)
76 {
77   unsigned long lRtnLength;
78   unsigned short iSamplesPerPixel, iPlanarConfig;
79   
80   DCM_ELEMENT elemPlanarConfig = {DCM_IMGPLANARCONFIGURATION, DCM_US, "", 1, sizeof(iPlanarConfig), 
81                       reinterpret_cast<char*>(&iPlanarConfig)};
82   DCM_ELEMENT elemSamplesPerPixel = {DCM_IMGSAMPLESPERPIXEL, DCM_US, "", 1, 
83                       sizeof(iSamplesPerPixel), reinterpret_cast<char*>(&iSamplesPerPixel)};
84
85   if (DCM_ParseObject (&m_pFile, &elemSamplesPerPixel, 1, NULL, 0, NULL) != DCM_NORMAL)
86     iSamplesPerPixel = 1;  // default value
87
88   if (iSamplesPerPixel > 1) {
89     void* ctx = NULL;
90     if (DCM_GetElementValue (&m_pFile, &elemPlanarConfig, &lRtnLength, &ctx) != DCM_NORMAL) {
91       m_bFail = true;
92       m_strFailMessage = "Planar Configuration not specified when iSamplesPerPixel > 1";
93     }
94   }
95
96   char szRescaleSlope[17];
97   char szRescaleIntercept[17];
98   double dRescaleSlope = 1;
99   double dRescaleIntercept = 0;
100   DCM_ELEMENT elemRescaleSlope = {DCM_IMGRESCALESLOPE, DCM_DS, "", 1, strlen(szRescaleSlope), szRescaleSlope};
101   DCM_ELEMENT elemRescaleIntercept = {DCM_IMGRESCALEINTERCEPT, DCM_DS, "", 1, strlen(szRescaleIntercept), szRescaleIntercept};
102   if (DCM_ParseObject (&m_pFile, &elemRescaleSlope, 1, NULL, 0, NULL) == DCM_NORMAL) {
103     if (sscanf (szRescaleSlope, "%lf", &dRescaleSlope) != 1)
104       sys_error (ERR_SEVERE, "Error parsing rescale slope");
105   }
106   if (DCM_ParseObject (&m_pFile, &elemRescaleIntercept, 1, NULL, 0, NULL) == DCM_NORMAL) {
107     if (sscanf (szRescaleIntercept, "%lf", &dRescaleIntercept) != 1)
108       sys_error (ERR_SEVERE, "Error parsing rescale intercept");
109   }
110
111   DCM_ELEMENT elemPixelData = {DCM_PXLPIXELDATA, DCM_OT, "", 1, 0, NULL};
112   // Get the actual  pixel data (the only other required element)
113   if (DCM_GetElementSize (&m_pFile, elemPixelData.tag, &lRtnLength) != DCM_NORMAL) {
114     m_bFail = true;
115     m_strFailMessage = "Can't get pixel data size";
116     return;
117   }
118
119   // Check the size of the pixel data to make sure we have the correct amount...
120   unsigned long lRealLength = lRtnLength;
121   unsigned long lCheckLengthInBytes = iNRows * iNCols * iSamplesPerPixel;
122   if (iBitsAllocated == 16)
123     lCheckLengthInBytes *= 2;
124   if (lCheckLengthInBytes > lRealLength) {
125     m_bFail = true;
126     m_strFailMessage = "Too little pixel data supplied";
127     return;
128   }
129   // Now allocate the memory to hold the pixel data and get it from the DICOM file...
130
131   unsigned char* pRawPixels = new unsigned char [lCheckLengthInBytes];
132   elemPixelData.length = lCheckLengthInBytes;
133   elemPixelData.d.ot = pRawPixels;
134   
135   void* ctx = NULL;
136   CONDITION cond = DCM_GetElementValue (&m_pFile, &elemPixelData, &lRtnLength, &ctx);
137   if ((cond != DCM_NORMAL) && (cond != DCM_GETINCOMPLETE)) {
138     m_bFail = true;
139     m_strFailMessage = "Can't read pixel data";
140     delete pRawPixels;
141     return;
142   }
143   if ((lCheckLengthInBytes < lRealLength) && (cond != DCM_GETINCOMPLETE)) {
144     m_bFail;
145     m_strFailMessage = "Should have gooten incomplete message reading pixel data";
146     delete pRawPixels;
147     return;
148   }
149
150   m_pImageFile = new ImageFile (iNCols, iNRows);
151   ImageFileArray v = m_pImageFile->getArray();
152   double dScale = 1 << iBitsStored;
153   unsigned int iMaskLength = iBitsStored;
154   if (iMaskLength > 8)
155     iMaskLength -= 8;
156   unsigned int iMask = (1 << iMaskLength) - 1;
157   for (int iy = 0; iy < iNRows; iy++) {
158     for (int ix = 0; ix < iNCols; ix++) {
159       if (iBitsAllocated == 8)  {
160         unsigned char cV = pRawPixels[iy * iNRows + ix];
161         v[ix][iy] = (cV & iMask) / dScale;
162       } else if (iBitsAllocated == 16) {
163         unsigned long lBase = (iy * iNRows + ix) * 2;
164         unsigned char cV1 = pRawPixels[lBase];
165         unsigned char cV2 = pRawPixels[lBase+1] & iMask;
166         int iV = cV1 + (cV2 << 8);
167         v[ix][iNRows - 1 - iy] = iV * dRescaleSlope + dRescaleIntercept;
168       }
169     }
170   }
171   m_iContents = DICOM_CONTENTS_IMAGE;
172 }
173
174 void
175 DicomImporter::loadProjections()
176 {
177   unsigned long lRtnLength;
178   void* ctx = NULL;
179
180   unsigned short iNViews, iNDets;
181   DCM_ELEMENT aElemRequired[] = {
182     {DCM_IMGROWS, DCM_US, "", 1, sizeof(iNViews), reinterpret_cast<char*>(&iNViews)},
183     {DCM_IMGCOLUMNS, DCM_US, "", 1, sizeof(iNDets), reinterpret_cast<char*>(&iNDets)},
184   };
185   int nElemRequired = sizeof (aElemRequired) / sizeof(DCM_ELEMENT);
186
187   if (DCM_ParseObject (&m_pFile, aElemRequired, nElemRequired, NULL, 0, NULL) != DCM_NORMAL) {
188     m_bFail = true;
189     m_strFailMessage = "Unable to read header for projections";
190     return;
191   }
192
193   DCM_TAG somatomTag = DCM_MAKETAG(TAG_GROUP_SOMATOM, TAG_MEMBER_SOMATOM_DATA);
194   DCM_ELEMENT elemProjections = {somatomTag, DCM_UNKNOWN, "", 1, 0, NULL};
195   if (DCM_GetElementSize (&m_pFile, elemProjections.tag, &lRtnLength) != DCM_NORMAL) {
196     m_bFail = true;
197     m_strFailMessage = "Can't find projection data";
198     return;
199   }
200
201   unsigned char* pRawProjections = new unsigned char [lRtnLength];
202   elemProjections.length = lRtnLength;
203   elemProjections.d.ot = pRawProjections;
204   
205   ctx = NULL;
206   CONDITION cond = DCM_GetElementValue (&m_pFile, &elemProjections, &lRtnLength, &ctx);
207   if ((cond != DCM_NORMAL) && (cond != DCM_GETINCOMPLETE)) {
208     m_bFail = true;
209     m_strFailMessage = "Can't read projections data";
210     delete pRawProjections;
211     return;
212   }
213   m_iContents = DICOM_CONTENTS_PROJECTIONS;
214   m_pProjections = new Projections;
215   if (! m_pProjections->initFromSomatomAR_STAR (iNViews, iNDets, pRawProjections, lRtnLength)) {
216     m_bFail = true;
217     m_strFailMessage = "Error converting raw projection data";
218     delete m_pProjections;
219     m_pProjections = NULL;
220   }
221
222   delete pRawProjections;
223 }
224
225
226 DicomImporter::~DicomImporter()
227 {
228   if (m_pFile)
229     DCM_CloseObject (&m_pFile);
230 }
231
232
233 DicomExporter::DicomExporter (ImageFile* pImageFile)
234 : m_pImageFile(pImageFile), m_pObject(NULL)
235 {
236   DCM_Debug (FALSE);
237   if (! pImageFile) {
238     m_bFail = true;
239     m_strFailMessage = "Initialized DicomExported with NULL imagefile";
240     return;
241   }
242   m_bFail = ! createDicomObject();
243 }
244
245 DicomExporter::~DicomExporter ()
246 {
247   if (m_pObject)
248     DCM_CloseObject (&m_pObject);
249 }
250
251 bool
252 DicomExporter::writeFile (const char* const pszFilename)
253 {
254   if (! m_pObject)
255     return false;
256
257   m_strFilename = pszFilename;
258
259   CONDITION cond = DCM_WriteFile (&m_pObject, DCM_ORDERLITTLEENDIAN, pszFilename);
260   if (cond != DCM_NORMAL) {
261     m_bFail = true;
262     m_strFailMessage = "Error writing DICOM file ";
263     m_strFailMessage += pszFilename;
264     return false;
265   }
266
267   return true;
268 }
269
270 bool
271 DicomExporter::createDicomObject()
272 {
273   CONDITION cond = DCM_CreateObject (&m_pObject, 0);
274   if (cond != DCM_NORMAL) {
275     m_bFail = true;
276     m_strFailMessage = "Error creating DICOM object";
277     return false;
278   }
279
280   double dMin, dMax;
281   m_pImageFile->getMinMax (dMin, dMax);
282   double dWidth = dMax - dMin;
283   if (dWidth == 0.)
284     dWidth = 1E-7;
285   double dScale = 65535. / dWidth;
286
287   double dRescaleIntercept = dMin;
288   double dRescaleSlope = 1 / dScale;
289   char szRescaleIntercept[17];
290   char szRescaleSlope[17];
291   snprintf (szRescaleIntercept, sizeof(szRescaleIntercept), "%e", dRescaleIntercept);
292   snprintf (szRescaleSlope, sizeof(szRescaleIntercept), "%e", dRescaleSlope);
293
294   char szCTSimRoot[] = "1.2.826.0.1.3680043.2.284.";
295   char szModality[] = "CT";
296   char szSOPClassUID[65] = "1.2.840.10008.5.4.1.1.2";
297   char szImgPhotometricInterp[] = "MONOCHROME2";
298   char szPixelSpacing[33] = "0\\0";
299   char szRelImageOrientationPatient[100] = "1\\0\\0\\0\\1\\0";
300   char szRelImagePositionPatient[49] = "0\\0\\0";
301   char szAcqKvp[] = "0";
302   char szRelAcquisitionNumber[] = "1";
303   char szRelImageNumber[] = "1";
304   char szIDSOPInstanceUID[65] = "";
305   char szIDManufacturer[] = "CTSim";
306   char szRelPositionRefIndicator[] = "0";
307   char szRelFrameOfReferenceUID[65] = "";
308   char szRelSeriesNumber[] = "1";
309   char szIDAccessionNumber[] = "0";
310   char szRelStudyID[] = "1";
311   char szIDReferringPhysician[] = "NONE";
312   char szIDStudyTime[] = "000000.0";
313   char szIDStudyDate[] = "00000000";
314   char szRelStudyInstanceUID[65] = "";
315   char szPatSex[] = "O";
316   char szPatBirthdate[] = "0000000";
317   char szPatID[] = "NONE";
318   char szPatName[] = "NONE";
319   char szIDImageType[] = "ORIGINAL";
320   char szIDManufacturerModel[65] = "";
321
322   std::ostringstream osPatComments;
323   m_pImageFile->printLabelsBrief (osPatComments);
324   size_t sizePatComments = osPatComments.str().length();
325   char* pszPatComments = new char [sizePatComments+1];
326   strncpy (pszPatComments, osPatComments.str().c_str(), sizePatComments);
327
328   snprintf (szIDSOPInstanceUID, sizeof(szIDSOPInstanceUID), "%s.2.1.6.1", szCTSimRoot);
329   snprintf (szRelStudyInstanceUID, sizeof(szRelStudyInstanceUID), "%s.2.1.6.1.1", szCTSimRoot);
330   snprintf (szRelFrameOfReferenceUID, sizeof(szRelFrameOfReferenceUID), "%s.99", szCTSimRoot);
331 #ifdef VERSION
332   snprintf (szIDManufacturerModel, sizeof(szIDManufacturerModel), "VERSION %s", VERSION);
333 #endif
334   snprintf (szPixelSpacing, sizeof(szPixelSpacing), "%e\\%e", m_pImageFile->axisIncrementX(), m_pImageFile->axisIncrementY());
335   double minX, maxX, minY, maxY;
336   if (m_pImageFile->getAxisExtent(minX, maxX, minY, maxY)) {
337     minX += m_pImageFile->axisIncrementX() / 2;
338     minY += m_pImageFile->axisIncrementY() / 2;
339     snprintf(szRelImagePositionPatient, sizeof(szRelImagePositionPatient), "%e\\%e\\0", minX, minY);
340   }
341
342   unsigned short iNRows = m_pImageFile->ny();
343   unsigned short iNCols = m_pImageFile->nx();
344   unsigned short iBitsAllocated = 16;
345   unsigned short iBitsStored = 16;
346   unsigned short iHighBit = 15;
347   unsigned short iPixRep = 0;
348   unsigned short iSamplesPerPixel = 1;
349   DCM_ELEMENT aElemRequired[] = {
350     {DCM_IMGROWS, DCM_US, "", 1, sizeof(iNRows), reinterpret_cast<char*>(&iNRows)},
351     {DCM_IMGCOLUMNS, DCM_US, "", 1, sizeof(iNCols), reinterpret_cast<char*>(&iNCols)},
352     {DCM_IMGBITSALLOCATED, DCM_US, "", 1, sizeof(iBitsAllocated), reinterpret_cast<char*>(&iBitsAllocated)},
353     {DCM_IMGBITSSTORED, DCM_US, "", 1, sizeof(iBitsStored), reinterpret_cast<char*>(&iBitsStored)},
354     {DCM_IMGHIGHBIT, DCM_US, "", 1, sizeof(iHighBit), reinterpret_cast<char*>(&iHighBit)},
355     {DCM_IMGPIXELREPRESENTATION, DCM_US, "", 1, sizeof(iPixRep), reinterpret_cast<char*>(&iPixRep)},
356     {DCM_IMGSAMPLESPERPIXEL, DCM_US, "", 1, sizeof(iSamplesPerPixel), reinterpret_cast<char*>(&iSamplesPerPixel)},
357     {DCM_IMGRESCALESLOPE, DCM_DS, "", 1, strlen(szRescaleSlope), szRescaleSlope},
358     {DCM_IMGRESCALEINTERCEPT, DCM_DS, "", 1, strlen(szRescaleIntercept), szRescaleIntercept},
359     {DCM_IMGPHOTOMETRICINTERP, DCM_CS, "", 1, strlen(szImgPhotometricInterp), szImgPhotometricInterp},
360     {DCM_IMGPIXELSPACING, DCM_DS, "", 1, strlen(szPixelSpacing), szPixelSpacing},
361     {DCM_RELIMAGEORIENTATIONPATIENT, DCM_DS, "", 1, strlen(szRelImageOrientationPatient), szRelImageOrientationPatient},
362     {DCM_RELIMAGEPOSITIONPATIENT, DCM_DS, "", 1, strlen(szRelImagePositionPatient), szRelImagePositionPatient},
363     {DCM_ACQKVP, DCM_DS, "", 1, strlen(szAcqKvp), szAcqKvp},
364     {DCM_RELACQUISITIONNUMBER, DCM_IS, "", 1, strlen(szRelAcquisitionNumber), szRelAcquisitionNumber},
365     {DCM_ACQSLICETHICKNESS, DCM_DS, "", 1, strlen(szRelAcquisitionNumber), szRelAcquisitionNumber},
366     {DCM_RELIMAGENUMBER, DCM_IS, "", 1, strlen(szRelImageNumber), szRelImageNumber},
367     {DCM_IDSOPINSTANCEUID, DCM_UI, "", 1, strlen(szIDSOPInstanceUID), szIDSOPInstanceUID},
368     {DCM_IDMANUFACTURER, DCM_LO, "", 1, strlen(szIDManufacturer), szIDManufacturer},
369     {DCM_RELPOSITIONREFINDICATOR, DCM_LO, "", 1, strlen(szRelPositionRefIndicator), szRelPositionRefIndicator},
370     {DCM_RELFRAMEOFREFERENCEUID, DCM_UI, "", 1, strlen(szRelFrameOfReferenceUID), szRelFrameOfReferenceUID},
371     {DCM_RELSERIESNUMBER, DCM_IS, "", 1, strlen(szRelSeriesNumber), szRelSeriesNumber},
372     {DCM_RELSERIESINSTANCEUID, DCM_UI, "", 1, strlen(szIDAccessionNumber), szIDAccessionNumber},
373     {DCM_IDACCESSIONNUMBER, DCM_SH, "", 1, strlen(szIDAccessionNumber), szIDAccessionNumber},
374     {DCM_RELSTUDYID, DCM_SH, "", 1, strlen(szRelStudyID), szRelStudyID},
375     {DCM_IDREFERRINGPHYSICIAN, DCM_PN, "", 1, strlen(szIDReferringPhysician), szIDReferringPhysician},
376     {DCM_IDSTUDYTIME, DCM_TM, "", 1, strlen(szIDStudyTime), szIDStudyTime},
377     {DCM_IDSTUDYDATE, DCM_DA, "", 1, strlen(szIDStudyDate), szIDStudyDate},
378     {DCM_RELSTUDYINSTANCEUID, DCM_UI, "", 1, strlen(szRelStudyInstanceUID), szRelStudyInstanceUID},
379     {DCM_PATSEX, DCM_CS, "", 1, strlen(szPatSex), szPatSex},
380     {DCM_PATBIRTHDATE, DCM_DA, "", 1, strlen(szPatBirthdate), szPatBirthdate},
381     {DCM_PATID, DCM_LO, "", 1, strlen(szPatID), szPatID},
382     {DCM_PATNAME, DCM_PN, "", 1, strlen(szPatName), szPatName},
383     {DCM_IDIMAGETYPE, DCM_CS, "", 1, strlen(szIDImageType), szIDImageType},
384     {DCM_IDMODALITY, DCM_CS, "", 1, strlen(szModality), szModality},
385     {DCM_IDSOPCLASSUID, DCM_UI, "", 1, strlen(szSOPClassUID), szSOPClassUID},
386     {DCM_IDMANUFACTURERMODEL, DCM_LO, "", 1, strlen(szIDManufacturerModel), szIDManufacturerModel},
387     {DCM_PATCOMMENTS, DCM_LT, "", 1, strlen(pszPatComments), pszPatComments},
388   };
389   int nElemRequired = sizeof (aElemRequired) / sizeof(DCM_ELEMENT);
390
391   int iUpdateCount;
392   cond = DCM_ModifyElements (&m_pObject, aElemRequired, nElemRequired, NULL, 0, &iUpdateCount);
393
394   DCM_ELEMENT elemPixelData = {DCM_PXLPIXELDATA, DCM_OT, "", 1, 0, NULL};
395
396   unsigned long lRealLength = 2 * m_pImageFile->nx() * m_pImageFile->ny();
397
398   unsigned char* pRawPixels = new unsigned char [lRealLength];
399   elemPixelData.length = lRealLength;
400   elemPixelData.d.ot = pRawPixels;
401   
402   ImageFileArray v = m_pImageFile->getArray();
403   for (int iy = 0; iy < iNRows; iy++) {
404     for (int ix = 0; ix < iNCols; ix++) {
405         unsigned long lBase = (iy * iNRows + ix) * 2;
406         unsigned int iValue = nearest<int>(dScale * (v[ix][iNRows - 1 - iy] - dMin));
407         pRawPixels[lBase] = iValue & 0xFF;
408         pRawPixels[lBase+1] = (iValue & 0xFF00) >> 8;
409     }
410   }
411   cond = DCM_ModifyElements (&m_pObject, &elemPixelData, 1, NULL, 0, &iUpdateCount);
412   delete pRawPixels;
413
414   delete pszPatComments;
415
416   if (cond != DCM_NORMAL || iUpdateCount != 1) {
417     m_bFail = true;
418     m_strFailMessage = "Error modifying pixel data";
419     return false;
420   }
421
422   return true;
423 }
424
425 #endif // HAVE_CTN_DICOM
426
427