32 #include "dcmtk/dcmdata/dcfilefo.h"
33 #include "dcmtk/dcmdata/dcdeftag.h"
34 #include "dcmtk/dcmdata/dcpixel.h"
35 #include "dcmtk/dcmdata/dcpxitem.h"
36 #include "dcmtk/dcmdata/dcpixseq.h"
37 #include "dcmtk/dcmrt/drtimage.h"
59 (
"Compression factor = " + std::to_string(fCompress)
60 +
" has to be a divisor of Number of voxels X = " + std::to_string(
fNoVoxelX)
61 +
" and Y " + std::to_string(
fNoVoxelY)).c_str());
66 for(
int ir = 0; ir <
fNoVoxelY; ir += fCompress ) {
67 for(
int ic = 0; ic <
fNoVoxelX; ic += fCompress ) {
69 int isumrMax =
std::min(ir+fCompress,fNoVoxelY);
70 int isumcMax =
std::min(ic+fCompress,fNoVoxelX);
71 for(
int isumr = ir; isumr < isumrMax; isumr ++ ) {
72 for(
int isumc = ic; isumc < isumcMax; isumc ++ ) {
77 meanHV /= (isumrMax-ir)*(isumcMax-ic);
81 fDensities.push_back(meanDens);
88 fMateIDs.push_back(mateID);
100 for(
int ir = 0; ir <
fNoVoxelY/fCompress; ir++ ) {
101 for(
int ic = 0; ic <
fNoVoxelX/fCompress; ic++ ) {
102 fout << fMateIDs[ic+ir*
fNoVoxelX/fCompress];
103 if( ic !=
fNoVoxelX/fCompress-1) fout <<
" ";
117 for(
int ir = 0; ir <
fNoVoxelY/fCompress; ir++ ) {
118 for(
int ic = 0; ic <
fNoVoxelX/fCompress; ic++ ) {
119 fout << fDensities[ic+ir*
fNoVoxelX/fCompress];
120 if( ic !=
fNoVoxelX/fCompress-1) fout <<
" ";
121 if( copyNo%8 == 7 ) fout <<
G4endl;
124 if( copyNo%8 != 0 ) fout <<
G4endl;
134 if( dfs.size() == 0 )
return;
137 G4int NMAXROI_real = std::log(
INT_MAX)/std::log(NMAXROI);
140 for(
int ir = 0; ir <
fNoVoxelY/fCompress; ir++ ) {
141 for(
int ic = 0; ic <
fNoVoxelX/fCompress; ic++ ) {
143 fStructure.push_back(0);
147 std::set<double> distInters;
151 for(
size_t ii = 0; ii < dfs.size(); ii++ ){
152 std::vector<DicomROI*> rois = dfs[ii]->GetROIs();
153 for(
size_t jj = 0; jj < rois.size(); jj++ ){
155 << rois[jj]->GetNumber() <<
" " << rois[jj]->GetName() <<
G4endl;
156 std::vector<DicomROIContour*> contours = rois[jj]->GetContours();
158 G4int roiID = rois[jj]->GetNumber();
159 for(
size_t kk = 0; kk < contours.size(); kk++ ){
173 std::vector<G4ThreeVector> points = roic->
GetPoints();
175 << points.size() <<
G4endl;
181 for(
size_t ll = 0; ll < points.size(); ll++ ){
184 minYc =
std::min(minYc,points[ll].y());
185 maxYc =
std::max(maxYc,points[ll].y());
187 if( minXc < fMinX || maxXc >
fMaxX || minYc < fMinY || maxYc >
fMaxY ){
189 <<
" maxXc " << maxXc <<
" > " <<
fMaxX
190 <<
" minYc " << minYc <<
" < " <<
fMinY
191 <<
" maxYc " << maxYc <<
" > " << fMaxY <<
G4endl;
195 "Contour limits exceed Z slice extent");
203 <<
" maxXc " << maxXc <<
" > " <<
fMaxX
204 <<
" minYc " << minYc <<
" < " <<
fMinY
205 <<
" maxYc " << maxYc <<
" > " << fMaxY <<
G4endl;
207 G4cout <<
" idMinX " << idMinX
208 <<
" idMaxX " << idMaxX
209 <<
" idMinY " << idMinY
210 <<
" idMaxY " << idMaxY <<
G4endl;
213 for(
int ix = idMinX; ix <= idMaxX; ix++ ) {
214 for(
int iy = idMinY; iy <= idMaxY; iy++ ) {
221 for(
int icx = 0; icx <= 1; icx++ ){
222 for(
int icy = 0; icy <= 1; icy++ ){
228 if( icx == 1 ) v0x = -1.;
231 << iy <<
" + " << icy <<
" CORNER (" << p0x <<
"," << p0y <<
") "
232 <<
" DIR= (" << v0x <<
"," << v0y <<
") " <<
G4endl;
234 for(
int ip = 0; ip < NTRIES; ip++) {
239 <<
" TRYING WITH DIRECTION (" <<
" DIR= (" << v0x <<
","
240 << v0y <<
") " << bOKs <<
G4endl;
241 for(
size_t ll = 0; ll < points.size(); ll++ ){
242 double d0x = points[ll].x() - p0x;
243 double d0y = points[ll].y() - p0y;
244 double w0x = dirs[ll].x();
245 double w0y = dirs[ll].y();
246 double fac1 = w0x*v0y - w0y*v0x;
250 double fac2 = d0x*v0y - d0y*v0x;
251 double fac3 = d0y*w0x - d0x*w0y;
252 double lambdaq = -fac2/fac1;
253 if( lambdaq < 0. || lambdaq >= 1. )
continue;
255 double lambdap = fac3/fac1;
257 distInters.insert(lambdap);
259 <<lambdaq <<
" (" << d0x+p0x+lambdaq*w0x <<
"," << d0y+p0y+lambdaq*w0y
260 <<
") = (" << p0x+lambdap*v0x <<
"," << p0y+lambdap*v0y <<
") "
261 <<
" N " << distInters.size() <<
G4endl;
264 << lambdaq <<
" P " << lambdap <<
G4endl;
267 << d0x+p0x+lambdaq*w0x <<
"," << d0y+p0y+lambdaq*w0y <<
") = ("
268 << p0x+lambdap*v0x <<
"," << p0y+lambdap*v0y <<
") " <<
G4endl;
270 if( distInters.size() % 2 == 1 ) {
271 if( *(distInters.begin() ) > halfDiagonal ) {
275 <<
" N INTERS " << distInters.size() <<
" " << *(distInters.begin())
276 <<
" > " << halfDiagonal << G4endl;
280 if(bOKs == NTRIES ) {
291 int roival = fStructure[ix+iy*
fNoVoxelX/fCompress];
293 if(roival != 0 && roival !=
int(roiID) ) {
294 std::set<G4int> roisVoxel;
295 int nRois = std::log10(roival)/std::log10(NMAXROI)+1;
296 if( nRois > NMAXROI_real ) {
300 G4String(
"Too many ROIs associated to a voxel: \
301 " + std::to_string(nRois) +
" > " + std::to_string(NMAXROI_real) +
" TRY CHAN\
302 GING -NStructureNMaxROI argument to a lower value").c_str());
304 for(
int inr = 0; inr < nRois; inr++ ) {
305 roisVoxel.insert(
int(roival/std::pow(NMAXROI,inr))%NMAXROI );
307 roisVoxel.insert(roiID);
310 for( std::set<G4int>::const_iterator ite = roisVoxel.begin();
311 ite != roisVoxel.end(); ite++, inr++ ) {
312 roival += (*ite)*std::pow(NMAXROI,inr);
314 fStructure[ix+iy*
fNoVoxelX/fCompress] = roival;
316 G4cout <<
" WITH PREVIOUS ROI IN VOXEL " << roival <<
G4endl;
319 fStructure[ix+iy*
fNoVoxelX/fCompress] = roiID;
334 for(
size_t ii = 0; ii < dfs.size(); ii++ ){
335 std::vector<DicomROI*> rois = dfs[ii]->GetROIs();
336 for(
size_t jj = 0; jj < rois.size(); jj++ ){
337 int roi = rois[jj]->GetNumber();
338 std::vector<DicomROIContour*> contours = rois[jj]->GetContours();
339 for(
size_t kk = 0; kk < contours.size(); kk++ ){
345 << kk <<
" INTERS CONTOUR " << roic->
GetZ() <<
" SLICE Z "
348 std::vector<G4ThreeVector> points = roic->
GetPoints();
352 for(
size_t ll = 0; ll < points.size(); ll++ ){
354 <<
" STRUCTURE POINT (" << points[ll].x() <<
"," << points[ll].y() <<
") "
355 <<
" (" << dirs[ll].x() <<
"," << dirs[ll].y() <<
") = " << roi <<
G4endl;
362 for(
int ir = 0; ir <
fNoVoxelY/fCompress; ir++ ) {
363 for(
int ic = 0; ic <
fNoVoxelX/fCompress; ic++ ) {
364 if( fStructure[ic+ir*
fNoVoxelX/fCompress] != 0 ) {
366 <<
" " << ir <<
" STRUCTURE VOXEL (" <<
fMinX +
fVoxelDimX*fCompress * (ic+0.5)
383 if( dfs.size() == 0 )
return;
385 for(
int ir = 0; ir <
fNoVoxelY/fCompress; ir++ ) {
386 for(
int ic = 0; ic <
fNoVoxelX/fCompress; ic++ ) {
387 fout << fStructure[ic+ir*
fNoVoxelX/fCompress];
388 if( ic !=
fNoVoxelX/fCompress-1) fout <<
" ";
G4bool IsMaterialsDensity() const
void DumpStructureIDsToTextFile(std::ofstream &fout)
G4double Hounsfield2density(Uint32 Hval)
void DumpMateIDsToTextFile(std::ofstream &fout)
G4int GetCompression() const
std::vector< DicomFileStructure * > GetStructFiles() const
G4GLOB_DLL std::ostream G4cout
std::vector< G4ThreeVector > GetDirections()
size_t GetMaterialIndexByDensity(G4double density)
std::vector< int > fHounsfieldV
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
DicomFileMgr * theFileMgr
OFString GetGeomType() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
std::vector< G4ThreeVector > GetPoints() const
size_t GetMaterialIndex(G4double Hval)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4int GetStructureNCheck() const
void DumpDensitiesToTextFile(std::ofstream &fout)
G4int GetStructureNMaxROI() const
static DicomFileMgr * GetInstance()
G4GLOB_DLL std::ostream G4cerr