X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fctndicom.cpp;h=156cebe2e342130f0fe139d5dcf4f1b874931765;hp=85ec88c44e56707dd2f3793033392c631d4631ca;hb=1a050c98763fbbc0662731b0b76953acede6f5d7;hpb=c358b8c8b5649f14e2b8203b999ba8549a244727 diff --git a/libctsim/ctndicom.cpp b/libctsim/ctndicom.cpp index 85ec88c..156cebe 100644 --- a/libctsim/ctndicom.cpp +++ b/libctsim/ctndicom.cpp @@ -1,15 +1,15 @@ /***************************************************************************** ** FILE IDENTIFICATION ** -** Name: ctndicomp.cpp +** Name: ctndicomp.cpp ** Purpose: Interface to CTN Dicom classes -** Programmer: Kevin Rosenberg -** Date Started: March 2001 +** Programmer: Kevin Rosenberg +** Date Started: March 2001 ** ** This is part of the CTSim program ** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: ctndicom.cpp,v 1.2 2001/03/01 20:02:18 kevin Exp $ +** $Id$ ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License (version 2) as @@ -30,7 +30,447 @@ #endif #ifdef HAVE_CTN_DICOM + #include "ctndicom.h" +#include "imagefile.h" +#include "projections.h" + + +DicomImporter::DicomImporter (const char* const pszFile) + : m_strFilename(pszFile), m_bFail(false), m_iContents(DICOM_CONTENTS_INVALID), + m_pImageFile(NULL), m_pProjections(NULL), m_pFile(NULL) +{ + unsigned long lOptions = DCM_ORDERLITTLEENDIAN | DCM_FORMATCONVERSION; + DCM_Debug (FALSE); + + CONDITION cond = DCM_OpenFile (pszFile, lOptions, &m_pFile); + if (cond != DCM_NORMAL) { + m_bFail = true; + char textbuf [2048]; + CONDITION cond2 = COND_TopCondition (&cond, textbuf, sizeof(textbuf)); + cond2 = DCM_NORMAL; // testing + if (cond2 != DCM_NORMAL) { + m_strFailMessage = "DCM_OpenFile failure: "; + m_strFailMessage += m_strFilename; + } else + m_strFailMessage = textbuf; + + return; + } + + unsigned short iNRows, iNCols, iBitsAllocated, iBitsStored, iHighBit, iPixRep; + DCM_ELEMENT aElemRequired[] = { + {DCM_IMGROWS, DCM_US, "", 1, sizeof(iNRows), + {reinterpret_cast(&iNRows)}}, + {DCM_IMGCOLUMNS, DCM_US, "", 1, sizeof(iNCols), + {reinterpret_cast(&iNCols)}}, + {DCM_IMGBITSALLOCATED, DCM_US, "", 1, sizeof(iBitsAllocated), + {reinterpret_cast(&iBitsAllocated)}}, + {DCM_IMGBITSSTORED, DCM_US, "", 1, sizeof(iBitsStored), + {reinterpret_cast(&iBitsStored)}}, + {DCM_IMGHIGHBIT, DCM_US, "", 1, sizeof(iHighBit), + {reinterpret_cast(&iHighBit)}}, + {DCM_IMGPIXELREPRESENTATION, DCM_US, "", 1, sizeof(iPixRep), + {reinterpret_cast(&iPixRep)}}, + }; + int nElemRequired = sizeof (aElemRequired) / sizeof(DCM_ELEMENT); + + if (DCM_ParseObject (&m_pFile, aElemRequired, nElemRequired, NULL, 0, NULL) == DCM_NORMAL) { + loadImage (iNRows, iNCols, iBitsAllocated, iBitsStored, iHighBit, iPixRep); + return; + } + U32 lRtnLength; + DCM_TAG somatomTag = DCM_MAKETAG(TAG_GROUP_SOMATOM, TAG_MEMBER_SOMATOM_DATA); + if (DCM_GetElementSize (&m_pFile, somatomTag, &lRtnLength) == DCM_NORMAL) + loadProjections(); +} + +void +DicomImporter::loadImage(unsigned short iNRows, unsigned short iNCols, unsigned short iBitsAllocated, + unsigned short iBitsStored, unsigned short iHighBit, unsigned short iPixRep) +{ + U32 lRtnLength; + unsigned short iSamplesPerPixel, iPlanarConfig; + + DCM_ELEMENT elemPlanarConfig = + {DCM_IMGPLANARCONFIGURATION, DCM_US, "", 1, sizeof(iPlanarConfig), + {reinterpret_cast(&iPlanarConfig)}}; + DCM_ELEMENT elemSamplesPerPixel = + {DCM_IMGSAMPLESPERPIXEL, DCM_US, "", 1, + sizeof(iSamplesPerPixel), {reinterpret_cast(&iSamplesPerPixel)}}; + + if (DCM_ParseObject (&m_pFile, &elemSamplesPerPixel, 1, NULL, 0, NULL) != DCM_NORMAL) + iSamplesPerPixel = 1; // default value + + if (iSamplesPerPixel > 1) { + void* ctx = NULL; + if (DCM_GetElementValue (&m_pFile, &elemPlanarConfig, &lRtnLength, &ctx) != DCM_NORMAL) { + m_bFail = true; + m_strFailMessage = "Planar Configuration not specified when iSamplesPerPixel > 1"; + } + } + + char szRescaleSlope[17]; + char szRescaleIntercept[17]; + double dRescaleSlope = 1; + double dRescaleIntercept = 0; + DCM_ELEMENT elemRescaleSlope = + {DCM_IMGRESCALESLOPE, DCM_DS, "", 1, strlen(szRescaleSlope), + {szRescaleSlope}}; + DCM_ELEMENT elemRescaleIntercept = + {DCM_IMGRESCALEINTERCEPT, DCM_DS, "", 1, strlen(szRescaleIntercept), + {szRescaleIntercept}}; + if (DCM_ParseObject (&m_pFile, &elemRescaleSlope, 1, NULL, 0, NULL) == DCM_NORMAL) { + if (sscanf (szRescaleSlope, "%lf", &dRescaleSlope) != 1) + sys_error (ERR_SEVERE, "Error parsing rescale slope"); + } + if (DCM_ParseObject (&m_pFile, &elemRescaleIntercept, 1, NULL, 0, NULL) == DCM_NORMAL) { + if (sscanf (szRescaleIntercept, "%lf", &dRescaleIntercept) != 1) + sys_error (ERR_SEVERE, "Error parsing rescale intercept"); + } + + DCM_ELEMENT elemPixelData = {DCM_PXLPIXELDATA, DCM_OT, "", 1, 0, {NULL}}; + // Get the actual pixel data (the only other required element) + if (DCM_GetElementSize (&m_pFile, elemPixelData.tag, &lRtnLength) != DCM_NORMAL) { + m_bFail = true; + m_strFailMessage = "Can't get pixel data size"; + return; + } + + // Check the size of the pixel data to make sure we have the correct amount... + unsigned long lRealLength = lRtnLength; + unsigned long lCheckLengthInBytes = iNRows * iNCols * iSamplesPerPixel; + if (iBitsAllocated == 16) + lCheckLengthInBytes *= 2; + if (lCheckLengthInBytes > lRealLength) { + m_bFail = true; + m_strFailMessage = "Too little pixel data supplied"; + return; + } + // Now allocate the memory to hold the pixel data and get it from the DICOM file... + + unsigned char* pRawPixels = new unsigned char [lCheckLengthInBytes]; + elemPixelData.length = lCheckLengthInBytes; + elemPixelData.d.ot = pRawPixels; + + void* ctx = NULL; + CONDITION cond = DCM_GetElementValue (&m_pFile, &elemPixelData, &lRtnLength, &ctx); + if ((cond != DCM_NORMAL) && (cond != DCM_GETINCOMPLETE)) { + m_bFail = true; + m_strFailMessage = "Can't read pixel data"; + delete pRawPixels; + return; + } + if ((lCheckLengthInBytes < lRealLength) && (cond != DCM_GETINCOMPLETE)) { + m_bFail = true; + m_strFailMessage = "Should have gooten incomplete message reading pixel data"; + delete pRawPixels; + return; + } + + m_pImageFile = new ImageFile (iNCols, iNRows); + ImageFileArray v = m_pImageFile->getArray(); + double dScale = 1 << iBitsStored; + unsigned int iMaskLength = iBitsStored; + if (iMaskLength > 8) + iMaskLength -= 8; + unsigned int iMask = (1 << iMaskLength) - 1; + for (int iy = 0; iy < iNRows; iy++) { + for (int ix = 0; ix < iNCols; ix++) { + if (iBitsAllocated == 8) { + unsigned char cV = pRawPixels[iy * iNRows + ix]; + v[ix][iy] = (cV & iMask) / dScale; + } else if (iBitsAllocated == 16) { + unsigned long lBase = (iy * iNRows + ix) * 2; + unsigned char cV1 = pRawPixels[lBase]; + unsigned char cV2 = pRawPixels[lBase+1] & iMask; + int iV = cV1 + (cV2 << 8); + v[ix][iNRows - 1 - iy] = iV * dRescaleSlope + dRescaleIntercept; + } + } + } + m_iContents = DICOM_CONTENTS_IMAGE; +} + +void +DicomImporter::loadProjections() +{ + U32 lRtnLength; + void* ctx = NULL; + + unsigned short iNViews, iNDets; + DCM_ELEMENT aElemRequired[] = { + {DCM_IMGROWS, DCM_US, "", 1, sizeof(iNViews), + {reinterpret_cast(&iNViews)}}, + {DCM_IMGCOLUMNS, DCM_US, "", 1, sizeof(iNDets), + {reinterpret_cast(&iNDets)}}, + }; + int nElemRequired = sizeof (aElemRequired) / sizeof(DCM_ELEMENT); + + if (DCM_ParseObject (&m_pFile, aElemRequired, nElemRequired, NULL, 0, NULL) != DCM_NORMAL) { + m_bFail = true; + m_strFailMessage = "Unable to read header for projections"; + return; + } + + DCM_TAG somatomTag = DCM_MAKETAG(TAG_GROUP_SOMATOM, TAG_MEMBER_SOMATOM_DATA); + DCM_ELEMENT elemProjections = {somatomTag, DCM_UN, "", 1, 0, {NULL}}; + if (DCM_GetElementSize (&m_pFile, elemProjections.tag, &lRtnLength) != DCM_NORMAL) { + m_bFail = true; + m_strFailMessage = "Can't find projection data"; + return; + } + + unsigned char* pRawProjections = new unsigned char [lRtnLength]; + elemProjections.length = lRtnLength; + elemProjections.d.ot = pRawProjections; + + ctx = NULL; + CONDITION cond = DCM_GetElementValue (&m_pFile, &elemProjections, &lRtnLength, &ctx); + if ((cond != DCM_NORMAL) && (cond != DCM_GETINCOMPLETE)) { + m_bFail = true; + m_strFailMessage = "Can't read projections data"; + delete pRawProjections; + return; + } + m_iContents = DICOM_CONTENTS_PROJECTIONS; + m_pProjections = new Projections; + if (! m_pProjections->initFromSomatomAR_STAR (iNViews, iNDets, pRawProjections, lRtnLength)) { + m_bFail = true; + m_strFailMessage = "Error converting raw projection data"; + delete m_pProjections; + m_pProjections = NULL; + } + + delete pRawProjections; +} + + +DicomImporter::~DicomImporter() +{ + if (m_pFile) + DCM_CloseObject (&m_pFile); +} + + +DicomExporter::DicomExporter (ImageFile* pImageFile) +: m_pImageFile(pImageFile), m_pObject(NULL) +{ + DCM_Debug (FALSE); + if (! pImageFile) { + m_bFail = true; + m_strFailMessage = "Initialized DicomExported with NULL imagefile"; + return; + } + m_bFail = ! createDicomObject(); +} + +DicomExporter::~DicomExporter () +{ + if (m_pObject) + DCM_CloseObject (&m_pObject); +} + +bool +DicomExporter::writeFile (const char* const pszFilename) +{ + if (! m_pObject) + return false; + + m_strFilename = pszFilename; + + CONDITION cond = DCM_WriteFile (&m_pObject, DCM_ORDERLITTLEENDIAN, pszFilename); + if (cond != DCM_NORMAL) { + m_bFail = true; + m_strFailMessage = "Error writing DICOM file "; + m_strFailMessage += pszFilename; + return false; + } + + return true; +} + +bool +DicomExporter::createDicomObject() +{ + CONDITION cond = DCM_CreateObject (&m_pObject, 0); + if (cond != DCM_NORMAL) { + m_bFail = true; + m_strFailMessage = "Error creating DICOM object"; + return false; + } + + double dMin, dMax; + m_pImageFile->getMinMax (dMin, dMax); + double dWidth = dMax - dMin; + if (dWidth == 0.) + dWidth = 1E-7; + double dScale = 65535. / dWidth; + + double dRescaleIntercept = dMin; + double dRescaleSlope = 1 / dScale; + char szRescaleIntercept[17]; + char szRescaleSlope[17]; + snprintf (szRescaleIntercept, sizeof(szRescaleIntercept), "%e", dRescaleIntercept); + snprintf (szRescaleSlope, sizeof(szRescaleIntercept), "%e", dRescaleSlope); + + char szCTSimRoot[] = "1.2.826.0.1.3680043.2.284."; + char szModality[] = "CT"; + char szSOPClassUID[65] = "1.2.840.10008.5.4.1.1.2"; + char szImgPhotometricInterp[] = "MONOCHROME2"; + char szPixelSpacing[33] = "0\\0"; + char szRelImageOrientationPatient[100] = "1\\0\\0\\0\\1\\0"; + char szRelImagePositionPatient[49] = "0\\0\\0"; + char szAcqKvp[] = "0"; + char szRelAcquisitionNumber[] = "1"; + char szRelImageNumber[] = "1"; + char szIDSOPInstanceUID[65] = ""; + char szIDManufacturer[] = "CTSim"; + char szRelPositionRefIndicator[] = "0"; + char szRelFrameOfReferenceUID[65] = ""; + char szRelSeriesNumber[] = "1"; + char szIDAccessionNumber[] = "0"; + char szRelStudyID[] = "1"; + char szIDReferringPhysician[] = "NONE"; + char szIDStudyTime[] = "000000.0"; + char szIDStudyDate[] = "00000000"; + char szRelStudyInstanceUID[65] = ""; + char szPatSex[] = "O"; + char szPatBirthdate[] = "0000000"; + char szPatID[] = "NONE"; + char szPatName[] = "NONE"; + char szIDImageType[] = "ORIGINAL"; + char szIDManufacturerModel[65] = ""; + + std::ostringstream osPatComments; + m_pImageFile->printLabelsBrief (osPatComments); + size_t sizePatComments = osPatComments.str().length(); + char* pszPatComments = new char [sizePatComments+1]; + strncpy (pszPatComments, osPatComments.str().c_str(), sizePatComments); + + snprintf (szIDSOPInstanceUID, sizeof(szIDSOPInstanceUID), "%s.2.1.6.1", szCTSimRoot); + snprintf (szRelStudyInstanceUID, sizeof(szRelStudyInstanceUID), "%s.2.1.6.1.1", szCTSimRoot); + snprintf (szRelFrameOfReferenceUID, sizeof(szRelFrameOfReferenceUID), "%s.99", szCTSimRoot); +#ifdef VERSION + snprintf (szIDManufacturerModel, sizeof(szIDManufacturerModel), "VERSION %s", VERSION); +#endif + snprintf (szPixelSpacing, sizeof(szPixelSpacing), "%e\\%e", m_pImageFile->axisIncrementX(), m_pImageFile->axisIncrementY()); + double minX, maxX, minY, maxY; + if (m_pImageFile->getAxisExtent(minX, maxX, minY, maxY)) { + minX += m_pImageFile->axisIncrementX() / 2; + minY += m_pImageFile->axisIncrementY() / 2; + snprintf(szRelImagePositionPatient, sizeof(szRelImagePositionPatient), "%e\\%e\\0", minX, minY); + } + + unsigned short iNRows = m_pImageFile->ny(); + unsigned short iNCols = m_pImageFile->nx(); + unsigned short iBitsAllocated = 16; + unsigned short iBitsStored = 16; + unsigned short iHighBit = 15; + unsigned short iPixRep = 0; + unsigned short iSamplesPerPixel = 1; + DCM_ELEMENT aElemRequired[] = { + {DCM_IMGROWS, DCM_US, "", 1, sizeof(iNRows), + {reinterpret_cast(&iNRows)}}, + {DCM_IMGCOLUMNS, DCM_US, "", 1, sizeof(iNCols), + {reinterpret_cast(&iNCols)}}, + {DCM_IMGBITSALLOCATED, DCM_US, "", 1, sizeof(iBitsAllocated), + {reinterpret_cast(&iBitsAllocated)}}, + {DCM_IMGBITSSTORED, DCM_US, "", 1, sizeof(iBitsStored), + {reinterpret_cast(&iBitsStored)}}, + {DCM_IMGHIGHBIT, DCM_US, "", 1, sizeof(iHighBit), + {reinterpret_cast(&iHighBit)}}, + {DCM_IMGPIXELREPRESENTATION, DCM_US, "", 1, sizeof(iPixRep), + {reinterpret_cast(&iPixRep)}}, + {DCM_IMGSAMPLESPERPIXEL, DCM_US, "", 1, sizeof(iSamplesPerPixel), + {reinterpret_cast(&iSamplesPerPixel)}}, + {DCM_IMGRESCALESLOPE, DCM_DS, "", 1, strlen(szRescaleSlope), + {szRescaleSlope}}, + {DCM_IMGRESCALEINTERCEPT, DCM_DS, "", 1, strlen(szRescaleIntercept), + {szRescaleIntercept}}, + {DCM_IMGPHOTOMETRICINTERP, DCM_CS, "", 1, strlen(szImgPhotometricInterp), + {szImgPhotometricInterp}}, + {DCM_IMGPIXELSPACING, DCM_DS, "", 1, strlen(szPixelSpacing), + {szPixelSpacing}}, + {DCM_RELIMAGEORIENTATIONPATIENT, DCM_DS, "", 1, + strlen(szRelImageOrientationPatient), {szRelImageOrientationPatient}}, + {DCM_RELIMAGEPOSITIONPATIENT, DCM_DS, "", 1, + strlen(szRelImagePositionPatient), {szRelImagePositionPatient}}, + {DCM_ACQKVP, DCM_DS, "", 1, strlen(szAcqKvp), {szAcqKvp}}, + {DCM_RELACQUISITIONNUMBER, DCM_IS, "", 1, strlen(szRelAcquisitionNumber), + {szRelAcquisitionNumber}}, + {DCM_ACQSLICETHICKNESS, DCM_DS, "", 1, strlen(szRelAcquisitionNumber), + {szRelAcquisitionNumber}}, + {DCM_RELIMAGENUMBER, DCM_IS, "", 1, strlen(szRelImageNumber), + {szRelImageNumber}}, + {DCM_IDSOPINSTANCEUID, DCM_UI, "", 1, strlen(szIDSOPInstanceUID), + {szIDSOPInstanceUID}}, + {DCM_IDMANUFACTURER, DCM_LO, "", 1, strlen(szIDManufacturer), + {szIDManufacturer}}, + {DCM_RELPOSITIONREFINDICATOR, DCM_LO, "", 1, + strlen(szRelPositionRefIndicator), {szRelPositionRefIndicator}}, + {DCM_RELFRAMEOFREFERENCEUID, DCM_UI, "", 1, + strlen(szRelFrameOfReferenceUID), {szRelFrameOfReferenceUID}}, + {DCM_RELSERIESNUMBER, DCM_IS, "", 1, strlen(szRelSeriesNumber), + {szRelSeriesNumber}}, + {DCM_RELSERIESINSTANCEUID, DCM_UI, "", 1, strlen(szIDAccessionNumber), + {szIDAccessionNumber}}, + {DCM_IDACCESSIONNUMBER, DCM_SH, "", 1, strlen(szIDAccessionNumber), + {szIDAccessionNumber}}, + {DCM_RELSTUDYID, DCM_SH, "", 1, strlen(szRelStudyID), {szRelStudyID}}, + {DCM_IDREFERRINGPHYSICIAN, DCM_PN, "", 1, strlen(szIDReferringPhysician), + {szIDReferringPhysician}}, + {DCM_IDSTUDYTIME, DCM_TM, "", 1, strlen(szIDStudyTime), {szIDStudyTime}}, + {DCM_IDSTUDYDATE, DCM_DA, "", 1, strlen(szIDStudyDate), {szIDStudyDate}}, + {DCM_RELSTUDYINSTANCEUID, DCM_UI, "", 1, strlen(szRelStudyInstanceUID), + {szRelStudyInstanceUID}}, + {DCM_PATSEX, DCM_CS, "", 1, strlen(szPatSex), {szPatSex}}, + {DCM_PATBIRTHDATE, DCM_DA, "", 1, strlen(szPatBirthdate), {szPatBirthdate}}, + {DCM_PATID, DCM_LO, "", 1, strlen(szPatID), {szPatID}}, + {DCM_PATNAME, DCM_PN, "", 1, strlen(szPatName), {szPatName}}, + {DCM_IDIMAGETYPE, DCM_CS, "", 1, strlen(szIDImageType), {szIDImageType}}, + {DCM_IDMODALITY, DCM_CS, "", 1, strlen(szModality), {szModality}}, + {DCM_IDSOPCLASSUID, DCM_UI, "", 1, strlen(szSOPClassUID), {szSOPClassUID}}, + {DCM_IDMANUFACTURERMODEL, DCM_LO, "", 1, strlen(szIDManufacturerModel), + {szIDManufacturerModel}}, + {DCM_PATCOMMENTS, DCM_LT, "", 1, strlen(pszPatComments), {pszPatComments}}, + }; + int nElemRequired = sizeof (aElemRequired) / sizeof(DCM_ELEMENT); + + int iUpdateCount; + cond = DCM_ModifyElements (&m_pObject, aElemRequired, nElemRequired, NULL, 0, &iUpdateCount); + + DCM_ELEMENT elemPixelData = {DCM_PXLPIXELDATA, DCM_OT, "", 1, 0, {NULL}}; + + unsigned long lRealLength = 2 * m_pImageFile->nx() * m_pImageFile->ny(); + + unsigned char* pRawPixels = new unsigned char [lRealLength]; + elemPixelData.length = lRealLength; + elemPixelData.d.ot = pRawPixels; + + ImageFileArray v = m_pImageFile->getArray(); + for (int iy = 0; iy < iNRows; iy++) { + for (int ix = 0; ix < iNCols; ix++) { + unsigned long lBase = (iy * iNRows + ix) * 2; + unsigned int iValue = nearest(dScale * (v[ix][iNRows - 1 - iy] - dMin)); + pRawPixels[lBase] = iValue & 0xFF; + pRawPixels[lBase+1] = (iValue & 0xFF00) >> 8; + } + } + cond = DCM_ModifyElements (&m_pObject, &elemPixelData, 1, NULL, 0, &iUpdateCount); + delete pRawPixels; + + delete pszPatComments; + + if (cond != DCM_NORMAL || iUpdateCount != 1) { + m_bFail = true; + m_strFailMessage = "Error modifying pixel data"; + return false; + } + + return true; +} #endif // HAVE_CTN_DICOM +