1 /*****************************************************************************
5 ** Purpose: Interface to CTN Dicom classes
6 ** Programmer: Kevin Rosenberg
7 ** Date Started: March 2001
9 ** This is part of the CTSim program
10 ** Copyright (c) 1983-2009 Kevin Rosenberg
12 ** This program is free software; you can redistribute it and/or modify
13 ** it under the terms of the GNU General Public License (version 2) as
14 ** published by the Free Software Foundation.
16 ** This program is distributed in the hope that it will be useful,
17 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
18 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 ** GNU General Public License for more details.
21 ** You should have received a copy of the GNU General Public License
22 ** along with this program; if not, write to the Free Software
23 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24 ******************************************************************************/
33 #include "imagefile.h"
34 #include "projections.h"
37 DicomImporter::DicomImporter (const char* const pszFile)
38 : m_strFilename(pszFile), m_bFail(false), m_iContents(DICOM_CONTENTS_INVALID),
39 m_pImageFile(NULL), m_pProjections(NULL), m_pFile(NULL)
41 unsigned long lOptions = DCM_ORDERLITTLEENDIAN | DCM_FORMATCONVERSION;
44 CONDITION cond = DCM_OpenFile (pszFile, lOptions, &m_pFile);
45 if (cond != DCM_NORMAL) {
48 CONDITION cond2 = COND_TopCondition (&cond, textbuf, sizeof(textbuf));
49 cond2 = DCM_NORMAL; // testing
50 if (cond2 != DCM_NORMAL) {
51 m_strFailMessage = "DCM_OpenFile failure: ";
52 m_strFailMessage += m_strFilename;
54 m_strFailMessage = textbuf;
59 unsigned short iNRows, iNCols, iBitsAllocated, iBitsStored, iHighBit, iPixRep;
60 DCM_ELEMENT aElemRequired[] = {
61 {DCM_IMGROWS, DCM_US, "", 1, sizeof(iNRows),
62 {reinterpret_cast<char*>(&iNRows)}},
63 {DCM_IMGCOLUMNS, DCM_US, "", 1, sizeof(iNCols),
64 {reinterpret_cast<char*>(&iNCols)}},
65 {DCM_IMGBITSALLOCATED, DCM_US, "", 1, sizeof(iBitsAllocated),
66 {reinterpret_cast<char*>(&iBitsAllocated)}},
67 {DCM_IMGBITSSTORED, DCM_US, "", 1, sizeof(iBitsStored),
68 {reinterpret_cast<char*>(&iBitsStored)}},
69 {DCM_IMGHIGHBIT, DCM_US, "", 1, sizeof(iHighBit),
70 {reinterpret_cast<char*>(&iHighBit)}},
71 {DCM_IMGPIXELREPRESENTATION, DCM_US, "", 1, sizeof(iPixRep),
72 {reinterpret_cast<char*>(&iPixRep)}},
74 int nElemRequired = sizeof (aElemRequired) / sizeof(DCM_ELEMENT);
76 if (DCM_ParseObject (&m_pFile, aElemRequired, nElemRequired, NULL, 0, NULL) == DCM_NORMAL) {
77 loadImage (iNRows, iNCols, iBitsAllocated, iBitsStored, iHighBit, iPixRep);
81 DCM_TAG somatomTag = DCM_MAKETAG(TAG_GROUP_SOMATOM, TAG_MEMBER_SOMATOM_DATA);
82 if (DCM_GetElementSize (&m_pFile, somatomTag, &lRtnLength) == DCM_NORMAL)
87 DicomImporter::loadImage(unsigned short iNRows, unsigned short iNCols, unsigned short iBitsAllocated,
88 unsigned short iBitsStored, unsigned short iHighBit, unsigned short iPixRep)
91 unsigned short iSamplesPerPixel, iPlanarConfig;
93 DCM_ELEMENT elemPlanarConfig =
94 {DCM_IMGPLANARCONFIGURATION, DCM_US, "", 1, sizeof(iPlanarConfig),
95 {reinterpret_cast<char*>(&iPlanarConfig)}};
96 DCM_ELEMENT elemSamplesPerPixel =
97 {DCM_IMGSAMPLESPERPIXEL, DCM_US, "", 1,
98 sizeof(iSamplesPerPixel), {reinterpret_cast<char*>(&iSamplesPerPixel)}};
100 if (DCM_ParseObject (&m_pFile, &elemSamplesPerPixel, 1, NULL, 0, NULL) != DCM_NORMAL)
101 iSamplesPerPixel = 1; // default value
103 if (iSamplesPerPixel > 1) {
105 if (DCM_GetElementValue (&m_pFile, &elemPlanarConfig, &lRtnLength, &ctx) != DCM_NORMAL) {
107 m_strFailMessage = "Planar Configuration not specified when iSamplesPerPixel > 1";
111 char szRescaleSlope[17];
112 char szRescaleIntercept[17];
113 double dRescaleSlope = 1;
114 double dRescaleIntercept = 0;
115 DCM_ELEMENT elemRescaleSlope =
116 {DCM_IMGRESCALESLOPE, DCM_DS, "", 1, strlen(szRescaleSlope),
118 DCM_ELEMENT elemRescaleIntercept =
119 {DCM_IMGRESCALEINTERCEPT, DCM_DS, "", 1, strlen(szRescaleIntercept),
120 {szRescaleIntercept}};
121 if (DCM_ParseObject (&m_pFile, &elemRescaleSlope, 1, NULL, 0, NULL) == DCM_NORMAL) {
122 if (sscanf (szRescaleSlope, "%lf", &dRescaleSlope) != 1)
123 sys_error (ERR_SEVERE, "Error parsing rescale slope");
125 if (DCM_ParseObject (&m_pFile, &elemRescaleIntercept, 1, NULL, 0, NULL) == DCM_NORMAL) {
126 if (sscanf (szRescaleIntercept, "%lf", &dRescaleIntercept) != 1)
127 sys_error (ERR_SEVERE, "Error parsing rescale intercept");
130 DCM_ELEMENT elemPixelData = {DCM_PXLPIXELDATA, DCM_OT, "", 1, 0, {NULL}};
131 // Get the actual pixel data (the only other required element)
132 if (DCM_GetElementSize (&m_pFile, elemPixelData.tag, &lRtnLength) != DCM_NORMAL) {
134 m_strFailMessage = "Can't get pixel data size";
138 // Check the size of the pixel data to make sure we have the correct amount...
139 unsigned long lRealLength = lRtnLength;
140 unsigned long lCheckLengthInBytes = iNRows * iNCols * iSamplesPerPixel;
141 if (iBitsAllocated == 16)
142 lCheckLengthInBytes *= 2;
143 if (lCheckLengthInBytes > lRealLength) {
145 m_strFailMessage = "Too little pixel data supplied";
148 // Now allocate the memory to hold the pixel data and get it from the DICOM file...
150 unsigned char* pRawPixels = new unsigned char [lCheckLengthInBytes];
151 elemPixelData.length = lCheckLengthInBytes;
152 elemPixelData.d.ot = pRawPixels;
155 CONDITION cond = DCM_GetElementValue (&m_pFile, &elemPixelData, &lRtnLength, &ctx);
156 if ((cond != DCM_NORMAL) && (cond != DCM_GETINCOMPLETE)) {
158 m_strFailMessage = "Can't read pixel data";
162 if ((lCheckLengthInBytes < lRealLength) && (cond != DCM_GETINCOMPLETE)) {
164 m_strFailMessage = "Should have gooten incomplete message reading pixel data";
169 m_pImageFile = new ImageFile (iNCols, iNRows);
170 ImageFileArray v = m_pImageFile->getArray();
171 double dScale = 1 << iBitsStored;
172 unsigned int iMaskLength = iBitsStored;
175 unsigned int iMask = (1 << iMaskLength) - 1;
176 for (int iy = 0; iy < iNRows; iy++) {
177 for (int ix = 0; ix < iNCols; ix++) {
178 if (iBitsAllocated == 8) {
179 unsigned char cV = pRawPixels[iy * iNRows + ix];
180 v[ix][iy] = (cV & iMask) / dScale;
181 } else if (iBitsAllocated == 16) {
182 unsigned long lBase = (iy * iNRows + ix) * 2;
183 unsigned char cV1 = pRawPixels[lBase];
184 unsigned char cV2 = pRawPixels[lBase+1] & iMask;
185 int iV = cV1 + (cV2 << 8);
186 v[ix][iNRows - 1 - iy] = iV * dRescaleSlope + dRescaleIntercept;
190 m_iContents = DICOM_CONTENTS_IMAGE;
194 DicomImporter::loadProjections()
199 unsigned short iNViews, iNDets;
200 DCM_ELEMENT aElemRequired[] = {
201 {DCM_IMGROWS, DCM_US, "", 1, sizeof(iNViews),
202 {reinterpret_cast<char*>(&iNViews)}},
203 {DCM_IMGCOLUMNS, DCM_US, "", 1, sizeof(iNDets),
204 {reinterpret_cast<char*>(&iNDets)}},
206 int nElemRequired = sizeof (aElemRequired) / sizeof(DCM_ELEMENT);
208 if (DCM_ParseObject (&m_pFile, aElemRequired, nElemRequired, NULL, 0, NULL) != DCM_NORMAL) {
210 m_strFailMessage = "Unable to read header for projections";
214 DCM_TAG somatomTag = DCM_MAKETAG(TAG_GROUP_SOMATOM, TAG_MEMBER_SOMATOM_DATA);
215 DCM_ELEMENT elemProjections = {somatomTag, DCM_UN, "", 1, 0, {NULL}};
216 if (DCM_GetElementSize (&m_pFile, elemProjections.tag, &lRtnLength) != DCM_NORMAL) {
218 m_strFailMessage = "Can't find projection data";
222 unsigned char* pRawProjections = new unsigned char [lRtnLength];
223 elemProjections.length = lRtnLength;
224 elemProjections.d.ot = pRawProjections;
227 CONDITION cond = DCM_GetElementValue (&m_pFile, &elemProjections, &lRtnLength, &ctx);
228 if ((cond != DCM_NORMAL) && (cond != DCM_GETINCOMPLETE)) {
230 m_strFailMessage = "Can't read projections data";
231 delete pRawProjections;
234 m_iContents = DICOM_CONTENTS_PROJECTIONS;
235 m_pProjections = new Projections;
236 if (! m_pProjections->initFromSomatomAR_STAR (iNViews, iNDets, pRawProjections, lRtnLength)) {
238 m_strFailMessage = "Error converting raw projection data";
239 delete m_pProjections;
240 m_pProjections = NULL;
243 delete pRawProjections;
247 DicomImporter::~DicomImporter()
250 DCM_CloseObject (&m_pFile);
254 DicomExporter::DicomExporter (ImageFile* pImageFile)
255 : m_pImageFile(pImageFile), m_pObject(NULL)
260 m_strFailMessage = "Initialized DicomExported with NULL imagefile";
263 m_bFail = ! createDicomObject();
266 DicomExporter::~DicomExporter ()
269 DCM_CloseObject (&m_pObject);
273 DicomExporter::writeFile (const char* const pszFilename)
278 m_strFilename = pszFilename;
280 CONDITION cond = DCM_WriteFile (&m_pObject, DCM_ORDERLITTLEENDIAN, pszFilename);
281 if (cond != DCM_NORMAL) {
283 m_strFailMessage = "Error writing DICOM file ";
284 m_strFailMessage += pszFilename;
292 DicomExporter::createDicomObject()
294 CONDITION cond = DCM_CreateObject (&m_pObject, 0);
295 if (cond != DCM_NORMAL) {
297 m_strFailMessage = "Error creating DICOM object";
302 m_pImageFile->getMinMax (dMin, dMax);
303 double dWidth = dMax - dMin;
306 double dScale = 65535. / dWidth;
308 double dRescaleIntercept = dMin;
309 double dRescaleSlope = 1 / dScale;
310 char szRescaleIntercept[17];
311 char szRescaleSlope[17];
312 snprintf (szRescaleIntercept, sizeof(szRescaleIntercept), "%e", dRescaleIntercept);
313 snprintf (szRescaleSlope, sizeof(szRescaleIntercept), "%e", dRescaleSlope);
315 char szCTSimRoot[] = "1.2.826.0.1.3680043.2.284.";
316 char szModality[] = "CT";
317 char szSOPClassUID[65] = "1.2.840.10008.5.4.1.1.2";
318 char szImgPhotometricInterp[] = "MONOCHROME2";
319 char szPixelSpacing[33] = "0\\0";
320 char szRelImageOrientationPatient[100] = "1\\0\\0\\0\\1\\0";
321 char szRelImagePositionPatient[49] = "0\\0\\0";
322 char szAcqKvp[] = "0";
323 char szRelAcquisitionNumber[] = "1";
324 char szRelImageNumber[] = "1";
325 char szIDSOPInstanceUID[65] = "";
326 char szIDManufacturer[] = "CTSim";
327 char szRelPositionRefIndicator[] = "0";
328 char szRelFrameOfReferenceUID[65] = "";
329 char szRelSeriesNumber[] = "1";
330 char szIDAccessionNumber[] = "0";
331 char szRelStudyID[] = "1";
332 char szIDReferringPhysician[] = "NONE";
333 char szIDStudyTime[] = "000000.0";
334 char szIDStudyDate[] = "00000000";
335 char szRelStudyInstanceUID[65] = "";
336 char szPatSex[] = "O";
337 char szPatBirthdate[] = "0000000";
338 char szPatID[] = "NONE";
339 char szPatName[] = "NONE";
340 char szIDImageType[] = "ORIGINAL";
341 char szIDManufacturerModel[65] = "";
343 std::ostringstream osPatComments;
344 m_pImageFile->printLabelsBrief (osPatComments);
345 size_t sizePatComments = osPatComments.str().length();
346 char* pszPatComments = new char [sizePatComments+1];
347 strncpy (pszPatComments, osPatComments.str().c_str(), sizePatComments);
349 snprintf (szIDSOPInstanceUID, sizeof(szIDSOPInstanceUID), "%s.2.1.6.1", szCTSimRoot);
350 snprintf (szRelStudyInstanceUID, sizeof(szRelStudyInstanceUID), "%s.2.1.6.1.1", szCTSimRoot);
351 snprintf (szRelFrameOfReferenceUID, sizeof(szRelFrameOfReferenceUID), "%s.99", szCTSimRoot);
353 snprintf (szIDManufacturerModel, sizeof(szIDManufacturerModel), "VERSION %s", VERSION);
355 snprintf (szPixelSpacing, sizeof(szPixelSpacing), "%e\\%e", m_pImageFile->axisIncrementX(), m_pImageFile->axisIncrementY());
356 double minX, maxX, minY, maxY;
357 if (m_pImageFile->getAxisExtent(minX, maxX, minY, maxY)) {
358 minX += m_pImageFile->axisIncrementX() / 2;
359 minY += m_pImageFile->axisIncrementY() / 2;
360 snprintf(szRelImagePositionPatient, sizeof(szRelImagePositionPatient), "%e\\%e\\0", minX, minY);
363 unsigned short iNRows = m_pImageFile->ny();
364 unsigned short iNCols = m_pImageFile->nx();
365 unsigned short iBitsAllocated = 16;
366 unsigned short iBitsStored = 16;
367 unsigned short iHighBit = 15;
368 unsigned short iPixRep = 0;
369 unsigned short iSamplesPerPixel = 1;
370 DCM_ELEMENT aElemRequired[] = {
371 {DCM_IMGROWS, DCM_US, "", 1, sizeof(iNRows),
372 {reinterpret_cast<char*>(&iNRows)}},
373 {DCM_IMGCOLUMNS, DCM_US, "", 1, sizeof(iNCols),
374 {reinterpret_cast<char*>(&iNCols)}},
375 {DCM_IMGBITSALLOCATED, DCM_US, "", 1, sizeof(iBitsAllocated),
376 {reinterpret_cast<char*>(&iBitsAllocated)}},
377 {DCM_IMGBITSSTORED, DCM_US, "", 1, sizeof(iBitsStored),
378 {reinterpret_cast<char*>(&iBitsStored)}},
379 {DCM_IMGHIGHBIT, DCM_US, "", 1, sizeof(iHighBit),
380 {reinterpret_cast<char*>(&iHighBit)}},
381 {DCM_IMGPIXELREPRESENTATION, DCM_US, "", 1, sizeof(iPixRep),
382 {reinterpret_cast<char*>(&iPixRep)}},
383 {DCM_IMGSAMPLESPERPIXEL, DCM_US, "", 1, sizeof(iSamplesPerPixel),
384 {reinterpret_cast<char*>(&iSamplesPerPixel)}},
385 {DCM_IMGRESCALESLOPE, DCM_DS, "", 1, strlen(szRescaleSlope),
387 {DCM_IMGRESCALEINTERCEPT, DCM_DS, "", 1, strlen(szRescaleIntercept),
388 {szRescaleIntercept}},
389 {DCM_IMGPHOTOMETRICINTERP, DCM_CS, "", 1, strlen(szImgPhotometricInterp),
390 {szImgPhotometricInterp}},
391 {DCM_IMGPIXELSPACING, DCM_DS, "", 1, strlen(szPixelSpacing),
393 {DCM_RELIMAGEORIENTATIONPATIENT, DCM_DS, "", 1,
394 strlen(szRelImageOrientationPatient), {szRelImageOrientationPatient}},
395 {DCM_RELIMAGEPOSITIONPATIENT, DCM_DS, "", 1,
396 strlen(szRelImagePositionPatient), {szRelImagePositionPatient}},
397 {DCM_ACQKVP, DCM_DS, "", 1, strlen(szAcqKvp), {szAcqKvp}},
398 {DCM_RELACQUISITIONNUMBER, DCM_IS, "", 1, strlen(szRelAcquisitionNumber),
399 {szRelAcquisitionNumber}},
400 {DCM_ACQSLICETHICKNESS, DCM_DS, "", 1, strlen(szRelAcquisitionNumber),
401 {szRelAcquisitionNumber}},
402 {DCM_RELIMAGENUMBER, DCM_IS, "", 1, strlen(szRelImageNumber),
404 {DCM_IDSOPINSTANCEUID, DCM_UI, "", 1, strlen(szIDSOPInstanceUID),
405 {szIDSOPInstanceUID}},
406 {DCM_IDMANUFACTURER, DCM_LO, "", 1, strlen(szIDManufacturer),
408 {DCM_RELPOSITIONREFINDICATOR, DCM_LO, "", 1,
409 strlen(szRelPositionRefIndicator), {szRelPositionRefIndicator}},
410 {DCM_RELFRAMEOFREFERENCEUID, DCM_UI, "", 1,
411 strlen(szRelFrameOfReferenceUID), {szRelFrameOfReferenceUID}},
412 {DCM_RELSERIESNUMBER, DCM_IS, "", 1, strlen(szRelSeriesNumber),
413 {szRelSeriesNumber}},
414 {DCM_RELSERIESINSTANCEUID, DCM_UI, "", 1, strlen(szIDAccessionNumber),
415 {szIDAccessionNumber}},
416 {DCM_IDACCESSIONNUMBER, DCM_SH, "", 1, strlen(szIDAccessionNumber),
417 {szIDAccessionNumber}},
418 {DCM_RELSTUDYID, DCM_SH, "", 1, strlen(szRelStudyID), {szRelStudyID}},
419 {DCM_IDREFERRINGPHYSICIAN, DCM_PN, "", 1, strlen(szIDReferringPhysician),
420 {szIDReferringPhysician}},
421 {DCM_IDSTUDYTIME, DCM_TM, "", 1, strlen(szIDStudyTime), {szIDStudyTime}},
422 {DCM_IDSTUDYDATE, DCM_DA, "", 1, strlen(szIDStudyDate), {szIDStudyDate}},
423 {DCM_RELSTUDYINSTANCEUID, DCM_UI, "", 1, strlen(szRelStudyInstanceUID),
424 {szRelStudyInstanceUID}},
425 {DCM_PATSEX, DCM_CS, "", 1, strlen(szPatSex), {szPatSex}},
426 {DCM_PATBIRTHDATE, DCM_DA, "", 1, strlen(szPatBirthdate), {szPatBirthdate}},
427 {DCM_PATID, DCM_LO, "", 1, strlen(szPatID), {szPatID}},
428 {DCM_PATNAME, DCM_PN, "", 1, strlen(szPatName), {szPatName}},
429 {DCM_IDIMAGETYPE, DCM_CS, "", 1, strlen(szIDImageType), {szIDImageType}},
430 {DCM_IDMODALITY, DCM_CS, "", 1, strlen(szModality), {szModality}},
431 {DCM_IDSOPCLASSUID, DCM_UI, "", 1, strlen(szSOPClassUID), {szSOPClassUID}},
432 {DCM_IDMANUFACTURERMODEL, DCM_LO, "", 1, strlen(szIDManufacturerModel),
433 {szIDManufacturerModel}},
434 {DCM_PATCOMMENTS, DCM_LT, "", 1, strlen(pszPatComments), {pszPatComments}},
436 int nElemRequired = sizeof (aElemRequired) / sizeof(DCM_ELEMENT);
439 cond = DCM_ModifyElements (&m_pObject, aElemRequired, nElemRequired, NULL, 0, &iUpdateCount);
441 DCM_ELEMENT elemPixelData = {DCM_PXLPIXELDATA, DCM_OT, "", 1, 0, {NULL}};
443 unsigned long lRealLength = 2 * m_pImageFile->nx() * m_pImageFile->ny();
445 unsigned char* pRawPixels = new unsigned char [lRealLength];
446 elemPixelData.length = lRealLength;
447 elemPixelData.d.ot = pRawPixels;
449 ImageFileArray v = m_pImageFile->getArray();
450 for (int iy = 0; iy < iNRows; iy++) {
451 for (int ix = 0; ix < iNCols; ix++) {
452 unsigned long lBase = (iy * iNRows + ix) * 2;
453 unsigned int iValue = nearest<int>(dScale * (v[ix][iNRows - 1 - iy] - dMin));
454 pRawPixels[lBase] = iValue & 0xFF;
455 pRawPixels[lBase+1] = (iValue & 0xFF00) >> 8;
458 cond = DCM_ModifyElements (&m_pObject, &elemPixelData, 1, NULL, 0, &iUpdateCount);
461 delete pszPatComments;
463 if (cond != DCM_NORMAL || iUpdateCount != 1) {
465 m_strFailMessage = "Error modifying pixel data";
472 #endif // HAVE_CTN_DICOM