Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /geant4/examples/extended/medical/DICOM/G4DicomReader/src/DicomFileCT.cc was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 #include "DicomFileCT.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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0042 DicomFileCT::DicomFileCT() {}
0043 
0044 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0045 DicomFileCT::DicomFileCT(DcmDataset* dset) : DicomVFileImage(dset) {}
0046 
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048 void DicomFileCT::BuildMaterials()
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   //  if( DicomVerb(debugVerb) ) G4cout << " BuildMaterials " << fFileName << G4endl;
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           // G4cout << isumr << " " << isumc << " GET mean " << meanHV << G4endl;
0070         }
0071       }
0072       meanHV /= (isumrMax - ir) * (isumcMax - ic);
0073       G4double meanDens = theFileMgr->Hounsfield2density(meanHV);
0074       //      G4cout << ir << " " << ic << " FINAL mean " << meanDens << G4endl;
0075 
0076       fDensities.push_back(meanDens);
0077       size_t mateID;
0078       if (theFileMgr->IsMaterialsDensity()) {
0079         mateID = theFileMgr->GetMaterialIndexByDensity(meanDens);
0080       }
0081       else {
0082         mateID = theFileMgr->GetMaterialIndex(meanHV);
0083       }
0084       fMateIDs.push_back(mateID);
0085     }
0086   }
0087 }
0088 
0089 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0090 void DicomFileCT::DumpMateIDsToTextFile(std::ofstream& fout)
0091 {
0092   G4int fCompress = theFileMgr->GetCompression();
0093   if (DicomFileMgr::verbose >= warningVerb)
0094     G4cout << fLocation << " DumpMateIDsToTextFile " << fFileName << " " << fMateIDs.size()
0095            << G4endl;
0096   for (int ir = 0; ir < fNoVoxelsY / fCompress; ir++) {
0097     for (int ic = 0; ic < fNoVoxelsX / fCompress; ic++) {
0098       fout << fMateIDs[ic + ir * fNoVoxelsX / fCompress];
0099       if (ic != fNoVoxelsX / fCompress - 1) fout << " ";
0100     }
0101     fout << G4endl;
0102   }
0103 }
0104 
0105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0106 void DicomFileCT::DumpDensitiesToTextFile(std::ofstream& fout)
0107 {
0108   G4int fCompress = theFileMgr->GetCompression();
0109   if (DicomFileMgr::verbose >= warningVerb)
0110     G4cout << fLocation << " DumpDensitiesToTextFile " << fFileName << " " << fDensities.size()
0111            << G4endl;
0112 
0113   G4int copyNo = 0;
0114   for (int ir = 0; ir < fNoVoxelsY / fCompress; ir++) {
0115     for (int ic = 0; ic < fNoVoxelsX / fCompress; ic++) {
0116       fout << fDensities[ic + ir * fNoVoxelsX / fCompress];
0117       if (ic != fNoVoxelsX / fCompress - 1) fout << " ";
0118       if (copyNo % 8 == 7) fout << G4endl;
0119       copyNo++;
0120     }
0121     if (copyNo % 8 != 0) fout << G4endl;
0122   }
0123 }
0124 
0125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0126 void DicomFileCT::BuildStructureIDs()
0127 {
0128   G4int fCompress = theFileMgr->GetCompression();
0129   std::vector<DicomFileStructure*> dfs = theFileMgr->GetStructFiles();
0130   if (dfs.size() == 0) return;
0131 
0132   G4int NMAXROI = DicomFileMgr::GetInstance()->GetStructureNMaxROI();
0133   G4int NMAXROI_real = std::log(INT_MAX) / std::log(NMAXROI);
0134 
0135   //  fStructure = new G4int(fNoVoxelsX*fNoVoxelsY);
0136   for (int ir = 0; ir < fNoVoxelsY / fCompress; ir++) {
0137     for (int ic = 0; ic < fNoVoxelsX / fCompress; ic++) {
0138       //      fStructure[ic+ir*fNoVoxelsX] = 0;
0139       fStructure.push_back(0);
0140     }
0141   }
0142 
0143   std::set<double> distInters;
0144 
0145   //  std::fill_n(fStructure,fNoVoxelsX*fNoVoxelsY,0);
0146   //
0147   for (size_t ii = 0; ii < dfs.size(); ii++) {
0148     std::vector<DicomROI*> rois = dfs[ii]->GetROIs();
0149     for (size_t jj = 0; jj < rois.size(); jj++) {
0150       if (DicomFileMgr::verbose >= debugVerb)
0151         std::cout << " BuildStructureIDs checking ROI " << rois[jj]->GetNumber() << " "
0152                   << rois[jj]->GetName() << G4endl;
0153       std::vector<DicomROIContour*> contours = rois[jj]->GetContours();
0154       //      G4int roi = jj+1;
0155       G4int roiID = rois[jj]->GetNumber();
0156       for (size_t kk = 0; kk < contours.size(); kk++) {
0157         // check contour corresponds to this CT slice
0158         DicomROIContour* roic = contours[kk];
0159         // if( DicomVerb(-debugVerb) ) G4cout << jj << " " << kk << " INTERS CONTOUR " << roic
0160         //                << " " << fLocation << G4endl;
0161         if (DicomFileMgr::verbose >= debugVerb)
0162           G4cout << " " << roiID << " " << kk << " INTERS CONTOUR " << roic->GetZ() << " SLICE Z "
0163                  << fMinZ << " " << fMaxZ << G4endl;
0164         // Check Contour correspond to slice
0165 
0166         if (roic->GetZ() > fMaxZ || roic->GetZ() < fMinZ) continue;
0167         if (DicomFileMgr::verbose >= debugVerb)
0168           G4cout << " INTERS CONTOUR WITH Z SLIZE " << fMinZ << " < " << roic->GetZ() << " < "
0169                  << fMaxZ << G4endl;
0170         if (roic->GetGeomType() == "CLOSED_PLANAR") {
0171           // Get min and max X and Y, and corresponding slice indexes
0172           std::vector<G4ThreeVector> points = roic->GetPoints();
0173           if (DicomFileMgr::verbose >= debugVerb)
0174             G4cout << jj << " " << kk << " NPOINTS " << points.size() << G4endl;
0175           std::vector<G4ThreeVector> dirs = roic->GetDirections();
0176           double minXc = DBL_MAX;
0177           double maxXc = -DBL_MAX;
0178           double minYc = DBL_MAX;
0179           double maxYc = -DBL_MAX;
0180           for (size_t ll = 0; ll < points.size(); ll++) {
0181             minXc = std::min(minXc, points[ll].x());
0182             maxXc = std::max(maxXc, points[ll].x());
0183             minYc = std::min(minYc, points[ll].y());
0184             maxYc = std::max(maxYc, points[ll].y());
0185           }
0186           if (minXc < fMinX || maxXc > fMaxX || minYc < fMinY || maxYc > fMaxY) {
0187             G4cerr << " minXc " << minXc << " < " << fMinX << " maxXc " << maxXc << " > " << fMaxX
0188                    << " minYc " << minYc << " < " << fMinY << " maxYc " << maxYc << " > " << fMaxY
0189                    << G4endl;
0190             G4Exception("DicomFileCT::BuildStructureIDs", "DFCT001", JustWarning,
0191                         "Contour limits exceed Z slice extent");
0192           }
0193           int idMinX = std::max(0, int((minXc - fMinX) / fVoxelDimX / fCompress));
0194           int idMaxX =
0195             std::min(fNoVoxelsX / fCompress - 1, int((maxXc - fMinX) / fVoxelDimX / fCompress + 1));
0196           int idMinY = std::max(0, int((minYc - fMinY) / fVoxelDimY / fCompress));
0197           int idMaxY =
0198             std::min(fNoVoxelsY / fCompress - 1, int((maxYc - fMinY) / fVoxelDimY / fCompress + 1));
0199           if (DicomFileMgr::verbose >= debugVerb)
0200             G4cout << " minXc " << minXc << " < " << fMinX << " maxXc " << maxXc << " > " << fMaxX
0201                    << " minYc " << minYc << " < " << fMinY << " maxYc " << maxYc << " > " << fMaxY
0202                    << G4endl;
0203           if (DicomFileMgr::verbose >= debugVerb)
0204             G4cout << " idMinX " << idMinX << " idMaxX " << idMaxX << " idMinY " << idMinY
0205                    << " idMaxY " << idMaxY << G4endl;
0206           // for each voxel: build 4 lines from the corner towards the center
0207           //  and check how many contour segments it crosses, and the minimum distance to a segment
0208           for (int ix = idMinX; ix <= idMaxX; ix++) {
0209             for (int iy = idMinY; iy <= idMaxY; iy++) {
0210               G4bool bOK = false;
0211               G4int bOKs;
0212               if (DicomFileMgr::verbose >= debugVerb)
0213                 G4cout << "@@@@@ TRYING POINT (" << fMinX + fVoxelDimX * fCompress * (ix + 0.5)
0214                        << "," << fMinY + fVoxelDimY * fCompress * (iy + 0.5) << ")" << G4endl;
0215               // four corners
0216               for (int icx = 0; icx <= 1; icx++) {
0217                 for (int icy = 0; icy <= 1; icy++) {
0218                   bOKs = 0;
0219                   if (bOK) continue;
0220                   double p0x = fMinX + fVoxelDimX * fCompress * (ix + icx);
0221                   double p0y = fMinY + fVoxelDimY * fCompress * (iy + icy);
0222                   double v0x = 1.;
0223                   if (icx == 1) v0x = -1.;
0224                   double v0y = 0.99 * fVoxelDimY / fVoxelDimX * std::pow(-1., icy);
0225                   if (DicomFileMgr::verbose >= testVerb)
0226                     G4cout << ix << " + " << icx << " " << iy << " + " << icy << " CORNER (" << p0x
0227                            << "," << p0y << ") "
0228                            << " DIR= (" << v0x << "," << v0y << ") " << G4endl;
0229                   int NTRIES = theFileMgr->GetStructureNCheck();
0230                   for (int ip = 0; ip < NTRIES; ip++) {
0231                     distInters.clear();
0232                     v0y -= ip * 0.1 * v0y;
0233                     G4double halfDiagonal = 0.5 * fVoxelDimX * fCompress;
0234                     if (DicomFileMgr::verbose >= testVerb)
0235                       G4cout << ip << " TRYING WITH DIRECTION ("
0236                              << " DIR= (" << v0x << "," << v0y << ") " << bOKs << G4endl;
0237                     for (size_t ll = 0; ll < points.size(); ll++) {
0238                       double d0x = points[ll].x() - p0x;
0239                       double d0y = points[ll].y() - p0y;
0240                       double w0x = dirs[ll].x();
0241                       double w0y = dirs[ll].y();
0242                       double fac1 = w0x * v0y - w0y * v0x;
0243                       if (fac1 == 0) {  // parallel lines
0244                         continue;
0245                       }
0246                       double fac2 = d0x * v0y - d0y * v0x;
0247                       double fac3 = d0y * w0x - d0x * w0y;
0248                       double lambdaq = -fac2 / fac1;
0249                       if (lambdaq < 0. || lambdaq >= 1.) continue;
0250                       // intersection further than segment length
0251                       double lambdap = fac3 / fac1;
0252                       if (lambdap > 0.) {
0253                         distInters.insert(lambdap);
0254                         if (DicomFileMgr::verbose >= testVerb)
0255                           G4cout << " !! GOOD INTERS " << lambdaq << "  ("
0256                                  << d0x + p0x + lambdaq * w0x << "," << d0y + p0y + lambdaq * w0y
0257                                  << ")  =  (" << p0x + lambdap * v0x << "," << p0y + lambdap * v0y
0258                                  << ") "
0259                                  << " N " << distInters.size() << G4endl;
0260                       }
0261                       if (DicomFileMgr::verbose >= testVerb)
0262                         G4cout << " INTERS LAMBDAQ " << lambdaq << " P " << lambdap << G4endl;
0263 
0264                       if (DicomFileMgr::verbose >= debugVerb)
0265                         G4cout << " INTERS POINT (" << d0x + p0x + lambdaq * w0x << ","
0266                                << d0y + p0y + lambdaq * w0y << ")  =  (" << p0x + lambdap * v0x
0267                                << "," << p0y + lambdap * v0y << ") " << G4endl;
0268                     }
0269                     if (distInters.size() % 2 == 1) {
0270                       if (*(distInters.begin()) > halfDiagonal) {
0271                         //                      bOK = true;
0272                         bOKs += 1;
0273                         if (DicomFileMgr::verbose >= debugVerb)
0274                           G4cout << " OK= " << bOKs << " N INTERS " << distInters.size() << " "
0275                                  << *(distInters.begin()) << " > " << halfDiagonal << G4endl;
0276                       }
0277                     }
0278                   }
0279                   if (bOKs == NTRIES) {
0280                     bOK = true;
0281                     if (DicomFileMgr::verbose >= debugVerb)
0282                       G4cout << "@@@@@ POINT OK (" << fMinX + fVoxelDimX * fCompress * (ix + 0.5)
0283                              << "," << fMinY + fVoxelDimY * fCompress * (iy + 0.5) << ")" << G4endl;
0284                   }
0285                 }
0286               }  // loop to four corners
0287               if (bOK) {
0288                 // extract previous ROI value
0289                 int roival = fStructure[ix + iy * fNoVoxelsX / fCompress];
0290                 //                roival = 2 + NMAXROI*3 + NMAXROI*NMAXROI*15;
0291                 if (roival != 0 && roival != int(roiID)) {
0292                   std::set<G4int> roisVoxel;
0293                   int nRois = std::log10(roival) / std::log10(NMAXROI) + 1;
0294                   if (nRois > NMAXROI_real) {  // another one cannot be added
0295                     G4Exception("DicomFileCT:BuildStructureIDs", "DFC0004", FatalException,
0296                                 G4String("Too many ROIs associated to a voxel: \
0297 " + std::to_string(nRois) + " > " + std::to_string(NMAXROI_real)
0298                                          + " TRY CHAN\
0299 GING -NStructureNMaxROI argument to a lower value")
0300                                   .c_str());
0301                   }
0302                   for (int inr = 0; inr < nRois; inr++) {
0303                     roisVoxel.insert(int(roival / std::pow(NMAXROI, inr)) % NMAXROI);
0304                   }
0305                   roisVoxel.insert(roiID);
0306                   roival = 0;
0307                   size_t inr = 0;
0308                   for (std::set<G4int>::const_iterator ite = roisVoxel.begin();
0309                        ite != roisVoxel.end(); ite++, inr++)
0310                   {
0311                     roival += (*ite) * std::pow(NMAXROI, inr);
0312                   }
0313                   fStructure[ix + iy * fNoVoxelsX / fCompress] = roival;
0314                   if (DicomFileMgr::verbose >= testVerb) {
0315                     G4cout << " WITH PREVIOUS ROI IN VOXEL " << roival << G4endl;
0316                   }
0317                 }
0318                 else {
0319                   fStructure[ix + iy * fNoVoxelsX / fCompress] = roiID;
0320                 }
0321               }
0322             }
0323           }
0324         }
0325       }
0326     }
0327   }
0328 
0329   if (DicomFileMgr::verbose >= 1) {
0330     //@@@@ PRINT structures
0331     //@@@ PRINT points of structures in this Z slice
0332     if (DicomFileMgr::verbose >= 0) G4cout << " STRUCTURES FOR SLICE " << fLocation << G4endl;
0333     for (size_t ii = 0; ii < dfs.size(); ii++) {
0334       std::vector<DicomROI*> rois = dfs[ii]->GetROIs();
0335       for (size_t jj = 0; jj < rois.size(); jj++) {
0336         int roi = rois[jj]->GetNumber();
0337         std::vector<DicomROIContour*> contours = rois[jj]->GetContours();
0338         for (size_t kk = 0; kk < contours.size(); kk++) {
0339           DicomROIContour* roic = contours[kk];
0340           // check contour corresponds to this CT slice
0341           if (roic->GetZ() > fMaxZ || roic->GetZ() < fMinZ) continue;
0342           if (roic->GetGeomType() == "CLOSED_PLANAR") {
0343             if (DicomFileMgr::verbose >= testVerb)
0344               G4cout << " PRINTING CONTOUR? " << roi << " " << kk << " INTERS CONTOUR "
0345                      << roic->GetZ() << " SLICE Z " << fMinZ << " " << fMaxZ << G4endl;
0346             if (roic->GetZ() > fMaxZ || roic->GetZ() < fMinZ) continue;
0347             std::vector<G4ThreeVector> points = roic->GetPoints();
0348             std::vector<G4ThreeVector> dirs = roic->GetDirections();
0349             if (DicomFileMgr::verbose >= 0)
0350               std::cout << " STRUCTURE Z " << roic->GetZ() << std::endl;
0351             for (size_t ll = 0; ll < points.size(); ll++) {
0352               if (DicomFileMgr::verbose >= 0)
0353                 G4cout << roi << " : " << ll << " STRUCTURE POINT (" << points[ll].x() << ","
0354                        << points[ll].y() << ") "
0355                        << " (" << dirs[ll].x() << "," << dirs[ll].y() << ") = " << roi << G4endl;
0356             }
0357           }
0358         }
0359       }
0360     }
0361     //@@@ PRINT points in slice inside structure
0362     for (int ir = 0; ir < fNoVoxelsY / fCompress; ir++) {
0363       for (int ic = 0; ic < fNoVoxelsX / fCompress; ic++) {
0364         if (fStructure[ic + ir * fNoVoxelsX / fCompress] != 0) {
0365           if (DicomFileMgr::verbose >= 0)
0366             G4cout << ic + ir * fNoVoxelsX / fCompress << " = " << ic << " " << ir
0367                    << " STRUCTURE VOXEL (" << fMinX + fVoxelDimX * fCompress * (ic + 0.5) << ","
0368                    << fMinY + fVoxelDimY * fCompress * (ir + 0.5)
0369                    << ") = " << fStructure[ic + ir * fNoVoxelsX / fCompress] << G4endl;
0370         }
0371       }
0372     }
0373   }
0374 }
0375 
0376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0377 void DicomFileCT::DumpStructureIDsToTextFile(std::ofstream& fout)
0378 {
0379   G4int fCompress = theFileMgr->GetCompression();
0380   if (DicomFileMgr::verbose >= 0)
0381     G4cout << fLocation << " DumpStructureIDsToTextFile " << fFileName << " " << fStructure.size()
0382            << G4endl;
0383   std::vector<DicomFileStructure*> dfs = theFileMgr->GetStructFiles();
0384   if (dfs.size() == 0) return;
0385 
0386   for (int ir = 0; ir < fNoVoxelsY / fCompress; ir++) {
0387     for (int ic = 0; ic < fNoVoxelsX / fCompress; ic++) {
0388       fout << fStructure[ic + ir * fNoVoxelsX / fCompress];
0389       if (ic != fNoVoxelsX / fCompress - 1) fout << " ";
0390     }
0391     fout << G4endl;
0392   }
0393 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!