From: jpr Date: Thu, 6 Jul 2006 09:49:15 +0000 (+0000) Subject: Add small utility to convert 'Dense' file format iamges into DICOM images X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=91e8773acafd98b328f10715d8c8964fded57449;p=gdcm.git Add small utility to convert 'Dense' file format iamges into DICOM images (will be usefull to anybody that have such images ...) --- diff --git a/Example/DenseToDicom.cxx b/Example/DenseToDicom.cxx new file mode 100755 index 00000000..3f6234d7 --- /dev/null +++ b/Example/DenseToDicom.cxx @@ -0,0 +1,371 @@ +/*========================================================================= + + Program: gdcm + Module: $RCSfile: DenseToDicom.cxx,v $ + Language: C++ + Date: $Date: 2006/07/06 09:49:15 $ + Version: $Revision: 1.1 $ + + Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de + l'Image). All rights reserved. See Doc/License.txt or + http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ + +#include +#include +//#include + +#include "gdcmFile.h" +#include "gdcmFileHelper.h" +#include "gdcmDebug.h" +#include "gdcmDirList.h" + +#include "gdcmArgMgr.h" + +/** + * \brief + * - explores recursively the given directory + * - examines the ".txt" files + * - Converts the files into 16 bits Dicom Files, + */ + + +void Load(std::ifstream &from, std::string imageName); +//std::ifstream& eatwhite(std::ifstream& is); + +/* +window center (level) and width, are defined in the DICOM +standard very precisely as follows (see PS 3.3 C.11.2.1.2): +> +>"These Attributes are applied according to the following pseudo-code, +>where : +x is the input value, +y is an output value with a range from ymin to ymax, +c is Window Center (0028,1050) +w is Window Width (0028,1051): +> +> if (x <= c - 0.5 - (w-1)/2), then y = ymin +> else if (x > c - 0.5 + (w-1)/2), then y = ymax, +> else y = ((x - (c - 0.5)) / (w-1) + 0.5) * (ymax - ymin)+ ymin + +*/ +/* +From: David Clunie - view profile +Date: Thurs, Jun 1 2006 3:03 pm +Email: David Clunie +Groups: comp.protocols.dicom + +The value of x is the output of the preceding step, the so-called +"modality LUT", which may either be: + +- identity (no rescale values or Modality LUT, or the SOP Class is + PET and rescale values are ignored), in which case x is the stored + pixel value in (7FE0,0010) + +- Rescale Slope and Intercept (as typically used in CT), in which + case x is the value obtained from applying the rescale values to + the stored pixel value + +- an actual LUT, in which case x is the value stored in the LUT + corresponding to the LUT index value that is the stored pixel + value + +The ymin and ymax are intended to represent the output range; for +example, if the hypothetical Presentation LUT step that follows +the VOI LUT (window) stage is an identity operation, then the +ymin and ymax represent P-Values, which might be the range of +digital driving levels for your display (calibrated to the GSDF), +in the 8-bit wide output case ranging from 0 to 255, for example. + +The terms brightness and contrast are not used in radiology imaging +- the window center and width are used instead. +*/ + +int main(int argc, char *argv[]) +{ + START_USAGE(usage) + " \n DenseToDicom :\n ", + " - explores recursively the given directory, ", + " - examines the '.txt' files ", + " - Converts the files into 16 bits Dicom Files, ", + " usage: ", + " Txt2Dcm dirin=rootDirectoryName ", + " [listonly] [verbose] [debug] ", + " ", + " verbose : user wants to run the program in 'verbose mode' ", + " debug : *developer* wants to run the program in 'debug mode' ", + FINISH_USAGE + + // ----- Initialize Arguments Manager ------ + + gdcm::ArgMgr *am = new gdcm::ArgMgr(argc, argv); + + if (argc == 1 || am->ArgMgrDefined("usage")) + { + am->ArgMgrUsage(usage); // Display 'usage' + delete am; + return 0; + } + + const char *dirNamein; + dirNamein = am->ArgMgrGetString("dirin","."); + + if (am->ArgMgrDefined("debug")) + gdcm::Debug::DebugOn(); + + int verbose = am->ArgMgrDefined("verbose"); + int listonly = am->ArgMgrDefined("listonly"); + + // if unused Param we give up + if ( am->ArgMgrPrintUnusedLabels() ) + { + am->ArgMgrUsage(usage); + delete am; + return 0; + } + delete am; // we don't need Argument Manager any longer + + // ----- Begin Processing ----- + + if ( ! gdcm::DirList::IsDirectory(dirNamein) ) + { + std::cout << "KO : [" << dirNamein << "] is not a Directory." << std::endl; + return 0; + } + else + { + std::cout << "OK : [" << dirNamein << "] is a Directory." << std::endl; + } + + std::string strDirNamein(dirNamein); + gdcm::DirList dirList(strDirNamein, true); // get recursively the list of files + + if (listonly) + { + std::cout << "------------List of found files ------------" << std::endl; + dirList.Print(); + std::cout << std::endl; + return 0; + } + + std::string filenameout; + + gdcm::DirListType fileNames; + fileNames = dirList.GetFilenames(); + for (gdcm::DirListType::iterator it = fileNames.begin(); + it != fileNames.end(); + ++it) + { + std::ifstream from( (*it).c_str() ); + if ( !from ) + { + std::cout << "Can't open file" << *it << std::endl; + //return 0; + } + else + { + filenameout = *it + ".dcm"; + std::cout << "Success in open file" << *it << std::endl; + Load(from, filenameout); + //return 0; + } + } +} + + +void Load(std::ifstream &from, std::string imageName) +{ + if (!from) + return; + /* + uint16_t group; + uint16_t elem; + VRKey vr; + TagName vm; + TagName name; +*/ + std::string str1; + int nx, ny; + // Get nx, ny + from >> str1; + from >> str1; + from >> str1; + from >> nx; + from >> str1; + from >> ny; + + std::cout << "nx " << nx << " ny " << ny << std::endl; + + // Skip 2 lines. + std::getline(from, str1); + std::cout << "[" << str1 << "]" << std::endl; + std::getline(from, str1); + std::cout << "[" << str1 << "]" << std::endl; + std::getline(from, str1); + std::cout << "[" << str1 << "]" << std::endl; + + //float *f = new float(nx*ny); + float *f = (float *) malloc(nx*ny*sizeof(float)); + // float mini = FLT_MAX, maxi = FLT_MIN; + float val; + for( int j=0; j> str1; + val = (float)atof(str1.c_str()); + *(f+j*nx+i) = val; + + // In our concern, just *100 all the values is enough! + /* + if (val < mini) + mini = val; + if (val > maxi) + maxi = val; + //std::cout << val << " "; + */ + if(from.eof()) + { + std::cout << "Missing values at [" << j <<"," << i << "]" + << std::endl; + break; + } + l++; + } + //std::cout << std::endl << " line nb : " << j + // << " line length : " << l << std::endl; + + } + + // std::cout << "mini : "<< mini << " maxi : " << maxi << std::endl; +/* +// values are expressed as %. +// It's useless to rescale them as uint16_t : just *100 + uint16_t *img = new uint16_t[ny*nx]; + uint16_t *ptr = img; + float *tmp = f; + float div = maxi-mini; + std::cout << "div : " << div << " mini : " << mini << std::endl; + for( int k=0; kpow(2,maxPow):maxi=pow(2,maxPow) + maxi=math.ceil(maxi)+1 + mini=math.ceil(mini-.5) + slope,intercept=(pow(2,16)-1)/(maxi-mini),-mini + print mini,maxi,"->",slope,intercept + sl=math.ceil(1e6/slope) + ok=1 + while(ok): + val=1/(sl/1e6) + print sl + if val==int(val):ok=0 + else:sl+=1 + print sl/1e6, 1/(sl/1e6) + slope=sl/1e6 + invSlope=1/slope + interc=-math.ceil(intercept*1e6)/1e6 + print invSlope,slope,interc + imageSC=vtkImageShiftScale() + imageSC.SetOutputScalarTypeToUnsignedShort() + imageSC.SetShift(-interc) + imageSC.SetScale(invSlope) + imageSC.SetInput(data) + imageSC.Update() + */ + // gdcm::Debug::DebugOn(); + + std::ostringstream str; + gdcm::File *file; + file = gdcm::File::New(); + + // Set the image size + str.str(""); + str << nx; + file->InsertEntryString(str.str(),0x0028,0x0011,"US"); // Columns + str.str(""); + str << ny; + file->InsertEntryString(str.str(),0x0028,0x0010,"US"); // Rows + // Set the pixel type + // 16; //8, 16, 32 + file->InsertEntryString("16",0x0028,0x0100,"US"); // Bits Allocated + str.str(""); + str << 16; // may be 12 or 16 if componentSize =16 + file->InsertEntryString("16",0x0028,0x0101,"US"); // Bits Stored + + file->InsertEntryString("15",0x0028,0x0102,"US"); // High Bit + // Set the pixel representation // 0/1 + + file->InsertEntryString("0",0x0028,0x0103, "US"); // Pixel Representation + // Set the samples per pixel // 1:Grey level, 3:RGB + file->InsertEntryString("1",0x0028,0x0002, "US"); // Samples per Pixel + +/* + // Set Rescale Intercept + str.str(""); + str << div; + file->InsertEntryString(str.str(),0x0028,0x1052,"DS"); + + // Set Rescale Slope + str.str(""); + str << mini; + file->InsertEntryString(str.str(),0x0028,0x1053,"DS"); +*/ + + file->Print(); + + gdcm::FileHelper *fileH; + fileH = gdcm::FileHelper::New(file); + // cast is just to avoid warnings (*no* conversion) + fileH->SetImageData((uint8_t *)img,nx*ny*sizeof(uint16_t)); + fileH->SetWriteModeToRaw(); + fileH->SetWriteTypeToDcmExplVR(); + + fileH->Write(imageName); + + if( !fileH->Write(imageName)) + std::cout << "Failed for [" << imageName << "]\n" + << " File is unwrittable" << std::endl; + + delete img; + from.close(); +}