Rev 1 | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
1 | f9daq | 1 | // |
2 | // ../bin/FBP2D FBP2D.par |
||
3 | // |
||
4 | // ******************************************************************** |
||
5 | // * License and Disclaimer * |
||
6 | // * * |
||
7 | // * The GAMOS software is copyright of the Copyright Holders of * |
||
8 | // * the GAMOS Collaboration. It is provided under the terms and * |
||
9 | // * conditions of the GAMOS Software License, included in the file * |
||
10 | // * LICENSE and available at http://fismed.ciemat.es/GAMOS/license .* |
||
11 | // * These include a list of copyright holders. * |
||
12 | // * * |
||
13 | // * Neither the authors of this software system, nor their employing * |
||
14 | // * institutes,nor the agencies providing financial support for this * |
||
15 | // * work make any representation or warranty, express or implied, * |
||
16 | // * regarding this software system or assume any liability for its * |
||
17 | // * use. Please see the license in the file LICENSE and URL above * |
||
18 | // * for the full disclaimer and the limitation of liability. * |
||
19 | // * * |
||
20 | // * This code implementation is the result of the scientific and * |
||
21 | // * technical work of the GAMOS collaboration. * |
||
22 | // * By using, copying, modifying or distributing the software (or * |
||
23 | // * any work based on the software) you agree to acknowledge its * |
||
24 | // * use in resulting scientific publications, and indicate your * |
||
25 | // * acceptance of all terms of the GAMOS Software license. * |
||
26 | // ******************************************************************** |
||
27 | // |
||
28 | #include "PETProjDataMgr.h" |
||
29 | #include "TH2F.h" |
||
30 | #include "TH3F.h" |
||
31 | #include "TMath.h" |
||
32 | #include "TRandom3.h" |
||
33 | |||
34 | #include <iomanip> |
||
35 | #include <iostream> |
||
36 | #include <cstdlib> |
||
37 | |||
38 | using namespace std; |
||
39 | |||
40 | //---------------------------------------------------------------------- |
||
41 | PETProjDataMgr* PETProjDataMgr::m_Instance = 0; |
||
42 | |||
43 | //---------------------------------------------------------------------- |
||
44 | PETProjDataMgr* PETProjDataMgr::GetInstance() |
||
45 | { |
||
46 | if( !m_Instance ) { |
||
47 | m_Instance = new PETProjDataMgr; |
||
48 | } |
||
49 | |||
50 | return m_Instance; |
||
51 | |||
52 | } |
||
53 | |||
54 | //----------------------------------------------------------------------- |
||
55 | PETProjDataMgr::PETProjDataMgr() |
||
56 | { |
||
57 | |||
58 | /* |
||
59 | // ARGUMENTS: |
||
60 | " cout << " -------------------------- \n" |
||
61 | " Arguments convention: \n" |
||
62 | " -a Axial FOV (mm), <theDist_axial=100.0> \n" |
||
63 | " -d Diameter Transaxial FOV (mm), <m_RingDiameter=300.0> \n" |
||
64 | " -i Type of the input file (by default: 0 = Arce_binary), <typeINfile=0> \n" |
||
65 | " \n" |
||
66 | " -m Maximum ring difference (by default: -1 = m_NOfPlanes), <m_MaxRingDifferenceiff> \n" |
||
67 | " -n Name of output file, <m_Filename> \n" |
||
68 | " -p Axial number of planes, <m_NOfPlanes> \n" |
||
69 | " -r Number of bins, \"distancias\", <m_NOfBins> \n" |
||
70 | // -s Span TO DO:span !!!!! |
||
71 | " -t Number of angular views, \"direcciones\", <m_NOfAngles> \n" |
||
72 | " -v Verbosity (by default: 0=silent, 3=debug), <verbos> \n" |
||
73 | " -x Maximum number of coincidences to be stored (by default: -1 = no limit), <Max_Coinci> \n" |
||
74 | " -o Output type (by default: 0 = mcc Interfile, 1 = STIR Interfile), <OutType> \n" |
||
75 | |||
76 | " \n" |
||
77 | " PET Reconstruction. CIEMAT 2009-11 \n" |
||
78 | " mario.canadas@ciemat.es \n" |
||
79 | " -------------------------- \n"; |
||
80 | */ |
||
81 | m_rnd= new TRandom3(); |
||
82 | m_AxialDistance = (9-1)*2.25; // Axial pixel dimension*NOfPlanes |
||
83 | m_RingDiameter = 120.0; // notranji premer peta |
||
84 | m_NOfPlanes = 9; // stevilo ravnin |
||
85 | m_NOfBins = 128; // stevilo binov v razdalji |
||
86 | m_nch = 128; // stevilo padov okoli in okoli |
||
87 | m_NOfAngles = TMath::Abs(m_nch)/2; // stevilo kotov = stevilo padov okoli in okoli /2 |
||
88 | |||
89 | m_MaxRingDifference = -1; |
||
90 | //m_MaxRingDifference = 3; // najvecja razdalja med padi |
||
91 | // toDo: theSpan = int(GmParameterMgr::GetInstance()->GetNumericValue("PET:ProjData:Span",1)); |
||
92 | |||
93 | |||
94 | m_OutFormat = 1; // 1.. projections 0 .. image |
||
95 | m_Filename = "sino3D"; |
||
96 | m_Debug=1; |
||
97 | |||
98 | if (m_MaxRingDifference==-1) m_MaxRingDifference=m_NOfPlanes-1; |
||
99 | |||
100 | |||
101 | m_TotalAxialPlanes=m_NOfPlanes*m_NOfPlanes; |
||
102 | if (m_OutFormat==1) m_TotalAxialPlanes= (2*m_NOfPlanes-1 - m_MaxRingDifference)*m_MaxRingDifference + m_NOfPlanes; // total number of Axial planes (segments*planes) in STIR format |
||
103 | |||
104 | /*--- Initialize sino3D ---*/ |
||
105 | m_projections = new SINO_TYPE**[m_NOfBins]; |
||
106 | for(int i=0;i<m_NOfBins;i++){ |
||
107 | m_projections[i] = new SINO_TYPE*[m_NOfAngles]; |
||
108 | for(int j=0;j<m_NOfAngles;j++){ |
||
109 | m_projections[i][j] = new SINO_TYPE[m_TotalAxialPlanes]; /// ! If m_OutFormat==1 (STIR output):Matrix size depends on the MAX_Ring_Difference |
||
110 | for(int k=0;k<m_TotalAxialPlanes;k++){ |
||
111 | m_projections[i][j][k]=0; |
||
112 | } |
||
113 | } |
||
114 | } |
||
115 | |||
116 | m_TotalProjectionCoincidences=0; |
||
117 | m_TotalCoincidences=0; |
||
118 | //OutputType = "pet"; |
||
119 | } |
||
120 | |||
121 | //----------------------------------------------------------------------- |
||
122 | void PETProjDataMgr::SetProjection( int axialplane, TH2F * h) |
||
123 | { |
||
124 | for(int i=0;i<m_NOfBins;i++){ |
||
125 | for(int j=0;j<m_NOfAngles;j++){ |
||
126 | m_projections[i][j][axialplane]=h->GetBinContent(i+1,j+1); |
||
127 | } |
||
128 | } |
||
129 | |||
130 | } |
||
131 | //----------------------------------------------------------------------- |
||
132 | |||
133 | //----------------------------------------------------------- |
||
134 | // from Gate |
||
135 | //------------------------------------------------------------ |
||
136 | double ComputeSinogramS(double X1, double Y1, double X2, double Y2) |
||
137 | { |
||
138 | double s; |
||
139 | |||
140 | double denom = (Y1-Y2) * (Y1-Y2) + (X2-X1) * (X2-X1); |
||
141 | |||
142 | if (denom!=0.) { |
||
143 | denom = sqrt(denom); |
||
144 | s = ( X1 * (Y2-Y1) + Y1 * (X1-X2) ) / denom; |
||
145 | } else { |
||
146 | s = 0.; |
||
147 | } |
||
148 | |||
149 | double theta; |
||
150 | if ((X1-X2)!=0.) { |
||
151 | theta=atan((X1-X2) /(Y1-Y2)); |
||
152 | } else { |
||
153 | theta=3.1416/2.; |
||
154 | } |
||
155 | if ((theta > 0.) && ((X1-X2) > 0.)) s = -s; |
||
156 | if ((theta < 0.) && ((X1-X2) < 0.)) s = -s; |
||
157 | if ( theta < 0.) { |
||
158 | theta = theta+3.1416; |
||
159 | s = -s; |
||
160 | } |
||
161 | return s; |
||
162 | } |
||
163 | |||
164 | |||
165 | void PETProjDataMgr::AddEvent( const TVector3& pos1 , const TVector3& pos2) |
||
166 | { |
||
167 | int z1_i, z2_i; |
||
168 | //for discretization on the crystal: int x1_i, x2_i, y1_i, y2_i; |
||
169 | |||
170 | double z1_abs=pos1.z()+m_AxialDistance/2; |
||
171 | double z2_abs=pos2.z()+m_AxialDistance/2; |
||
172 | double a, b, phi, dis; |
||
173 | int phi_i, dis_i; |
||
174 | int ring_diff; |
||
175 | |||
176 | double _PI=2*asin(1); |
||
177 | |||
178 | m_TotalCoincidences++; |
||
179 | |||
180 | z1_i=(int)(m_NOfPlanes* z1_abs/m_AxialDistance); //round --> m_NOfPlanes+1 ... |
||
181 | z2_i=(int)(m_NOfPlanes* z2_abs/m_AxialDistance); |
||
182 | |||
183 | // control; if z_i out of range: return |
||
184 | |||
185 | if ( (pos1.x()==pos2.x()) && (pos1.y()==pos2.y()) ) { |
||
186 | #ifndef GAMOS_NO_VERBOSE |
||
187 | if( m_Debug ) { |
||
188 | cout << "PETProjDataMgr::AddEvent:WARNING! Event_1 == Event_2 ; x= " << pos2.x() << ", y= " << pos2.y() << ", z= " << pos2.z() << endl; |
||
189 | } |
||
190 | #endif |
||
191 | return; |
||
192 | } |
||
193 | |||
194 | if ( (z1_i<0) || (z2_i<0) || (z1_i>= m_NOfPlanes) || (z2_i>= m_NOfPlanes) ) { |
||
195 | #ifndef GAMOS_NO_VERBOSE |
||
196 | if( m_Debug ) { |
||
197 | cout << "PETProjDataMgr::AddEvent:WARNING! Event out of bounds (Axial): x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl; |
||
198 | } |
||
199 | #endif |
||
200 | return; |
||
201 | } |
||
202 | |||
203 | ring_diff = (int)fabs(z1_i-z2_i); |
||
204 | |||
205 | // max ring difference; control: |
||
206 | if (ring_diff > m_MaxRingDifference) { |
||
207 | #ifndef GAMOS_NO_VERBOSE |
||
208 | if( m_Debug ) { |
||
209 | cout <<"PETProjDataMgr::AddEvent:WARNING! Event out of bounds (Max. Ring Diff.): " << ring_diff << ">" << m_MaxRingDifference << " x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl; |
||
210 | } |
||
211 | #endif |
||
212 | return; |
||
213 | } |
||
214 | |||
215 | a=(double)(pos2.y()- pos1.y()); |
||
216 | b=(double)(pos2.x()- pos1.x()); |
||
217 | |||
218 | if (a==0.0){ |
||
219 | phi=_PI*0.5; |
||
220 | } |
||
221 | else{ |
||
222 | phi=atan(b/a); |
||
223 | } |
||
224 | |||
225 | if (phi<0) phi = phi +_PI; |
||
226 | |||
227 | dis=pos1.x()*cos(phi) - pos1.y()*sin(phi); |
||
228 | //dis=ComputeSinogramS(pos1.x(), pos1.y(), pos2.x(), pos2.x()); |
||
229 | // control; transaxial FOV |
||
230 | if ( fabs(dis) > m_RingDiameter*0.5 ) { |
||
231 | #ifndef GAMOS_NO_VERBOSE |
||
232 | if( m_Debug ) { |
||
233 | cout << "PETProjDataMgr::AddEvent:WARNING! Event out of bounds (Transaxial): x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl; |
||
234 | } |
||
235 | #endif |
||
236 | return; |
||
237 | } |
||
238 | |||
239 | dis = dis + m_RingDiameter*0.5; |
||
240 | |||
241 | // discret values: |
||
242 | phi_i=(int)round( (double)(m_NOfAngles-1)*phi/_PI ); |
||
243 | dis_i=(int)round( (double)(m_NOfBins-1)*dis/(double)m_RingDiameter ); |
||
244 | |||
245 | if ((phi_i>=m_NOfAngles) || (dis_i>=m_NOfBins)) return; // only possible "=" because 'round' check it.. |
||
246 | |||
247 | // OLD: (SRRB included) sino3D[dis_i][phi_i][ (z1_i+z2_i)+ring_diff*(m_NOfPlanes-1) ]++; |
||
248 | |||
249 | int Zpos; |
||
250 | |||
251 | if (m_OutFormat==0) { |
||
252 | Zpos = (z1_i*m_NOfPlanes + z2_i); |
||
253 | } |
||
254 | else{ |
||
255 | |||
256 | if (z1_i>=z2_i) { // SIN Max Ring_Diff: Zpos= ( ((m_NOfPlanes-ring_diff)*(m_NOfPlanes-1-ring_diff))/2 + z2_i ); |
||
257 | |||
258 | Zpos= ( ((2*m_NOfPlanes-1 - m_MaxRingDifference - ring_diff)*(m_MaxRingDifference - ring_diff))/2 + z2_i); |
||
259 | |||
260 | }else{ |
||
261 | Zpos= ( (m_TotalAxialPlanes) - ((2*m_NOfPlanes-1 - m_MaxRingDifference - ring_diff +1)*(m_MaxRingDifference - ring_diff +1))/2 + z1_i ); |
||
262 | |||
263 | } |
||
264 | } |
||
265 | |||
266 | m_projections[dis_i][phi_i][ Zpos ]++; |
||
267 | m_TotalProjectionCoincidences++; |
||
268 | |||
269 | #ifndef GAMOS_NO_VERBOSE |
||
270 | if( m_Debug >1) { |
||
271 | cout << "PETProjDataMgr::AddEvent: x1= " << pos1.x() << ", y1= " << pos1.y() << ", z1= " << pos1.z() << " ; x2= " << pos2.x() << ", y2= " << pos2.y() << ", z2= " << pos2.z() << endl; |
||
272 | cout << "PETProjDataMgr::AddEvent: Sinogram pos.: distance(s)= " << dis_i << ", angular view(phi)= " << phi_i << " ; Zpos=" << Zpos <<"; Segment (Ring diff.) = " << ring_diff << endl; |
||
273 | } |
||
274 | #endif |
||
275 | |||
276 | |||
277 | } |
||
278 | |||
279 | //----------------------------------------------------------------------- |
||
280 | PETProjDataMgr::~PETProjDataMgr() |
||
281 | { |
||
282 | int i,j; |
||
283 | |||
284 | for(i=0;i<m_NOfBins;i++){ |
||
285 | for(j=0;j<m_NOfAngles;j++){ |
||
286 | free(m_projections[i][j]); |
||
287 | |||
288 | } |
||
289 | free(m_projections[i]); |
||
290 | |||
291 | } |
||
292 | free(m_projections); |
||
293 | |||
294 | } |
||
295 | |||
296 | /* TO DO: call lm_to_sino3D program |
||
297 | //----------------------------------------------------------------------- |
||
298 | void PETIOMgr::ReadFile() |
||
299 | { |
||
300 | if( !theFileIn ) OpenFileIn(); |
||
301 | |||
302 | PETOutput po; |
||
303 | G4bool bEof; |
||
304 | for(;;) { |
||
305 | po = ReadEvent( bEof ); |
||
306 | // theFileIn->read(reinterpret_cast<char *>(&po),sizeof(PetOutput)); |
||
307 | if(bDumpCout) PrintEvent(" PETOutput: ", po, bCartesian); |
||
308 | if( bEof ) break; |
||
309 | } |
||
310 | } |
||
311 | //----------------------------------------------------------------------- |
||
312 | PETOutput PETProjDataMgr::ReadEvent( G4bool& bEof ) |
||
313 | { |
||
314 | if( theFileIn == 0 ){ |
||
315 | G4Exception("PETIOMgr::ReadEvent, file not opened, call OpenFileIn() first "); |
||
316 | } |
||
317 | |||
318 | PETOutput po; |
||
319 | fread (&po, sizeof(struct PETOutput),1,theFileIn); |
||
320 | if ( feof( theFileIn ) ) { |
||
321 | bEof = TRUE; |
||
322 | } else { |
||
323 | bEof = FALSE; |
||
324 | } |
||
325 | |||
326 | return po; |
||
327 | |||
328 | } |
||
329 | */ |
||
330 | |||
331 | |||
332 | //----------------------------------------------------------------------- |
||
333 | void PETProjDataMgr::WriteInterfile() |
||
334 | { |
||
335 | |||
336 | char name_hv[512]; |
||
337 | char name_v[512]; |
||
338 | |||
339 | if (m_OutFormat==0){ |
||
340 | |||
341 | strcpy(name_hv, m_Filename); |
||
342 | strcpy(name_v,m_Filename); |
||
343 | strcat(name_hv, ".hv"); |
||
344 | strcat(name_v, ".v"); |
||
345 | |||
346 | fp=fopen(name_hv, "w"); |
||
347 | |||
348 | fprintf (fp, "!INTERFILE := \n"); |
||
349 | fprintf (fp, "name of data file := %s\n", name_v); |
||
350 | fprintf (fp, "!GENERAL DATA := \n"); |
||
351 | fprintf (fp, "!GENERAL IMAGE DATA :=\n"); |
||
352 | fprintf (fp, "!type of data := tomographic\n"); |
||
353 | fprintf (fp, "!version of keys := 3.3\n"); |
||
354 | fprintf (fp, "!data offset in bytes := 0\n"); |
||
355 | fprintf (fp, "imagedata byte order := littleendian\n"); |
||
356 | fprintf (fp, "!PET STUDY (General) :=\n"); |
||
357 | fprintf (fp, "!PET data type := 3D-Sinogram\n"); |
||
358 | fprintf (fp, "process status := Reconstructed\n"); |
||
359 | fprintf (fp, "!number format := unsigned short\n"); |
||
360 | fprintf (fp, "!number of bytes per pixel := 2\n"); |
||
361 | fprintf (fp, "number of dimensions := 3\n"); |
||
362 | fprintf (fp, "matrix axis label [1] := x\n"); |
||
363 | fprintf (fp, "!matrix size [1] := %i\n",m_NOfBins); |
||
364 | fprintf (fp, "scaling factor (mm/pixel) [1] := %f\n",(float)(m_RingDiameter/(m_NOfBins-1.0))); |
||
365 | |||
366 | fprintf (fp, "matrix axis label [2] := y\n"); |
||
367 | fprintf (fp, "!matrix size [2] := %i\n",m_NOfAngles); |
||
368 | |||
369 | fprintf (fp, "scaling factor (degree/pixel) [2] := %f\n",(float)(360./m_NOfAngles)); |
||
370 | |||
371 | fprintf (fp, "matrix axis label [3] := z\n"); |
||
372 | fprintf (fp, "!matrix size [3] := %i\n",m_NOfPlanes*m_NOfPlanes); |
||
373 | fprintf (fp, "scaling factor (mm/pixel) [3] := %f\n",(float)(m_AxialDistance/(m_NOfPlanes-1.0))); |
||
374 | |||
375 | fprintf (fp, "number of slices := %i\n",m_NOfPlanes*m_NOfPlanes); |
||
376 | fprintf (fp, "number of time frames := 1\n"); |
||
377 | fprintf (fp, "image scaling factor[1] := 1\n"); |
||
378 | fprintf (fp, "data offset in bytes[1] := 0\n"); |
||
379 | fprintf (fp, "quantification units := 1\n"); |
||
380 | fprintf (fp, "!END OF INTERFILE := \n"); |
||
381 | |||
382 | fclose(fp); |
||
383 | //(size_t)(m_NOfBins*m_NOfAngles*m_NOfPlanes*m_NOfPlanes); |
||
384 | |||
385 | }else{ |
||
386 | |||
387 | strcpy(name_hv, m_Filename); |
||
388 | strcpy(name_v,m_Filename); |
||
389 | |||
390 | strcat(name_hv, ".hs"); // STIR extension: .hs .s |
||
391 | strcat(name_v, ".s"); |
||
392 | |||
393 | fp=fopen(name_hv, "w"); |
||
394 | |||
395 | fprintf (fp, "!INTERFILE := \n"); |
||
396 | fprintf (fp, "name of data file := %s\n",name_v); |
||
397 | fprintf (fp, "!GENERAL DATA := \n"); |
||
398 | fprintf (fp, "!GENERAL IMAGE DATA :=\n"); |
||
399 | fprintf (fp, "!type of data := PET\n"); |
||
400 | // fprintf (fp, "!version of keys := 3.3\n"); STIR format is not 3.3 (almost but not completely), ERROR in STIR if it is not removed |
||
401 | // fprintf (fp, "!data offset in bytes := 0\n"); |
||
402 | fprintf (fp, "imagedata byte order := littleendian\n"); |
||
403 | fprintf (fp, "!PET STUDY (General) :=\n"); |
||
404 | fprintf (fp, "!PET data type := Emission\n"); |
||
405 | fprintf (fp, "applied corrections := {arc correction}\n"); // {none}\n"); |
||
406 | // fprintf (fp, "process status := Reconstructed\n"); |
||
407 | fprintf (fp, "!number format := unsigned integer\n"); |
||
408 | fprintf (fp, "!number of bytes per pixel := 2\n"); |
||
409 | |||
410 | fprintf (fp, "number of dimensions := 4\n"); |
||
411 | fprintf (fp, "matrix axis label [4] := segment\n"); |
||
412 | fprintf (fp, "!matrix size [4] := %i\n",m_MaxRingDifference*2 + 1); |
||
413 | // fprintf (fp, "scaling factor (mm/pixel) [1] := %f\n",(float)(d_FOV/(m_NOfBins-1))); |
||
414 | fprintf (fp, "matrix axis label [3] := axial coordinate\n"); |
||
415 | fprintf (fp, "!matrix size [3] := { "); |
||
416 | if (m_MaxRingDifference==0) |
||
417 | { |
||
418 | fprintf (fp, "%i}\n", m_NOfPlanes); |
||
419 | }else{ |
||
420 | for(int m=m_NOfPlanes-m_MaxRingDifference;m<=m_NOfPlanes;m++) fprintf (fp, "%i,", m); |
||
421 | for(int m=m_NOfPlanes-1;m>m_NOfPlanes-m_MaxRingDifference;m--) fprintf (fp, "%i,", m); |
||
422 | fprintf (fp, "%i}\n", m_NOfPlanes-m_MaxRingDifference); |
||
423 | } |
||
424 | fprintf (fp, "matrix axis label [2] := view\n"); |
||
425 | fprintf (fp, "!matrix size [2] := %i\n",m_NOfAngles); |
||
426 | fprintf (fp, "matrix axis label [1] := tangential coordinate\n"); |
||
427 | fprintf (fp, "!matrix size [1] := %i\n",m_NOfBins); |
||
428 | |||
429 | fprintf (fp, "minimum ring difference per segment := {"); // TO DO : add SPAN (m_MaxRingDifferenceiff per seg. variable) |
||
430 | fprintf (fp, "%i", -m_MaxRingDifference); |
||
431 | for(int m=-m_MaxRingDifference+1;m<=m_MaxRingDifference;m++) fprintf (fp, ",%i", m); |
||
432 | fprintf (fp, "}\n"); |
||
433 | fprintf (fp, "maximum ring difference per segment := {"); // TO DO : add SPAN (m_MaxRingDifferenceiff per seg. variable) |
||
434 | fprintf (fp, "%i", -m_MaxRingDifference); |
||
435 | for(int m=-m_MaxRingDifference+1;m<=m_MaxRingDifference;m++) fprintf (fp, ",%i", m); |
||
436 | fprintf (fp, "}\n"); |
||
437 | |||
438 | fprintf (fp, "inner ring diameter (cm) := %f\n", m_RingDiameter/10); // STIR Required parameter, now assigned to FOV (not detectors) |
||
439 | fprintf (fp, "average depth of interaction (cm) := 0.0001\n"); |
||
440 | fprintf (fp, "default bin size (cm) := %f\n",0.1*((float)m_RingDiameter/((float)m_NOfBins-1.0)) ); |
||
441 | fprintf (fp, "number of rings := %i\n",m_NOfPlanes ); |
||
442 | fprintf (fp, "distance between rings (cm) := %f\n", 0.1*((float)m_AxialDistance/(float)(m_NOfPlanes-1)) ); // Axial pixel dimension |
||
443 | |||
444 | fprintf (fp, "number of detectors per ring := %i\n",m_NOfAngles*2 ); |
||
445 | // fprintf (fp, "number of slices := %i\n",m_NOfPlanes*m_NOfPlanes); |
||
446 | fprintf (fp, "number of time frames := 1\n"); |
||
447 | fprintf (fp, "image scaling factor[1] := 1\n"); |
||
448 | fprintf (fp, "data offset in bytes[1] := 0\n"); |
||
449 | fprintf (fp, "quantification units := 1\n"); |
||
450 | fprintf (fp, "!END OF INTERFILE := \n"); |
||
451 | |||
452 | fclose(fp); |
||
453 | |||
454 | } |
||
455 | m_Buffer = (SINO_TYPE*) malloc( m_NOfBins*m_NOfAngles*m_TotalAxialPlanes*sizeof(SINO_TYPE)); |
||
456 | |||
457 | long unsigned int cont=0; |
||
458 | int i,j,k; |
||
459 | |||
460 | for(k=0;k<m_TotalAxialPlanes;k++){ |
||
461 | for(j=0;j<m_NOfAngles;j++){ |
||
462 | for(i=0;i<m_NOfBins;i++){ |
||
463 | m_Buffer[cont]=m_projections[i][j][k]; |
||
464 | cont++; |
||
465 | } |
||
466 | } |
||
467 | } |
||
468 | |||
469 | fp=fopen(name_v, "wb"); |
||
470 | |||
471 | //cout << 4096*sizeof(SINO_TYPE) << endl; |
||
472 | int nb=fwrite(m_Buffer,1,m_NOfBins*m_NOfAngles*m_TotalAxialPlanes*sizeof(SINO_TYPE), fp); |
||
473 | fclose(fp); |
||
474 | |||
475 | #ifndef GAMOS_NO_VERBOSE |
||
476 | cout << "PETProjDataMgr::WriteInterfile: File name: " << m_Filename << endl; |
||
477 | cout << "PETProjDataMgr::WriteInterfile: Numer of bytes written: " << nb << endl; |
||
478 | cout << "PETProjDataMgr::WriteInterfile: Planes = " << m_NOfPlanes << "; bins = " << m_NOfBins << "; ang_views = " << m_NOfAngles << endl; |
||
479 | cout << "PETProjDataMgr::WriteInterfile: Dimensions (mm): Transaxial FOV = " << m_RingDiameter << "; Axial FOV = " << m_AxialDistance << " ; Transaxial_pix = " << m_RingDiameter/(m_NOfBins-1) <<"; Plane width = " << m_AxialDistance/(m_NOfPlanes-1) << endl; // Image Axial Pixel(ssrb) == 0.5*(Plane_Width); |
||
480 | cout << "... " << endl; |
||
481 | |||
482 | cout << "PETProjDataMgr::WriteInterfile: Total Coinci: " << m_TotalCoincidences << endl; |
||
483 | cout << "PETProjDataMgr::WriteInterfile: Sino3D Coinci: " << m_TotalProjectionCoincidences << endl; |
||
484 | #endif |
||
485 | |||
486 | } |
||
487 | |||
488 | |||
489 | TVector3 PETProjDataMgr::Hits2Digits(const TVector3 &r){ |
||
490 | if (!m_nch) return r; |
||
491 | float smear=0.5; |
||
492 | |||
493 | if (m_nch<0) smear=m_rnd->Rndm(); |
||
494 | |||
495 | double angle = TMath::ATan2(r.X(),r.Y()); // vrne kot med -pi in pi |
||
496 | if (angle<0) angle+=TMath::TwoPi(); |
||
497 | |||
498 | angle= (int(angle/TMath::TwoPi()*TMath::Abs(m_nch))+smear)*TMath::TwoPi()/TMath::Abs(m_nch); |
||
499 | //(m_rnd->Rndm()-0.5)*m_AxialDistance; |
||
500 | return TVector3(sin(angle), cos(angle),0); // z coordinata ni cisto v redu |
||
501 | |||
502 | } |
||
503 | |||
504 | int PETProjDataMgr::FwdProject(double x,double y, double z, int nmax, TH1*h){ |
||
505 | TVector3 r(x,y,z); |
||
506 | int h2d=h->InheritsFrom("TH2F"); |
||
507 | double tfac=m_RingDiameter*m_RingDiameter/4-r.Mag2(); |
||
508 | double rfac= m_AxialDistance/m_RingDiameter; |
||
509 | for (int i=0;i<nmax;i++){ |
||
510 | |||
511 | double phi= m_rnd->Rndm()*TMath::Pi(); |
||
512 | TVector3 s(1,0,0); |
||
513 | s.SetPhi(phi); |
||
514 | double sign = (m_rnd->Rndm()>0.5)? 1 : 0; |
||
515 | double theta = TMath::ACos(m_rnd->Rndm()*rfac); |
||
516 | theta+=sign*TMath::Pi(); |
||
517 | |||
518 | s.SetTheta(theta); |
||
519 | double t=r*s; |
||
520 | TVector3 rx=r-t*s; |
||
521 | |||
522 | double d=TMath::Sqrt(t*t+tfac); |
||
523 | |||
524 | TVector3 r1=rx+d*s; |
||
525 | TVector3 r2=rx-d*s; |
||
526 | |||
527 | //r1=Hits2Digits(r1); |
||
528 | //r2=Hits2Digits(r2); |
||
529 | |||
530 | AddEvent( r1 , r2); |
||
531 | if (h!=NULL){ |
||
532 | TVector3 s1=r2-r1; |
||
533 | double s1len= s1.Mag(); |
||
534 | int niter=int (100*s1len/m_RingDiameter); |
||
535 | for (int j=0;j<niter;j++){ |
||
536 | r2=r1+m_rnd->Rndm()*s1; |
||
537 | if (h2d) ((TH2F *) h)->Fill(r2.X(),r2.Y()); |
||
538 | else ((TH3F *) h)->Fill(r2.X(),r2.Y(),r2.Z()); |
||
539 | } |
||
540 | } |
||
541 | } |
||
542 | return 0; |
||
543 | } |
||
544 | |||
545 | int PETProjDataMgr::FwdProject(TH2F *img, TH2F *h){ |
||
546 | |||
547 | for (int i=0;i<img->GetNbinsX();i++) { |
||
548 | double x_=img->GetXaxis()->GetBinCenter( i+1 ); |
||
549 | for (int j=0;j<img->GetNbinsY();j++) { |
||
550 | double y_=img->GetYaxis()->GetBinCenter( j+1 ); |
||
551 | double density= img->GetBinContent(i+1,j+1); |
||
552 | if (density>0) FwdProject(x_,y_,m_AxialDistance*(m_rnd->Rndm()-0.5), density,h); |
||
553 | } |
||
554 | } |
||
555 | return 0; |
||
556 | } |
||
557 | |||
558 | int PETProjDataMgr::FwdProject(TH3F *img, TH3F *h){ |
||
559 | |||
560 | for (int i=0;i<img->GetNbinsX();i++) { |
||
561 | double x_=img->GetXaxis()->GetBinCenter( i+1 ); |
||
562 | for (int j=0;j<img->GetNbinsY();j++) { |
||
563 | double y_=img->GetYaxis()->GetBinCenter( j+1 ); |
||
564 | for (int k=0;k<img->GetNbinsZ();k++) { |
||
565 | double z_=img->GetZaxis()->GetBinCenter( k+1 ); |
||
566 | double density= img->GetBinContent(i+1,j+1,k+1); |
||
567 | if (density>0) FwdProject(x_,y_,z_, density,h); |
||
568 | } |
||
569 | } |
||
570 | } |
||
571 | return 0; |
||
572 | } |
||
573 | |||
574 | TH2F *PETProjDataMgr::Phantom(int kaj){ |
||
575 | TH2F *img= new TH2F("img","Original Image",100,-50,50,100,-50,50); |
||
576 | |||
577 | // izberi sliko 0: kroglice, 1: point source 2: central ball |
||
578 | switch (kaj){ |
||
579 | |||
580 | case 0: |
||
581 | for (int i=0;i<img->GetNbinsX();i++) { |
||
582 | for (int j=0;j<img->GetNbinsY();j++) { |
||
583 | double x_=img->GetXaxis()->GetBinCenter( i+1 ); |
||
584 | double y_=img->GetYaxis()->GetBinCenter( j+1 ); |
||
585 | double density=1000; |
||
586 | if ((x_*x_+y_*y_)<6) img->SetBinContent(i+1,j+1,density); |
||
587 | |||
588 | density=500; if ((x_-25)*(x_-25)+y_*y_<12) img->SetBinContent(i+1,j+1,density); |
||
589 | density=2000; if ((y_-25)*(y_-25)+x_*x_<2) img->SetBinContent(i+1,j+1,density); |
||
590 | } |
||
591 | } |
||
592 | break; |
||
593 | |||
594 | case 2: |
||
595 | for (int i=0;i<img->GetNbinsX();i++) { |
||
596 | for (int j=0;j<img->GetNbinsY();j++) { |
||
597 | double x_=img->GetXaxis()->GetBinCenter( i+1 ); |
||
598 | double y_=img->GetYaxis()->GetBinCenter( j+1 ); |
||
599 | double density=1000; |
||
600 | if ((x_*x_+y_*y_)<12.5) img->SetBinContent(i+1,j+1,density); |
||
601 | } |
||
602 | } |
||
603 | break; |
||
604 | |||
605 | case 1: |
||
606 | img->Fill(25,25,10000); |
||
607 | break; |
||
608 | |||
609 | } |
||
610 | |||
611 | return img; |
||
612 | } |