88 if( params.size() < 8 ) {
89 G4Exception(
"DicomIntersectVolume::SetNewValue",
90 " There must be at least 8 parameter: SOLID_TYPE POS_X POS_Y POS_Z "
91 "ANG_X ANG_Y ANG_Z SOLID_PARAM_1 (SOLID_PARAM_2 ...)",
93 G4String(
"Number of parameters given = " +
94 G4UIcommand::ConvertToString( G4int(params.size()) )).c_str());
102 G4ThreeVector pos = G4ThreeVector( G4UIcommand::ConvertToDouble(params[0]),
103 G4UIcommand::ConvertToDouble(params[1]),
104 G4UIcommand::ConvertToDouble(params[2]) );
105 G4RotationMatrix* rotmat =
new G4RotationMatrix;
106 std::vector<G4double> angles;
107 rotmat->rotateX( G4UIcommand::ConvertToDouble(params[3]) );
108 rotmat->rotateY( G4UIcommand::ConvertToDouble(params[4]) );
109 rotmat->rotateY( G4UIcommand::ConvertToDouble(params[5]) );
116 if( params.size() !=1 )
117 G4Exception(
"DicomIntersectVolume::SetNewValue",
119 FatalErrorInArgument,
120 G4String(
"Command: "+ command->GetCommandPath() +
121 "/" + command->GetCommandName() +
123 " needs 1 argument: VOLUME_NAME").c_str());
138 G4RotationMatrix* rotph =
new G4RotationMatrix();
145 G4ThreeVector axisX( 1., 0., 0. );
146 G4ThreeVector axisY( 0., 1., 0. );
147 G4ThreeVector axisZ( 0., 0., 1. );
148 theTransform.ApplyAxisTransform(axisX);
149 theTransform.ApplyAxisTransform(axisY);
150 theTransform.ApplyAxisTransform(axisZ);
153 G4String thePhantomFileName =
"phantom.g4pdcm";
154 fout.open(thePhantomFileName);
155 std::vector<G4Material*> materials = thePhantomParam->GetMaterials();
156 fout << materials.size() << G4endl;
157 for(
unsigned int ii = 0; ii < materials.size(); ++ii )
159 fout << ii <<
" " << materials[ii]->GetName() << G4endl;
163 G4int nx = G4int(thePhantomParam->GetNoVoxelsX());
164 G4int ny = G4int(thePhantomParam->GetNoVoxelsY());
165 G4int nz = G4int(thePhantomParam->GetNoVoxelsZ());
168 G4double voxelHalfWidthX = thePhantomParam->GetVoxelHalfX();
169 G4double voxelHalfWidthY = thePhantomParam->GetVoxelHalfY();
170 G4double voxelHalfWidthZ = thePhantomParam->GetVoxelHalfZ();
173 fout << nx <<
" " << ny <<
" " << nz << G4endl;
174 fout << -voxelHalfWidthX*nx+thePhantomTransform.NetTranslation().x() <<
" "
175 << voxelHalfWidthX*nx+thePhantomTransform.NetTranslation().x() << G4endl;
176 fout << -voxelHalfWidthY*ny+thePhantomTransform.NetTranslation().y() <<
" "
177 << voxelHalfWidthY*ny+thePhantomTransform.NetTranslation().y() << G4endl;
178 fout << -voxelHalfWidthZ*nz+thePhantomTransform.NetTranslation().z() <<
" "
179 << voxelHalfWidthZ*nz+thePhantomTransform.NetTranslation().z() << G4endl;
181 for( G4int iz = 0; iz < nz; ++iz ){
182 for( G4int iy = 0; iy < ny; ++iy) {
184 G4bool bPrevVoxelInside =
true;
185 G4bool b1VoxelFoundInside =
false;
186 G4int firstVoxel = -1;
187 G4int lastVoxel = -1;
188 for(G4int ix = 0; ix < nx; ++ix ){
189 G4ThreeVector voxelCentre( (-nx+ix*2+1)*voxelHalfWidthX,
190 (-ny+iy*2+1)*voxelHalfWidthY, (-nz+iz*2+1)*voxelHalfWidthZ);
191 theTransform.ApplyPointTransform(voxelCentre);
192 G4bool bVoxelIsInside =
true;
193 for( G4int ivx = -1; ivx <= 1; ivx+=2 ) {
194 for( G4int ivy = -1; ivy <= 1; ivy+=2 ){
195 for( G4int ivz = -1; ivz <= 1; ivz+=2 ) {
196 G4ThreeVector voxelPoint = voxelCentre
197 + ivx*voxelHalfWidthX*axisX +
198 ivy*voxelHalfWidthY*axisY + ivz*voxelHalfWidthZ*axisZ;
199 if(
fSolid->Inside( voxelPoint ) == kOutside ) {
200 bVoxelIsInside =
false;
205 if( !bVoxelIsInside )
break;
207 if( !bVoxelIsInside )
break;
210 G4int copyNo = ix + nx*iy + nxy*iz;
211 if( bVoxelIsInside ) {
212 if( !bPrevVoxelInside ) {
213 G4Exception(
"DicomIntersectVolume::SetNewValue",
216 "Volume cannot intersect phantom in discontiguous voxels, "
217 "please use other voxel");
219 if( !b1VoxelFoundInside ) {
221 b1VoxelFoundInside =
true;
230 fout << firstVoxel <<
" " << lastVoxel << G4endl;
235 for( G4int iz = 0; iz < nz; ++iz ){
236 for( G4int iy = 0; iy < ny; ++iy) {
237 G4bool b1xFound =
false;
238 for(G4int ix = 0; ix < nx; ++ix ){
239 size_t copyNo = ix + ny*iy + nxy*iz;
241 fout << thePhantomParam->GetMaterialIndex(copyNo)<<
" ";
245 if(b1xFound )
fout << G4endl;
250 for( G4int iz = 0; iz < nz; ++iz ){
251 for( G4int iy = 0; iy < ny; ++iy) {
252 G4bool b1xFound =
false;
253 for(G4int ix = 0; ix < nx; ++ix ){
254 size_t copyNo = ix + ny*iy + nxy*iz;
256 fout <<thePhantomParam->GetMaterial(copyNo)->GetDensity()/g*cm3<<
" ";
260 if(b1xFound )
fout << G4endl;
319 const G4String& name, G4bool exists, G4int nVols )
321 std::vector<G4VPhysicalVolume*> vvolu;
322 std::string::size_type ial = name.rfind(
":");
325 if( ial != std::string::npos )
327 std::string::size_type ial2 = name.rfind(
":",ial-1);
328 if( ial2 != std::string::npos )
330 G4Exception(
"DicomIntersectVolume::GetPhysicalVolumes",
332 FatalErrorInArgument,
333 G4String(
"Name corresponds to a touchable " + name).c_str());
337 volname = name.substr( 0, ial );
338 volcopy = G4UIcommand::ConvertToInt( name.substr(ial+1, name.length()).c_str() );
348 for(
auto citepv = pvs->cbegin(); citepv != pvs->cend(); ++citepv )
350 if( volname == (*citepv)->GetName()
351 && ( (*citepv)->GetCopyNo() == volcopy || -1 == volcopy ) )
353 vvolu.push_back( *citepv );
358 if( vvolu.size() == 0 ) {
360 G4Exception(
" DicomIntersectVolume::GetPhysicalVolumes",
362 FatalErrorInArgument,
363 G4String(
"No physical volume found with name " + name).c_str());
365 G4cerr <<
"!!WARNING: DicomIntersectVolume::GetPhysicalVolumes: "<<
366 " no physical volume found with name " << name << G4endl;
370 if( nVols != -1 && G4int(vvolu.size()) != nVols ) {
371 G4Exception(
"DicomIntersectVolume::GetLogicalVolumes:",
372 "Wrong number of physical volumes found",
373 FatalErrorInArgument,
374 (
"Number of physical volumes " +
375 G4UIcommand::ConvertToString(G4int(vvolu.size())) +
376 ", requesting " + G4UIcommand::ConvertToString(nVols)).c_str());
401const G4String& name, G4bool exists, G4int nVols )
404 std::vector<G4LogicalVolume*> vvolu;
405 G4int ial = G4int(name.rfind(
":"));
407 G4Exception(
"DicomIntersectVolume::GetLogicalVolumes",
409 FatalErrorInArgument,
410 G4String(
"Name corresponds to a touchable or physcal volume"
415 for(
auto citelv = lvs->cbegin(); citelv != lvs->cend(); ++citelv )
417 if( name == (*citelv)->GetName() ) {
418 vvolu.push_back( *citelv );
423 if( vvolu.size() == 0 ) {
425 G4Exception(
"DicomIntersectVolume::GetLogicalVolumes:",
"",
426 FatalErrorInArgument,(
"no logical volume found with name "
429 G4Exception(
"DicomIntersectVolume::GetLogicalVolumes:",
"",
430 JustWarning,(
"no logical volume found with name " + name).c_str());
434 if( nVols != -1 && G4int(vvolu.size()) != nVols ) {
435 G4Exception(
"DicomIntersectVolume::GetLogicalVolumes:",
436 "Wrong number of logical volumes found",
437 FatalErrorInArgument,
438 (
"Number of logical volumes " +
439 G4UIcommand::ConvertToString(G4int(vvolu.size())) +
440 ", requesting " + G4UIcommand::ConvertToString(nVols)).c_str());