File indexing completed on 2025-04-10 08:06:20
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026 #include "DicomFilePET.hh"
0027
0028 #include "DicomFileStructure.hh"
0029 #include "DicomROI.hh"
0030 #include "dcmtk/dcmdata/dcdeftag.h"
0031 #include "dcmtk/dcmdata/dcfilefo.h"
0032 #include "dcmtk/dcmdata/dcpixel.h"
0033 #include "dcmtk/dcmdata/dcpixseq.h"
0034 #include "dcmtk/dcmdata/dcpxitem.h"
0035 #include "dcmtk/dcmrt/drtimage.h"
0036
0037 #include "G4GeometryTolerance.hh"
0038
0039 #include <set>
0040
0041
0042 DicomFilePET::DicomFilePET() {}
0043
0044
0045 DicomFilePET::DicomFilePET(DcmDataset* dset) : DicomVFileImage(dset) {}
0046
0047
0048 void DicomFilePET::BuildActivities()
0049 {
0050 G4int fCompress = theFileMgr->GetCompression();
0051 if (fNoVoxelsX % fCompress != 0 || fNoVoxelsY % fCompress != 0) {
0052 G4Exception("DicompFileMgr.:BuildMaterials", "DFC004", FatalException,
0053 ("Compression factor = " + std::to_string(fCompress)
0054 + " has to be a divisor of Number of voxels X = " + std::to_string(fNoVoxelsX)
0055 + " and Y " + std::to_string(fNoVoxelsY))
0056 .c_str());
0057 }
0058
0059
0060 double meanHV = 0.;
0061 for (int ir = 0; ir < fNoVoxelsY; ir += fCompress) {
0062 for (int ic = 0; ic < fNoVoxelsX; ic += fCompress) {
0063 meanHV = 0.;
0064 int isumrMax = std::min(ir + fCompress, fNoVoxelsY);
0065 int isumcMax = std::min(ic + fCompress, fNoVoxelsX);
0066 for (int isumr = ir; isumr < isumrMax; isumr++) {
0067 for (int isumc = ic; isumc < isumcMax; isumc++) {
0068 meanHV += fHounsfieldV[isumc + isumr * fNoVoxelsX];
0069 }
0070 }
0071 meanHV /= (isumrMax - ir) * (isumcMax - ic);
0072 fActivities.push_back(meanHV);
0073 }
0074 }
0075 }
0076
0077
0078 void DicomFilePET::DumpActivitiesToTextFile(std::ofstream& fout)
0079 {
0080 G4int fCompress = theFileMgr->GetCompression();
0081 if (DicomFileMgr::verbose >= warningVerb)
0082 G4cout << fLocation << " DumpDensitiesToTextFile " << fFileName << " " << fActivities.size()
0083 << G4endl;
0084
0085 G4int copyNo = 0;
0086 for (int ir = 0; ir < fNoVoxelsY / fCompress; ir++) {
0087 for (int ic = 0; ic < fNoVoxelsX / fCompress; ic++) {
0088 fout << fActivities[ic + ir * fNoVoxelsX / fCompress];
0089 if (ic != fNoVoxelsX / fCompress - 1) fout << " ";
0090 if (copyNo % 8 == 7) fout << G4endl;
0091 copyNo++;
0092 }
0093 if (copyNo % 8 != 0) fout << G4endl;
0094 }
0095 }