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-2001 Kevin Rosenberg
12 ** $Id: ctndicom.cpp,v 1.9 2001/03/07 17:33:38 kevin Exp $
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.
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.
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 ******************************************************************************/
35 #include "imagefile.h"
36 #include "projections.h"
39 DicomImporter::DicomImporter (const char* const pszFile)
40 : m_strFilename(pszFile), m_bFail(false), m_iContents(DICOM_CONTENTS_INVALID),
41 m_pImageFile(NULL), m_pProjections(NULL)
43 unsigned long lOptions = DCM_ORDERLITTLEENDIAN;
45 if (DCM_OpenFile (pszFile, lOptions, &m_pFile) != DCM_NORMAL) {
47 m_strFailMessage = "Can't open file ";
48 m_strFailMessage += m_strFilename;
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)}
61 int nElemRequired = sizeof (aElemRequired) / sizeof(DCM_ELEMENT);
63 if (DCM_ParseObject (&m_pFile, aElemRequired, nElemRequired, NULL, 0, NULL) == DCM_NORMAL) {
64 loadImage (iNRows, iNCols, iBitsAllocated, iBitsStored, iHighBit, iPixRep);
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)
74 DicomImporter::loadImage(unsigned short iNRows, unsigned short iNCols, unsigned short iBitsAllocated,
75 unsigned short iBitsStored, unsigned short iHighBit, unsigned short iPixRep)
77 unsigned long lRtnLength;
78 unsigned short iSamplesPerPixel, iPlanarConfig;
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)};
85 if (DCM_ParseObject (&m_pFile, &elemSamplesPerPixel, 1, NULL, 0, NULL) != DCM_NORMAL)
86 iSamplesPerPixel = 1; // default value
88 if (iSamplesPerPixel > 1) {
90 if (DCM_GetElementValue (&m_pFile, &elemPlanarConfig, &lRtnLength, &ctx) != DCM_NORMAL) {
92 m_strFailMessage = "Planar Configuration not specified when iSamplesPerPixel > 1";
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");
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");
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) {
115 m_strFailMessage = "Can't get pixel data size";
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) {
126 m_strFailMessage = "Too little pixel data supplied";
129 // Now allocate the memory to hold the pixel data and get it from the DICOM file...
131 unsigned char* pRawPixels = new unsigned char [lCheckLengthInBytes];
132 elemPixelData.length = lCheckLengthInBytes;
133 elemPixelData.d.ot = pRawPixels;
136 CONDITION cond = DCM_GetElementValue (&m_pFile, &elemPixelData, &lRtnLength, &ctx);
137 if ((cond != DCM_NORMAL) && (cond != DCM_GETINCOMPLETE)) {
139 m_strFailMessage = "Can't read pixel data";
143 if ((lCheckLengthInBytes < lRealLength) && (cond != DCM_GETINCOMPLETE)) {
145 m_strFailMessage = "Should have gooten incomplete message reading pixel data";
150 m_pImageFile = new ImageFile (iNCols, iNRows);
151 ImageFileArray v = m_pImageFile->getArray();
152 double dScale = 1 << iBitsStored;
153 unsigned int iMaskLength = iBitsStored;
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;
171 m_iContents = DICOM_CONTENTS_IMAGE;
175 DicomImporter::loadProjections()
177 unsigned long lRtnLength;
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)},
185 int nElemRequired = sizeof (aElemRequired) / sizeof(DCM_ELEMENT);
187 if (DCM_ParseObject (&m_pFile, aElemRequired, nElemRequired, NULL, 0, NULL) != DCM_NORMAL) {
189 m_strFailMessage = "Unable to read header for projections";
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) {
197 m_strFailMessage = "Can't find projection data";
201 unsigned char* pRawProjections = new unsigned char [lRtnLength];
202 elemProjections.length = lRtnLength;
203 elemProjections.d.ot = pRawProjections;
206 CONDITION cond = DCM_GetElementValue (&m_pFile, &elemProjections, &lRtnLength, &ctx);
207 if ((cond != DCM_NORMAL) && (cond != DCM_GETINCOMPLETE)) {
209 m_strFailMessage = "Can't read projections data";
210 delete pRawProjections;
213 m_iContents = DICOM_CONTENTS_PROJECTIONS;
214 m_pProjections = new Projections;
215 if (! m_pProjections->initFromSomatomAR_STAR (iNViews, iNDets, pRawProjections, lRtnLength)) {
217 m_strFailMessage = "Error converting raw projection data";
218 delete m_pProjections;
219 m_pProjections = NULL;
222 delete pRawProjections;
226 DicomImporter::~DicomImporter()
229 DCM_CloseObject (&m_pFile);
233 DicomExporter::DicomExporter (ImageFile* pImageFile)
234 : m_pImageFile(pImageFile), m_pObject(NULL)
239 m_strFailMessage = "Initialized DicomExported with NULL imagefile";
242 m_bFail = ! createDicomObject();
245 DicomExporter::~DicomExporter ()
248 DCM_CloseObject (&m_pObject);
252 DicomExporter::writeFile (const char* const pszFilename)
257 m_strFilename = pszFilename;
259 CONDITION cond = DCM_WriteFile (&m_pObject, DCM_ORDERLITTLEENDIAN, pszFilename);
260 if (cond != DCM_NORMAL) {
262 m_strFailMessage = "Error writing DICOM file ";
263 m_strFailMessage += pszFilename;
271 DicomExporter::createDicomObject()
273 CONDITION cond = DCM_CreateObject (&m_pObject, 0);
274 if (cond != DCM_NORMAL) {
276 m_strFailMessage = "Error creating DICOM object";
281 m_pImageFile->getMinMax (dMin, dMax);
282 double dWidth = dMax - dMin;
285 double dScale = 65535. / dWidth;
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);
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] = "";
322 snprintf (szIDSOPInstanceUID, sizeof(szIDSOPInstanceUID), "%s.2.1.6.1", szCTSimRoot);
323 snprintf (szRelStudyInstanceUID, sizeof(szRelStudyInstanceUID), "%s.2.1.6.1.1", szCTSimRoot);
324 snprintf (szRelFrameOfReferenceUID, sizeof(szRelFrameOfReferenceUID), "%s.99", szCTSimRoot);
326 snprintf (szIDManufacturerModel, sizeof(szIDManufacturerModel), "VERSION %s", VERSION);
328 snprintf (szPixelSpacing, sizeof(szPixelSpacing), "%e\\%e", m_pImageFile->axisIncrementX(), m_pImageFile->axisIncrementY());
329 double minX, maxX, minY, maxY;
330 if (m_pImageFile->getAxisExtent(minX, maxX, minY, maxY)) {
331 minX += m_pImageFile->axisIncrementX() / 2;
332 minY += m_pImageFile->axisIncrementY() / 2;
333 snprintf(szRelImagePositionPatient, sizeof(szRelImagePositionPatient), "%e\\%e\\0", minX, minY);
336 unsigned short iNRows = m_pImageFile->ny();
337 unsigned short iNCols = m_pImageFile->nx();
338 unsigned short iBitsAllocated = 16;
339 unsigned short iBitsStored = 16;
340 unsigned short iHighBit = 15;
341 unsigned short iPixRep = 0;
342 unsigned short iSamplesPerPixel = 1;
343 DCM_ELEMENT aElemRequired[] = {
344 {DCM_IMGROWS, DCM_US, "", 1, sizeof(iNRows), reinterpret_cast<char*>(&iNRows)},
345 {DCM_IMGCOLUMNS, DCM_US, "", 1, sizeof(iNCols), reinterpret_cast<char*>(&iNCols)},
346 {DCM_IMGBITSALLOCATED, DCM_US, "", 1, sizeof(iBitsAllocated), reinterpret_cast<char*>(&iBitsAllocated)},
347 {DCM_IMGBITSSTORED, DCM_US, "", 1, sizeof(iBitsStored), reinterpret_cast<char*>(&iBitsStored)},
348 {DCM_IMGHIGHBIT, DCM_US, "", 1, sizeof(iHighBit), reinterpret_cast<char*>(&iHighBit)},
349 {DCM_IMGPIXELREPRESENTATION, DCM_US, "", 1, sizeof(iPixRep), reinterpret_cast<char*>(&iPixRep)},
350 {DCM_IMGSAMPLESPERPIXEL, DCM_US, "", 1, sizeof(iSamplesPerPixel), reinterpret_cast<char*>(&iSamplesPerPixel)},
351 {DCM_IMGRESCALESLOPE, DCM_DS, "", 1, strlen(szRescaleSlope), szRescaleSlope},
352 {DCM_IMGRESCALEINTERCEPT, DCM_DS, "", 1, strlen(szRescaleIntercept), szRescaleIntercept},
353 {DCM_IMGPHOTOMETRICINTERP, DCM_CS, "", 1, strlen(szImgPhotometricInterp), szImgPhotometricInterp},
354 {DCM_IMGPIXELSPACING, DCM_DS, "", 1, strlen(szPixelSpacing), szPixelSpacing},
355 {DCM_RELIMAGEORIENTATIONPATIENT, DCM_DS, "", 1, strlen(szRelImageOrientationPatient), szRelImageOrientationPatient},
356 {DCM_RELIMAGEPOSITIONPATIENT, DCM_DS, "", 1, strlen(szRelImagePositionPatient), szRelImagePositionPatient},
357 {DCM_ACQKVP, DCM_DS, "", 1, strlen(szAcqKvp), szAcqKvp},
358 {DCM_RELACQUISITIONNUMBER, DCM_IS, "", 1, strlen(szRelAcquisitionNumber), szRelAcquisitionNumber},
359 {DCM_ACQSLICETHICKNESS, DCM_DS, "", 1, strlen(szRelAcquisitionNumber), szRelAcquisitionNumber},
360 {DCM_RELIMAGENUMBER, DCM_IS, "", 1, strlen(szRelImageNumber), szRelImageNumber},
361 {DCM_IDSOPINSTANCEUID, DCM_UI, "", 1, strlen(szIDSOPInstanceUID), szIDSOPInstanceUID},
362 {DCM_IDMANUFACTURER, DCM_LO, "", 1, strlen(szIDManufacturer), szIDManufacturer},
363 {DCM_RELPOSITIONREFINDICATOR, DCM_LO, "", 1, strlen(szRelPositionRefIndicator), szRelPositionRefIndicator},
364 {DCM_RELFRAMEOFREFERENCEUID, DCM_UI, "", 1, strlen(szRelFrameOfReferenceUID), szRelFrameOfReferenceUID},
365 {DCM_RELSERIESNUMBER, DCM_IS, "", 1, strlen(szRelSeriesNumber), szRelSeriesNumber},
366 {DCM_RELSERIESINSTANCEUID, DCM_UI, "", 1, strlen(szIDAccessionNumber), szIDAccessionNumber},
367 {DCM_IDACCESSIONNUMBER, DCM_SH, "", 1, strlen(szIDAccessionNumber), szIDAccessionNumber},
368 {DCM_RELSTUDYID, DCM_SH, "", 1, strlen(szRelStudyID), szRelStudyID},
369 {DCM_IDREFERRINGPHYSICIAN, DCM_PN, "", 1, strlen(szIDReferringPhysician), szIDReferringPhysician},
370 {DCM_IDSTUDYTIME, DCM_TM, "", 1, strlen(szIDStudyTime), szIDStudyTime},
371 {DCM_IDSTUDYDATE, DCM_DA, "", 1, strlen(szIDStudyDate), szIDStudyDate},
372 {DCM_RELSTUDYINSTANCEUID, DCM_UI, "", 1, strlen(szRelStudyInstanceUID), szRelStudyInstanceUID},
373 {DCM_PATSEX, DCM_CS, "", 1, strlen(szPatSex), szPatSex},
374 {DCM_PATBIRTHDATE, DCM_DA, "", 1, strlen(szPatBirthdate), szPatBirthdate},
375 {DCM_PATID, DCM_LO, "", 1, strlen(szPatID), szPatID},
376 {DCM_PATNAME, DCM_PN, "", 1, strlen(szPatName), szPatName},
377 {DCM_IDIMAGETYPE, DCM_CS, "", 1, strlen(szIDImageType), szIDImageType},
378 {DCM_IDMODALITY, DCM_CS, "", 1, strlen(szModality), szModality},
379 {DCM_IDSOPCLASSUID, DCM_UI, "", 1, strlen(szSOPClassUID), szSOPClassUID},
380 {DCM_IDMANUFACTURERMODEL, DCM_LO, "", 1, strlen(szIDManufacturerModel), szIDManufacturerModel},
382 int nElemRequired = sizeof (aElemRequired) / sizeof(DCM_ELEMENT);
385 cond = DCM_ModifyElements (&m_pObject, aElemRequired, nElemRequired, NULL, 0, &iUpdateCount);
388 DCM_ELEMENT elemPixelData = {DCM_PXLPIXELDATA, DCM_OT, "", 1, 0, NULL};
390 unsigned long lRealLength = 2 * m_pImageFile->nx() * m_pImageFile->ny();
392 unsigned char* pRawPixels = new unsigned char [lRealLength];
393 elemPixelData.length = lRealLength;
394 elemPixelData.d.ot = pRawPixels;
396 ImageFileArray v = m_pImageFile->getArray();
397 for (int iy = 0; iy < iNRows; iy++) {
398 for (int ix = 0; ix < iNCols; ix++) {
399 unsigned long lBase = (iy * iNRows + ix) * 2;
400 unsigned int iValue = nearest<int>(dScale * (v[ix][iNRows - 1 - iy] - dMin));
401 pRawPixels[lBase] = iValue & 0xFF;
402 pRawPixels[lBase+1] = (iValue & 0xFF00) >> 8;
405 cond = DCM_ModifyElements (&m_pObject, &elemPixelData, 1, NULL, 0, &iUpdateCount);
408 if (cond != DCM_NORMAL || iUpdateCount != 1) {
410 m_strFailMessage = "Error modifying pixel data";
417 #endif // HAVE_CTN_DICOM