2 ***********************************************************
3 $SNARK_Header: S N A R K 1 4 - A PICTURE RECONSTRUCTION PROGRAM $
4 $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/projfile.cpp $
5 $LastChangedRevision: 118 $
6 $Date: 2014-07-09 14:22:29 -0400 (Wed, 09 Jul 2014) $
8 ***********************************************************
19 SnarkPrjFile ProjFile;
21 SnarkPrjFile::SnarkPrjFile()
26 INTEGER SnarkPrjFile::Open(const char* FileName)
29 if (DIGFileSnarkProj::Open(FileName) != 0)
38 INTEGER SnarkPrjFile::ReadHeader(char* name, Noise_class* Noise,
39 spctrm_class* Spct, anglst_class* Anglst)
42 double treal1, treal2;
52 strncpy(name, tname, 80);
55 /////////////////////////////
56 // get GeoPar's attributes //
57 /////////////////////////////
60 NumberOfProjections = GeoPar.prjnum = tuint;
63 GetSampling(&tDIGSamp);
64 GeoPar.pinc = tDIGSamp;
67 GetNoOfRays(&tDIGDim);
68 NumberOfRays = GeoPar.usrays = tDIGDim;
70 GeometryTypeEnum tGeomType;
71 GetGeometryType(&tGeomType);
73 DetectorTypeEnum tDetType;
74 ProjTypeEnum tProjType;
76 if (tGeomType == GT_DIVERGENT)
80 GeoPar.strip = false; // bug 207 - swr - 3/3/06
83 GetGeometry(&treal1, &treal2, &tDetType);
85 if (tDetType == DT_TANGENT)
96 GeoPar.radius = treal1;
104 GetGeometry(&tDetType, &tProjType);
106 if (tDetType == DT_UNIFORM)
117 if (tProjType == PT_STRIP)
124 GeoPar.strip = false;
129 ////////////////////////////
130 // get Noise's attributes //
131 ////////////////////////////
133 //bug221 start - swr - 03/03/07
134 if (GetNoiseQuantum(&treal1, &treal2, &tint) == 0)
136 Noise->setQuantumNoise(tint, treal1, treal2);
140 Noise->setQuantumNoise(0, 1.0, 1.0);
143 if (GetNoiseScatter(&treal1, &treal2) == 0)
145 Noise->setScatter(treal1, treal2);
148 if (GetNoiseAdditive(&treal1, &treal2) == 0)
150 Noise->setAdditiveNoise(treal1, treal2);
153 if (GetNoiseMultiplicative(&treal1, &treal2) == 0)
155 Noise->setMultiplicativeNoise(treal1, treal2);
157 //bug221 end - swr - 03/03/07
159 ///////////////////////////
160 // get Spct's attributes //
161 ///////////////////////////
163 GetNoOfEnergies(&tuint);
166 for (i = 0; i < tuint; i++)
168 GetSpectrum(i, &tint, &treal1, &treal2);
170 Spct->energy[i] = tint;
171 Spct->engwt[i] = treal1;
172 Spct->backgr[i] = treal2;
175 ///////////////////////////////////
176 // calculate GeoPar's attributes //
177 ///////////////////////////////////
179 GeoPar.snrays = GeoPar.numray(GeoPar.nelem, GeoPar.pixsiz);
181 if (GeoPar.usrays > GeoPar.snrays)
183 GeoPar.nrays = GeoPar.usrays;
186 { // bug194 - swr - 1/6/06
187 GeoPar.nrays = GeoPar.snrays;
191 GeoPar.midray = GeoPar.nrays / 2;
192 GeoPar.fsnray = (GeoPar.nrays - GeoPar.snrays) / 2;
193 GeoPar.lsnray = GeoPar.fsnray + GeoPar.snrays - 1;
194 GeoPar.fusray = (GeoPar.nrays - GeoPar.usrays) / 2;
195 GeoPar.lusray = GeoPar.fusray + GeoPar.usrays - 1;
197 //////////////////////////////
198 // get Anglst's attributes //
199 //////////////////////////////
201 Anglst->Init(GeoPar.prjnum);
203 for (i = 0; i < GeoPar.prjnum; i++)
206 if (SelectProj(i) != 0)
211 if (GetProjAngle(&(Anglst->bth[i])) != 0)
216 Anglst->bsin[i] = (REAL) sin(Anglst->bth[i]);
217 Anglst->bcos[i] = (REAL) cos(Anglst->bth[i]);
220 Anglst->pbase = new REAL[GeoPar.nrays];
231 INTEGER SnarkPrjFile::ReadProj(INTEGER np, REAL* data, unsigned INTEGER nsize)
238 if (SelectProj(np) != 0)
244 tArray = new REAL[tNum];
247 if (nsize <= NumberOfRays)
250 // NSIZE IS FEWER THAN OR EQUAL TO USRAYS
251 tSkip = (NumberOfRays - nsize) / 2;
252 for (i = 0; i < nsize; i++)
254 data[i] = tArray[i + tSkip];
260 // NSIZE IS GREATER THAN USRAYS
261 for (i = 0; i < nsize; i++)
266 tSkip = (nsize - NumberOfRays) / 2;
268 for (i = 0; i < NumberOfRays; i++)
270 data[i + tSkip] = tArray[i];
279 INTEGER SnarkPrjFile::Open(const char* FileName, INTEGER usrays, INTEGER prjnum)
281 NumberOfRays = usrays;
282 NumberOfProjections = prjnum;
284 AnglesBuff = new REAL[prjnum];
286 if (DIGFileSnarkProj::Open(FileName, "unknown", usrays, GeoPar.pinc,
287 "This projection file has been created by snark14") != 0)
295 INTEGER SnarkPrjFile::WriteHeader(char* name, Geometry_class* Geo,
296 Noise_class* Noise, spctrm_class* Spct)
304 /////////////////////////////
305 // set GeoPar's attributes //
306 /////////////////////////////
312 SetGeometry(Geo->radius, Geo->stod, DT_TANGENT);
316 SetGeometry(Geo->radius, Geo->stod, DT_ARC);
326 SetGeometry(DT_UNIFORM, PT_STRIP);
330 SetGeometry(DT_UNIFORM, PT_LINE);
337 SetGeometry(DT_VARIABLE, PT_STRIP);
341 SetGeometry(DT_VARIABLE, PT_LINE);
346 ////////////////////////////
347 // set Noise's attributes //
348 ////////////////////////////
350 if (Noise->quanin > 0)
352 SetNoiseQuantum(Noise->quanmn, Noise->quancm, Noise->quanin);
355 if (Noise->sctnfl == true)
357 SetNoiseScatter(Noise->sctnpk, Noise->sctnwd);
360 if (Noise->addnfl == true)
362 SetNoiseAdditive(Noise->addnmn, Noise->addnsd);
365 if (Noise->ultnfl == true)
367 SetNoiseMultiplicative(Noise->ultnmn, Noise->ultnsd);
370 ///////////////////////////
371 // set Spct's attributes //
372 ///////////////////////////
375 double* treal1 = new double[Spct->nergy];
376 double* treal2 = new double[Spct->nergy];
378 for (i = 0; i < Spct->nergy; i++)
380 treal1[i] = Spct->engwt[i];
381 treal2[i] = Spct->backgr[i];
384 SetSpectrum(Spct->nergy, Spct->energy, treal1, treal2);
392 INTEGER SnarkPrjFile::WriteAngles(REAL* bth)
398 // if angles already stored report error
404 // store angles into buffer
405 for (i = 0; i < NumberOfProjections; i++)
407 AnglesBuff[i] = bth[i];
414 INTEGER SnarkPrjFile::WriteProj(INTEGER np, REAL* iprjec)
419 fprintf(output, "\n **** PRJFile error: projection can not "
420 "be stored while angles are not given !!!\n");
424 if (AppendProj(np, AnglesBuff[np], "Projection created by snark14", iprjec)
433 INTEGER SnarkPrjFile::Close()
436 if (AnglesBuff != NULL)
438 delete[] AnglesBuff; // bug 92 - Lajos - 03/02/2005
444 if (DIGFileSnarkProj::Close() != 0)