File indexing completed on 2025-02-23 09:21:45
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
0027
0028
0029
0030
0031 #include "DicomIntersectVolume.hh"
0032
0033 #include "G4LogicalVolume.hh"
0034 #include "G4LogicalVolumeStore.hh"
0035 #include "G4Material.hh"
0036 #include "G4PVParameterised.hh"
0037 #include "G4PhantomParameterisation.hh"
0038 #include "G4PhysicalVolumeStore.hh"
0039 #include "G4SystemOfUnits.hh"
0040 #include "G4UIcmdWithAString.hh"
0041 #include "G4UIcommand.hh"
0042 #include "G4VPhysicalVolume.hh"
0043 #include "G4VSolid.hh"
0044 #include "G4tgbVolume.hh"
0045 #include "G4tgrSolid.hh"
0046
0047
0048 DicomIntersectVolume::DicomIntersectVolume()
0049 : G4UImessenger(), fG4VolumeCmd(0), fSolid(0), fPhantomMinusCorner(), fVoxelIsInside(0)
0050 {
0051 fUserVolumeCmd = new G4UIcmdWithAString("/dicom/intersectWithUserVolume", this);
0052 fUserVolumeCmd->SetGuidance(
0053 "Intersects a phantom with a user-defined volume "
0054 "and outputs the voxels that are totally inside the intersection as"
0055 " a new phantom file. It must have the parameters: POS_X POS_Y POS_Z "
0056 "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_1 (SOLID_PARAM_2 ...)");
0057 fUserVolumeCmd->SetParameterName("choice", true);
0058 fUserVolumeCmd->AvailableForStates(G4State_Idle);
0059
0060 fG4VolumeCmd = new G4UIcmdWithAString("/dicom/intersectWithUserVolume", this);
0061 fG4VolumeCmd->SetGuidance(
0062 "Intersects a phantom with a user-defined volume"
0063 " and outputs the voxels that are totally inside the intersection as "
0064 " a new phantom file. It must have the parameters: POS_X POS_Y POS_Z "
0065 "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_1 (SOLID_PARAM_2 ...)");
0066 fG4VolumeCmd->SetParameterName("choice", true);
0067 fG4VolumeCmd->AvailableForStates(G4State_Idle);
0068 }
0069
0070
0071 DicomIntersectVolume::~DicomIntersectVolume()
0072 {
0073 if (fUserVolumeCmd) delete fUserVolumeCmd;
0074 if (fG4VolumeCmd) delete fG4VolumeCmd;
0075 }
0076
0077
0078 void DicomIntersectVolume::SetNewValue(G4UIcommand* command, G4String newValues)
0079 {
0080 G4AffineTransform theVolumeTransform;
0081
0082 if (command == fUserVolumeCmd) {
0083 std::vector<G4String> params = GetWordsInString(newValues);
0084 if (params.size() < 8) {
0085 G4Exception("DicomIntersectVolume::SetNewValue",
0086 " There must be at least 8 parameter: SOLID_TYPE POS_X POS_Y POS_Z "
0087 "ANG_X ANG_Y ANG_Z SOLID_PARAM_1 (SOLID_PARAM_2 ...)",
0088 FatalErrorInArgument,
0089 G4String("Number of parameters given = "
0090 + G4UIcommand::ConvertToString(G4int(params.size())))
0091 .c_str());
0092 }
0093
0094
0095 BuildUserSolid(params);
0096
0097
0098 G4ThreeVector pos = G4ThreeVector(G4UIcommand::ConvertToDouble(params[0]),
0099 G4UIcommand::ConvertToDouble(params[1]),
0100 G4UIcommand::ConvertToDouble(params[2]));
0101 G4RotationMatrix* rotmat = new G4RotationMatrix;
0102 std::vector<G4double> angles;
0103 rotmat->rotateX(G4UIcommand::ConvertToDouble(params[3]));
0104 rotmat->rotateY(G4UIcommand::ConvertToDouble(params[4]));
0105 rotmat->rotateY(G4UIcommand::ConvertToDouble(params[5]));
0106 theVolumeTransform = G4AffineTransform(rotmat, pos).Invert();
0107 }
0108 else if (command == fG4VolumeCmd) {
0109 std::vector<G4String> params = GetWordsInString(newValues);
0110 if (params.size() != 1)
0111 G4Exception("DicomIntersectVolume::SetNewValue", "", FatalErrorInArgument,
0112 G4String("Command: " + command->GetCommandPath() + "/" + command->GetCommandName()
0113 + " " + newValues + " needs 1 argument: VOLUME_NAME")
0114 .c_str());
0115
0116
0117 BuildG4Solid(params);
0118
0119
0120 G4VPhysicalVolume* pv = GetPhysicalVolumes(params[0], 1, 1)[0];
0121
0122 theVolumeTransform = G4AffineTransform(pv->GetFrameRotation(), pv->GetFrameTranslation());
0123 }
0124
0125
0126 G4PhantomParameterisation* thePhantomParam = GetPhantomParam(true);
0127
0128 G4RotationMatrix* rotph = new G4RotationMatrix();
0129
0130 G4AffineTransform thePhantomTransform(rotph, G4ThreeVector());
0131
0132
0133 G4AffineTransform theTransform = theVolumeTransform * thePhantomTransform;
0134
0135 G4ThreeVector axisX(1., 0., 0.);
0136 G4ThreeVector axisY(0., 1., 0.);
0137 G4ThreeVector axisZ(0., 0., 1.);
0138 theTransform.ApplyAxisTransform(axisX);
0139 theTransform.ApplyAxisTransform(axisY);
0140 theTransform.ApplyAxisTransform(axisZ);
0141
0142
0143 G4String thePhantomFileName = "phantom.g4pdcm";
0144 fout.open(thePhantomFileName);
0145 std::vector<G4Material*> materials = thePhantomParam->GetMaterials();
0146 fout << materials.size() << G4endl;
0147 for (unsigned int ii = 0; ii < materials.size(); ++ii) {
0148 fout << ii << " " << materials[ii]->GetName() << G4endl;
0149 }
0150
0151
0152 G4int nx = G4int(thePhantomParam->GetNoVoxelsX());
0153 G4int ny = G4int(thePhantomParam->GetNoVoxelsY());
0154 G4int nz = G4int(thePhantomParam->GetNoVoxelsZ());
0155 G4int nxy = nx * ny;
0156 fVoxelIsInside = new G4bool[nx * ny * nz];
0157 G4double voxelHalfWidthX = thePhantomParam->GetVoxelHalfX();
0158 G4double voxelHalfWidthY = thePhantomParam->GetVoxelHalfY();
0159 G4double voxelHalfWidthZ = thePhantomParam->GetVoxelHalfZ();
0160
0161
0162 fout << nx << " " << ny << " " << nz << G4endl;
0163 fout << -voxelHalfWidthX * nx + thePhantomTransform.NetTranslation().x() << " "
0164 << voxelHalfWidthX * nx + thePhantomTransform.NetTranslation().x() << G4endl;
0165 fout << -voxelHalfWidthY * ny + thePhantomTransform.NetTranslation().y() << " "
0166 << voxelHalfWidthY * ny + thePhantomTransform.NetTranslation().y() << G4endl;
0167 fout << -voxelHalfWidthZ * nz + thePhantomTransform.NetTranslation().z() << " "
0168 << voxelHalfWidthZ * nz + thePhantomTransform.NetTranslation().z() << G4endl;
0169
0170 for (G4int iz = 0; iz < nz; ++iz) {
0171 for (G4int iy = 0; iy < ny; ++iy) {
0172 G4bool bPrevVoxelInside = true;
0173 G4bool b1VoxelFoundInside = false;
0174 G4int firstVoxel = -1;
0175 G4int lastVoxel = -1;
0176 for (G4int ix = 0; ix < nx; ++ix) {
0177 G4ThreeVector voxelCentre((-nx + ix * 2 + 1) * voxelHalfWidthX,
0178 (-ny + iy * 2 + 1) * voxelHalfWidthY,
0179 (-nz + iz * 2 + 1) * voxelHalfWidthZ);
0180 theTransform.ApplyPointTransform(voxelCentre);
0181 G4bool bVoxelIsInside = true;
0182 for (G4int ivx = -1; ivx <= 1; ivx += 2) {
0183 for (G4int ivy = -1; ivy <= 1; ivy += 2) {
0184 for (G4int ivz = -1; ivz <= 1; ivz += 2) {
0185 G4ThreeVector voxelPoint = voxelCentre + ivx * voxelHalfWidthX * axisX
0186 + ivy * voxelHalfWidthY * axisY
0187 + ivz * voxelHalfWidthZ * axisZ;
0188 if (fSolid->Inside(voxelPoint) == kOutside) {
0189 bVoxelIsInside = false;
0190 break;
0191 }
0192 else {
0193 }
0194 }
0195 if (!bVoxelIsInside) break;
0196 }
0197 if (!bVoxelIsInside) break;
0198 }
0199
0200 G4int copyNo = ix + nx * iy + nxy * iz;
0201 if (bVoxelIsInside) {
0202 if (!bPrevVoxelInside) {
0203 G4Exception("DicomIntersectVolume::SetNewValue", "", FatalException,
0204 "Volume cannot intersect phantom in discontiguous voxels, "
0205 "please use other voxel");
0206 }
0207 if (!b1VoxelFoundInside) {
0208 firstVoxel = ix;
0209 b1VoxelFoundInside = true;
0210 }
0211 lastVoxel = ix;
0212 fVoxelIsInside[copyNo] = true;
0213 }
0214 else {
0215 fVoxelIsInside[copyNo] = false;
0216 }
0217 }
0218 fout << firstVoxel << " " << lastVoxel << G4endl;
0219 }
0220 }
0221
0222
0223 for (G4int iz = 0; iz < nz; ++iz) {
0224 for (G4int iy = 0; iy < ny; ++iy) {
0225 G4bool b1xFound = false;
0226 for (G4int ix = 0; ix < nx; ++ix) {
0227 size_t copyNo = ix + ny * iy + nxy * iz;
0228 if (fVoxelIsInside[copyNo]) {
0229 fout << thePhantomParam->GetMaterialIndex(copyNo) << " ";
0230 b1xFound = true;
0231 }
0232 }
0233 if (b1xFound) fout << G4endl;
0234 }
0235 }
0236
0237
0238 for (G4int iz = 0; iz < nz; ++iz) {
0239 for (G4int iy = 0; iy < ny; ++iy) {
0240 G4bool b1xFound = false;
0241 for (G4int ix = 0; ix < nx; ++ix) {
0242 size_t copyNo = ix + ny * iy + nxy * iz;
0243 if (fVoxelIsInside[copyNo]) {
0244 fout << thePhantomParam->GetMaterial(copyNo)->GetDensity() / g * cm3 << " ";
0245 b1xFound = true;
0246 }
0247 }
0248 if (b1xFound) fout << G4endl;
0249 }
0250 }
0251 }
0252
0253
0254 void DicomIntersectVolume::BuildUserSolid(std::vector<G4String> params)
0255 {
0256 for (G4int ii = 0; ii < 6; ++ii)
0257 params.erase(params.begin());
0258
0259 params.insert(params.begin(), ":SOLID");
0260 params.insert(params.begin(), params[1]);
0261 G4tgrSolid* tgrSolid = new G4tgrSolid(params);
0262 G4tgbVolume* tgbVolume = new G4tgbVolume();
0263 fSolid = tgbVolume->FindOrConstructG4Solid(tgrSolid);
0264 }
0265
0266
0267 void DicomIntersectVolume::BuildG4Solid(std::vector<G4String> params)
0268 {
0269 fSolid = GetLogicalVolumes(params[0], 1, 1)[0]->GetSolid();
0270 }
0271
0272
0273 G4PhantomParameterisation* DicomIntersectVolume::GetPhantomParam(G4bool bMustExist)
0274 {
0275 G4PhantomParameterisation* paramreg = 0;
0276
0277 G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
0278 for (auto cite = pvs->cbegin(); cite != pvs->cend(); ++cite) {
0279 if (IsPhantomVolume(*cite)) {
0280 const G4PVParameterised* pvparam = static_cast<const G4PVParameterised*>(*cite);
0281 G4VPVParameterisation* param = pvparam->GetParameterisation();
0282
0283
0284
0285
0286 paramreg = static_cast<G4PhantomParameterisation*>(param);
0287 }
0288 }
0289
0290 if (!paramreg && bMustExist)
0291 G4Exception("DicomIntersectVolume::GetPhantomParam", "", FatalErrorInArgument,
0292 " No G4PhantomParameterisation found ");
0293
0294 return paramreg;
0295 }
0296
0297
0298 std::vector<G4VPhysicalVolume*> DicomIntersectVolume::GetPhysicalVolumes(const G4String& name,
0299 G4bool exists, G4int nVols)
0300 {
0301 std::vector<G4VPhysicalVolume*> vvolu;
0302 std::string::size_type ial = name.rfind(":");
0303 G4String volname = "";
0304 G4int volcopy = 0;
0305 if (ial != std::string::npos) {
0306 std::string::size_type ial2 = name.rfind(":", ial - 1);
0307 if (ial2 != std::string::npos) {
0308 G4Exception("DicomIntersectVolume::GetPhysicalVolumes", "", FatalErrorInArgument,
0309 G4String("Name corresponds to a touchable " + name).c_str());
0310 }
0311 else {
0312 volname = name.substr(0, ial);
0313 volcopy = G4UIcommand::ConvertToInt(name.substr(ial + 1, name.length()).c_str());
0314 }
0315 }
0316 else {
0317 volname = name;
0318 volcopy = -1;
0319 }
0320
0321 G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
0322 for (auto citepv = pvs->cbegin(); citepv != pvs->cend(); ++citepv) {
0323 if (volname == (*citepv)->GetName() && ((*citepv)->GetCopyNo() == volcopy || -1 == volcopy)) {
0324 vvolu.push_back(*citepv);
0325 }
0326 }
0327
0328
0329 if (vvolu.size() == 0) {
0330 if (exists) {
0331 G4Exception(" DicomIntersectVolume::GetPhysicalVolumes", "", FatalErrorInArgument,
0332 G4String("No physical volume found with name " + name).c_str());
0333 }
0334 else {
0335 G4cerr << "!!WARNING: DicomIntersectVolume::GetPhysicalVolumes: "
0336 << " no physical volume found with name " << name << G4endl;
0337 }
0338 }
0339
0340 if (nVols != -1 && G4int(vvolu.size()) != nVols) {
0341 G4Exception("DicomIntersectVolume::GetLogicalVolumes:",
0342 "Wrong number of physical volumes found", FatalErrorInArgument,
0343 ("Number of physical volumes " + G4UIcommand::ConvertToString(G4int(vvolu.size()))
0344 + ", requesting " + G4UIcommand::ConvertToString(nVols))
0345 .c_str());
0346 }
0347
0348 return vvolu;
0349 }
0350
0351
0352 G4bool DicomIntersectVolume::IsPhantomVolume(G4VPhysicalVolume* pv)
0353 {
0354 EAxis axis;
0355 G4int nReplicas;
0356 G4double width, offset;
0357 G4bool consuming;
0358 pv->GetReplicationData(axis, nReplicas, width, offset, consuming);
0359 EVolume type = (consuming) ? kReplica : kParameterised;
0360 if (type == kParameterised && pv->GetRegularStructureId() == 1) {
0361 return TRUE;
0362 }
0363 else {
0364 return FALSE;
0365 }
0366 }
0367
0368
0369 std::vector<G4LogicalVolume*> DicomIntersectVolume::GetLogicalVolumes(const G4String& name,
0370 G4bool exists, G4int nVols)
0371 {
0372
0373 std::vector<G4LogicalVolume*> vvolu;
0374 G4int ial = G4int(name.rfind(":"));
0375 if (ial != -1) {
0376 G4Exception("DicomIntersectVolume::GetLogicalVolumes", "", FatalErrorInArgument,
0377 G4String("Name corresponds to a touchable or physcal volume" + name).c_str());
0378 }
0379
0380 G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
0381 for (auto citelv = lvs->cbegin(); citelv != lvs->cend(); ++citelv) {
0382 if (name == (*citelv)->GetName()) {
0383 vvolu.push_back(*citelv);
0384 }
0385 }
0386
0387
0388 if (vvolu.size() == 0) {
0389 if (exists) {
0390 G4Exception("DicomIntersectVolume::GetLogicalVolumes:", "", FatalErrorInArgument,
0391 ("no logical volume found with name " + name).c_str());
0392 }
0393 else {
0394 G4Exception("DicomIntersectVolume::GetLogicalVolumes:", "", JustWarning,
0395 ("no logical volume found with name " + name).c_str());
0396 }
0397 }
0398
0399 if (nVols != -1 && G4int(vvolu.size()) != nVols) {
0400 G4Exception("DicomIntersectVolume::GetLogicalVolumes:", "Wrong number of logical volumes found",
0401 FatalErrorInArgument,
0402 ("Number of logical volumes " + G4UIcommand::ConvertToString(G4int(vvolu.size()))
0403 + ", requesting " + G4UIcommand::ConvertToString(nVols))
0404 .c_str());
0405 }
0406
0407 return vvolu;
0408 }
0409
0410
0411 std::vector<G4String> DicomIntersectVolume::GetWordsInString(const G4String& stemp)
0412 {
0413 std::vector<G4String> wordlist;
0414
0415
0416
0417 G4int ii;
0418 const char* cstr = stemp.c_str();
0419 G4int siz = G4int(stemp.length());
0420 G4int istart = 0;
0421 G4int nQuotes = 0;
0422 G4bool lastIsBlank = false;
0423 G4bool lastIsQuote = false;
0424 for (ii = 0; ii < siz; ++ii) {
0425 if (cstr[ii] == '\"') {
0426 if (lastIsQuote) {
0427 G4Exception("GmGenUtils:GetWordsFromString", "", FatalException,
0428 ("There cannot be two quotes together " + stemp).c_str());
0429 }
0430 if (nQuotes % 2 == 1) {
0431
0432 wordlist.push_back(stemp.substr(istart, ii - istart));
0433
0434
0435
0436 }
0437 else {
0438 istart = ii + 1;
0439 }
0440 ++nQuotes;
0441 lastIsQuote = true;
0442 lastIsBlank = false;
0443 }
0444 else if (cstr[ii] == ' ') {
0445
0446
0447 if (nQuotes % 2 == 0) {
0448 if (!lastIsBlank && !lastIsQuote) {
0449 wordlist.push_back(stemp.substr(istart, ii - istart));
0450
0451
0452
0453 }
0454
0455 istart = ii + 1;
0456 lastIsQuote = false;
0457 lastIsBlank = true;
0458 }
0459 }
0460 else {
0461 if (ii == siz - 1) {
0462 wordlist.push_back(stemp.substr(istart, ii - istart + 1));
0463
0464
0465 }
0466 lastIsQuote = false;
0467 lastIsBlank = false;
0468 }
0469 }
0470 if (nQuotes % 2 == 1) {
0471 G4Exception("GmGenUtils:GetWordsFromString", "", FatalException,
0472 ("unbalanced quotes in line " + stemp).c_str());
0473 }
0474
0475 return wordlist;
0476 }
