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