Context Navigation

DicomDetectorConstruction.cc @ 1337

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago
tag geant4.9.4 beta 1 + modifs locales
File size: 17.7 KB

Line
1	//
2	// ********************************************************************
3	// * License and Disclaimer *
4	// * *
5	// * The Geant4 software is copyright of the Copyright Holders of *
6	// * the Geant4 Collaboration. It is provided under the terms and *
7	// * conditions of the Geant4 Software License, included in the file *
8	// * LICENSE and available at http://cern.ch/geant4/license . These *
9	// * include a list of copyright holders. *
10	// * *
11	// * Neither the authors of this software system, nor their employing *
12	// * institutes,nor the agencies providing financial support for this *
13	// * work make any representation or warranty, express or implied, *
14	// * regarding this software system or assume any liability for its *
15	// * use. Please see the license in the file LICENSE and URL above *
16	// * for the full disclaimer and the limitation of liability. *
17	// * *
18	// * This code implementation is the result of the scientific and *
19	// * technical work of the GEANT4 collaboration. *
20	// * By using, copying, modifying or distributing the software (or *
21	// * any work based on the software) you agree to acknowledge its *
22	// * use in resulting scientific publications, and indicate your *
23	// * acceptance of all terms of the Geant4 Software license. *
24	// ********************************************************************
25	//
26
27	#include "globals.hh"
28
29	#include "G4Box.hh"
30	#include "G4LogicalVolume.hh"
31	#include "G4VPhysicalVolume.hh"
32	#include "G4PVPlacement.hh"
33	#include "G4Material.hh"
34	#include "G4Element.hh"
35
36	#include "DicomDetectorConstruction.hh"
37	#include "DicomPatientZSliceHeader.hh"
38
39	//-------------------------------------------------------------
40	DicomDetectorConstruction::DicomDetectorConstruction()
41	{
42	}
43
44	//-------------------------------------------------------------
45	DicomDetectorConstruction::~DicomDetectorConstruction()
46	{
47	}
48
49	//-------------------------------------------------------------
50	G4VPhysicalVolume* DicomDetectorConstruction::Construct()
51	{
52	InitialisationOfMaterials();
53
54	//----- Build world
55	G4double worldXDimension = 1.*m;
56	G4double worldYDimension = 1.*m;
57	G4double worldZDimension = 1.*m;
58
59	world_solid = new G4Box( "WorldSolid",
60	worldXDimension,
61	worldYDimension,
62	worldZDimension );
63
64	world_logic = new G4LogicalVolume( world_solid,
65	air,
66	"WorldLogical",
67	0, 0, 0 );
68
69	world_phys = new G4PVPlacement( 0,
70	G4ThreeVector(0,0,0),
71	"World",
72	world_logic,
73	0,
74	false,
75	0 );
76
77	ReadPatientData();
78
79	ConstructPatientContainer();
80	ConstructPatient();
81
82	return world_phys;
83	}
84
85
86	//-------------------------------------------------------------
87	void DicomDetectorConstruction::InitialisationOfMaterials()
88	{
89	// Creating elements :
90	G4double z, a, density;
91	G4String name, symbol;
92
93	G4Element* elC = new G4Element( name = "Carbon",
94	symbol = "C",
95	z = 6.0, a = 12.011 * g/mole );
96	G4Element* elH = new G4Element( name = "Hydrogen",
97	symbol = "H",
98	z = 1.0, a = 1.008 * g/mole );
99	G4Element* elN = new G4Element( name = "Nitrogen",
100	symbol = "N",
101	z = 7.0, a = 14.007 * g/mole );
102	G4Element* elO = new G4Element( name = "Oxygen",
103	symbol = "O",
104	z = 8.0, a = 16.00 * g/mole );
105	G4Element* elNa = new G4Element( name = "Sodium",
106	symbol = "Na",
107	z= 11.0, a = 22.98977* g/mole );
108	G4Element* elS = new G4Element( name = "Sulfur",
109	symbol = "S",
110	z = 16.0,a = 32.065* g/mole );
111	G4Element* elCl = new G4Element( name = "Chlorine",
112	symbol = "P",
113	z = 17.0, a = 35.453* g/mole );
114	G4Element* elK = new G4Element( name = "Potassium",
115	symbol = "P",
116	z = 19.0, a = 30.0983* g/mole );
117	G4Element* elP = new G4Element( name = "Phosphorus",
118	symbol = "P",
119	z = 30.0, a = 30.973976* g/mole );
120	G4Element* elFe = new G4Element( name = "Iron",
121	symbol = "Fe",
122	z = 26, a = 56.845* g/mole );
123	G4Element* elMg = new G4Element( name = "Magnesium",
124	symbol = "Mg",
125	z = 12.0, a = 24.3050* g/mole );
126	G4Element* elCa = new G4Element( name="Calcium",
127	symbol = "Ca",
128	z = 20.0, a = 40.078* g/mole );
129
130	// Creating Materials :
131	G4int numberofElements;
132
133	// Air
134	air = new G4Material( "Air",
135	1.290*mg/cm3,
136	numberofElements = 2 );
137	air->AddElement(elN, 0.7);
138	air->AddElement(elO, 0.3);
139
140	// Lung Inhale
141	G4Material* lunginhale = new G4Material( "LungInhale",
142	density = 0.217*g/cm3,
143	numberofElements = 9);
144	lunginhale->AddElement(elH,0.103);
145	lunginhale->AddElement(elC,0.105);
146	lunginhale->AddElement(elN,0.031);
147	lunginhale->AddElement(elO,0.749);
148	lunginhale->AddElement(elNa,0.002);
149	lunginhale->AddElement(elP,0.002);
150	lunginhale->AddElement(elS,0.003);
151	lunginhale->AddElement(elCl,0.002);
152	lunginhale->AddElement(elK,0.003);
153
154	// Lung exhale
155	G4Material* lungexhale = new G4Material( "LungExhale",
156	density = 0.508*g/cm3,
157	numberofElements = 9 );
158	lungexhale->AddElement(elH,0.103);
159	lungexhale->AddElement(elC,0.105);
160	lungexhale->AddElement(elN,0.031);
161	lungexhale->AddElement(elO,0.749);
162	lungexhale->AddElement(elNa,0.002);
163	lungexhale->AddElement(elP,0.002);
164	lungexhale->AddElement(elS,0.003);
165	lungexhale->AddElement(elCl,0.002);
166	lungexhale->AddElement(elK,0.003);
167
168	// Adipose tissue
169	G4Material* adiposeTissue = new G4Material( "AdiposeTissue",
170	density = 0.967*g/cm3,
171	numberofElements = 7);
172	adiposeTissue->AddElement(elH,0.114);
173	adiposeTissue->AddElement(elC,0.598);
174	adiposeTissue->AddElement(elN,0.007);
175	adiposeTissue->AddElement(elO,0.278);
176	adiposeTissue->AddElement(elNa,0.001);
177	adiposeTissue->AddElement(elS,0.001);
178	adiposeTissue->AddElement(elCl,0.001);
179
180	// Breast
181	G4Material* breast = new G4Material( "Breast",
182	density = 0.990*g/cm3,
183	numberofElements = 8 );
184	breast->AddElement(elH,0.109);
185	breast->AddElement(elC,0.506);
186	breast->AddElement(elN,0.023);
187	breast->AddElement(elO,0.358);
188	breast->AddElement(elNa,0.001);
189	breast->AddElement(elP,0.001);
190	breast->AddElement(elS,0.001);
191	breast->AddElement(elCl,0.001);
192
193	// Water
194	G4Material* water = new G4Material( "Water",
195	density = 1.0*g/cm3,
196	numberofElements = 2 );
197	water->AddElement(elH,0.112);
198	water->AddElement(elO,0.888);
199
200	// Muscle
201	G4Material* muscle = new G4Material( "Muscle",
202	density = 1.061*g/cm3,
203	numberofElements = 9 );
204	muscle->AddElement(elH,0.102);
205	muscle->AddElement(elC,0.143);
206	muscle->AddElement(elN,0.034);
207	muscle->AddElement(elO,0.710);
208	muscle->AddElement(elNa,0.001);
209	muscle->AddElement(elP,0.002);
210	muscle->AddElement(elS,0.003);
211	muscle->AddElement(elCl,0.001);
212	muscle->AddElement(elK,0.004);
213
214	// Liver
215	G4Material* liver = new G4Material( "Liver",
216	density = 1.071*g/cm3,
217	numberofElements = 9);
218	liver->AddElement(elH,0.102);
219	liver->AddElement(elC,0.139);
220	liver->AddElement(elN,0.030);
221	liver->AddElement(elO,0.716);
222	liver->AddElement(elNa,0.002);
223	liver->AddElement(elP,0.003);
224	liver->AddElement(elS,0.003);
225	liver->AddElement(elCl,0.002);
226	liver->AddElement(elK,0.003);
227
228	// Trabecular Bone
229	G4Material* trabecularBone = new G4Material( "TrabecularBone",
230	density = 1.159*g/cm3,
231	numberofElements = 12 );
232	trabecularBone->AddElement(elH,0.085);
233	trabecularBone->AddElement(elC,0.404);
234	trabecularBone->AddElement(elN,0.058);
235	trabecularBone->AddElement(elO,0.367);
236	trabecularBone->AddElement(elNa,0.001);
237	trabecularBone->AddElement(elMg,0.001);
238	trabecularBone->AddElement(elP,0.034);
239	trabecularBone->AddElement(elS,0.002);
240	trabecularBone->AddElement(elCl,0.002);
241	trabecularBone->AddElement(elK,0.001);
242	trabecularBone->AddElement(elCa,0.044);
243	trabecularBone->AddElement(elFe,0.001);
244
245	// Dense Bone
246	G4Material* denseBone = new G4Material( "DenseBone",
247	density = 1.575*g/cm3,
248	numberofElements = 11 );
249	denseBone->AddElement(elH,0.056);
250	denseBone->AddElement(elC,0.235);
251	denseBone->AddElement(elN,0.050);
252	denseBone->AddElement(elO,0.434);
253	denseBone->AddElement(elNa,0.001);
254	denseBone->AddElement(elMg,0.001);
255	denseBone->AddElement(elP,0.072);
256	denseBone->AddElement(elS,0.003);
257	denseBone->AddElement(elCl,0.001);
258	denseBone->AddElement(elK,0.001);
259	denseBone->AddElement(elCa,0.146);
260
261	//----- Put the materials in a vector
262	fOriginalMaterials.push_back(air); // rho = 0.00129
263	fOriginalMaterials.push_back(lunginhale); // rho = 0.217
264	fOriginalMaterials.push_back(lungexhale); // rho = 0.508
265	fOriginalMaterials.push_back(adiposeTissue); // rho = 0.967
266	fOriginalMaterials.push_back(breast ); // rho = 0.990
267	fOriginalMaterials.push_back(water); // rho = 1.018
268	fOriginalMaterials.push_back(muscle); // rho = 1.061
269	fOriginalMaterials.push_back(liver); // rho = 1.071
270	fOriginalMaterials.push_back(trabecularBone); // rho = 1.159
271	fOriginalMaterials.push_back(denseBone); // rho = 1.575
272
273	}
274
275
276	//-------------------------------------------------------------
277	void DicomDetectorConstruction::ReadPatientData()
278	{
279	std::ifstream finDF("Data.dat");
280	G4String fname;
281	if(finDF.good() != 1 ) {
282	G4Exception(" DicomDetectorConstruction::ReadPatientData. Problem reading data file: Data.dat");
283	}
284
285	G4int compression;
286	finDF >> compression; // not used here
287
288	finDF >> fNoFiles;
289	for(G4int i = 0; i < fNoFiles; i++ ) {
290	finDF >> fname;
291	//--- Read one data file
292	fname += ".g4dcm";
293	ReadPatientDataFile(fname);
294	}
295
296	//----- Merge data headers
297	MergeZSliceHeaders();
298
299	finDF.close();
300
301	}
302
303	//-------------------------------------------------------------
304	void DicomDetectorConstruction::ReadPatientDataFile(const G4String& fname)
305	{
306	#ifdef G4VERBOSE
307	G4cout << " DicomDetectorConstruction::ReadPatientDataFile opening file " << fname << G4endl;
308	#endif
309	std::ifstream fin(fname.c_str(), std::ios_base::in);
310	if( !fin.is_open() ) {
311	G4Exception("DicomDetectorConstruction::ReadPatientDataFil. File not found " + fname );
312	}
313	//----- Define density differences (maximum density difference to create a new material)
314	G4double densityDiff = 0.1;
315	std::map<G4int,G4double> fDensityDiffs; // to be able to use a different densityDiff for each material
316	for( unsigned int ii = 0; ii < fOriginalMaterials.size(); ii++ ){
317	fDensityDiffs[ii] = densityDiff; //currently all materials with same difference
318	}
319
320	//----- Read data header
321	DicomPatientZSliceHeader* sliceHeader = new DicomPatientZSliceHeader( fin );
322	fZSliceHeaders.push_back( sliceHeader );
323
324	//----- Read material indices
325	G4int nVoxels = sliceHeader->GetNoVoxels();
326
327	//--- If first slice, initiliaze fMateIDs
328	if( fZSliceHeaders.size() == 1 ) {
329	//fMateIDs = new unsigned int[fNoFiles*nVoxels];
330	fMateIDs = new size_t[fNoFiles*nVoxels];
331
332	}
333
334	unsigned int mateID;
335	G4int voxelCopyNo = (fZSliceHeaders.size()-1)*nVoxels; // number of voxels from previously read slices
336	for( G4int ii = 0; ii < nVoxels; ii++, voxelCopyNo++ ){
337	fin >> mateID;
338	fMateIDs[voxelCopyNo] = mateID;
339	}
340
341	//----- Read material densities and build new materials if two voxels have same material but its density is in a different density interval (size of density intervals defined by densityDiff)
342	G4double density;
343	voxelCopyNo = (fZSliceHeaders.size()-1)*nVoxels; // number of voxels from previously read slices
344	for( G4int ii = 0; ii < nVoxels; ii++, voxelCopyNo++ ){
345	fin >> density;
346
347	//-- Get material from list of original materials
348	int mateID = fMateIDs[voxelCopyNo];
349	G4Material* mateOrig = fOriginalMaterials[mateID];
350
351	//-- Get density bin: middle point of the bin in which the current density is included
352	float densityBin = fDensityDiffs[mateID] * (G4int(density/fDensityDiffs[mateID])+0.5);
353	//-- Build the new material name
354	G4String newMateName = mateOrig->GetName()+"__"+ftoa(densityBin);
355
356	//-- Look if a material with this name is already created (because a previous voxel was already in this density bin)
357	unsigned int im;
358	for( im = 0; im < fMaterials.size(); im++ ){
359	if( fMaterials[im]->GetName() == newMateName ) {
360	break;
361	}
362	}
363	//-- If material is already created use index of this material
364	if( im != fMaterials.size() ) {
365	fMateIDs[voxelCopyNo] = im;
366	//-- else, create the material
367	} else {
368	fMaterials.push_back( BuildMaterialWithChangingDensity( mateOrig, densityBin, newMateName ) );
369	fMateIDs[voxelCopyNo] = fMaterials.size()-1;
370	}
371	}
372
373	}
374
375
376	//-------------------------------------------------------------
377	void DicomDetectorConstruction::MergeZSliceHeaders()
378	{
379	//----- Images must have the same dimension ...
380	fZSliceHeaderMerged = new DicomPatientZSliceHeader( *fZSliceHeaders[0] );
381	for( unsigned int ii = 1; ii < fZSliceHeaders.size(); ii++ ) {
382	fZSliceHeaderMerged += fZSliceHeaders[ii];
383	};
384	}
385
386	//-------------------------------------------------------------
387	G4Material* DicomDetectorConstruction::BuildMaterialWithChangingDensity( const G4Material* origMate, float density, G4String newMateName )
388	{
389	//----- Copy original material, but with new density
390	G4int nelem = origMate->GetNumberOfElements();
391	G4Material* mate = new G4Material( newMateName, density*g/cm3, nelem, kStateUndefined, STP_Temperature );
392
393	for( G4int ii = 0; ii < nelem; ii++ ){
394	G4double frac = origMate->GetFractionVector()[ii];
395	G4Element* elem = const_cast<G4Element*>(origMate->GetElement(ii));
396	mate->AddElement( elem, frac );
397	}
398
399	return mate;
400	}
401
402	//-----------------------------------------------------------------------
403	G4String DicomDetectorConstruction::ftoa(float flo)
404	{
405	char ctmp[100];
406	gcvt( flo, 10, ctmp );
407	return G4String(ctmp);
408	}
409
410
411	//-------------------------------------------------------------
412	void DicomDetectorConstruction::ConstructPatientContainer()
413	{
414	//---- Extract number of voxels and voxel dimensions
415	nVoxelX = fZSliceHeaderMerged->GetNoVoxelX();
416	nVoxelY = fZSliceHeaderMerged->GetNoVoxelY();
417	nVoxelZ = fZSliceHeaderMerged->GetNoVoxelZ();
418
419	voxelHalfDimX = fZSliceHeaderMerged->GetVoxelHalfX();
420	voxelHalfDimY = fZSliceHeaderMerged->GetVoxelHalfY();
421	voxelHalfDimZ = fZSliceHeaderMerged->GetVoxelHalfZ();
422	#ifdef G4VERBOSE
423	G4cout << " nVoxelX " << nVoxelX << " voxelHalfDimX " << voxelHalfDimX <<G4endl;
424	G4cout << " nVoxelY " << nVoxelY << " voxelHalfDimY " << voxelHalfDimY <<G4endl;
425	G4cout << " nVoxelZ " << nVoxelZ << " voxelHalfDimZ " << voxelHalfDimZ <<G4endl;
426	G4cout << " totalPixels " << nVoxelXnVoxelYnVoxelZ << G4endl;
427	#endif
428
429	//----- Define the volume that contains all the voxels
430	container_solid = new G4Box("PhantomContainer",nVoxelXvoxelHalfDimX,nVoxelYvoxelHalfDimY,nVoxelZ*voxelHalfDimZ);
431	container_logic =
432	new G4LogicalVolume( container_solid,
433	fMaterials[0], //the material is not important, it will be fully filled by the voxels
434	"PhantomContainer",
435	0, 0, 0 );
436	//--- Place it on the world
437	G4double offsetX = (fZSliceHeaderMerged->GetMaxX() + fZSliceHeaderMerged->GetMinX() ) /2.;
438	G4double offsetY = (fZSliceHeaderMerged->GetMaxY() + fZSliceHeaderMerged->GetMinY() ) /2.;
439	G4double offsetZ = (fZSliceHeaderMerged->GetMaxZ() + fZSliceHeaderMerged->GetMinZ() ) /2.;
440	G4ThreeVector posCentreVoxels(offsetX,offsetY,offsetZ);
441	#ifdef G4VERBOSE
442	G4cout << " placing voxel container volume at " << posCentreVoxels << G4endl;
443	#endif
444	container_phys =
445	new G4PVPlacement(0, // rotation
446	posCentreVoxels,
447	container_logic, // The logic volume
448	"PhantomContainer", // Name
449	world_logic, // Mother
450	false, // No op. bool.
451	1); // Copy number
452
453	}
454
455
456	#include "G4SDManager.hh"
457	#include "G4MultiFunctionalDetector.hh"
458	#include "G4PSDoseDeposit_RegNav.hh"
459
460	//-------------------------------------------------------------
461	void DicomDetectorConstruction::SetScorer(G4LogicalVolume* voxel_logic)
462	{
463
464	G4SDManager* SDman = G4SDManager::GetSDMpointer();
465	//
466	// Sensitive Detector Name
467	G4String concreteSDname = "PatientSD";
468
469	//------------------------
470	// MultiFunctionalDetector
471	//------------------------
472	//
473	// Define MultiFunctionalDetector with name.
474	G4MultiFunctionalDetector* MFDet = new G4MultiFunctionalDetector(concreteSDname);
475	SDman->AddNewDetector( MFDet ); // Register SD to SDManager
476
477	voxel_logic->SetSensitiveDetector(MFDet);
478
479	G4PSDoseDeposit_RegNav* scorer = new G4PSDoseDeposit_RegNav("DoseDeposit");
480	MFDet->RegisterPrimitive(scorer);
481
482	}
483

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: trunk/examples/extended/medical/DICOM/src/DicomDetectorConstruction.cc @ 1337

Download in other formats: