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/bh_correction.cpp $
5 $LastChangedRevision: 118 $
6 $Date: 2014-07-09 14:22:29 -0400 (Wed, 09 Jul 2014) $
8 ***********************************************************
12 * The class bh_correction is responsible for performing beam hardening
13 * correction for polychromatic input data. Given the prjfil that was calculated
14 * based on weighted average of attenuations at different energy levels and
15 * corrections coefficients [1,2] it approximates monochromatic projection data
16 * which is then used by the reconstruction algorithms.
21 * [1] G. T. Herman, Image Reconstruction from Projections, 1980
31 #include <sys/types.h>
36 #include "bh_correction.h"
40 int bh_correction::readInputFilePoly()
43 BOOLEAN eol, eol1, eol2;
44 int i, j; //loop counter variables
46 //read in the number of iterations
47 noIter = InFile.getint(TRUE, &eol);
48 if (eol == TRUE) //too few values provided
50 fprintf(output, "\n **** ERROR: number of iterations missing\n"
51 " **** program aborted\n\n");
54 fprintf(output, "\n Number of iterations for beam hardening "
55 "correction set to r = %i\n\n", noIter);
57 //read in the degree of the polynomial that approximates monochromatic
58 //projections given polychromatic projections
59 noDegree = InFile.getint(FALSE, &eol);
60 if (eol == TRUE) //too few values provided
62 fprintf(output, "\n **** ERROR: the degree is missing\n"
63 " **** program aborted\n\n");
66 fprintf(output, "\n The degree of the polynomial for beam hardening "
67 "correction set to %i\n\n", noDegree);
69 //read in the coefficients of polynomial that approximates monochromatic
70 //projections given polychromatic projections
71 a = new REAL[noDegree + 1];
72 a[0] = InFile.getnum(TRUE, &eol);
73 if (eol == TRUE) //too few values provided
75 fprintf(output, "\n **** ERROR: too few polynomial coefficients \n"
76 " **** program aborted\n\n");
79 for (i = 1; i <= noDegree; i++)
81 a[i] = InFile.getnum(FALSE, &eol);
82 if (eol == TRUE) //too few values provided
84 fprintf(output, "\n **** ERROR: too few polynomial coefficients \n"
85 " **** program aborted\n\n");
92 //read in the number of energy levels
93 //this should be the same as Spctrm.nergy, otherwise exit
94 nerg = InFile.getint(TRUE, &eol);
95 if (eol == TRUE) //too few values provided
97 fprintf(output, "\n **** ERROR: the nerg is missing\n"
98 " **** program aborted\n\n");
102 if (nerg != Spctrm.nergy) //these two should be the same
104 fprintf(output, "\n **** ERROR: the value of nerg does not agree\n"
105 "with the nomber of energy levels written on file11\n"
106 " **** program aborted\n\n");
109 fprintf(output, "\n The number of energy levels is %i\n\n", nerg);
111 //read in the number of points for piecewise linear conversion of
112 //phantom data from energy[0] to other energy levels
113 noPointsConv = InFile.getint(TRUE, &eol);
114 if (eol == TRUE) //too few values provided
116 fprintf(output, "\n **** ERROR: the noPointsConv is missing\n"
117 " **** program aborted\n\n");
121 "\n The number of points for conversion of the polynomial "
122 "for beam hardening correction set to %i\n\n",
125 //increase noPointsConv by 1 since we need point (0,0)
127 //read in the points used for constant linear conversion of
128 //phantom data from energy[0] to other energy levels
129 for (j = 1; j < Spctrm.nergy; j++)
131 //skip j=0 since there is no conversion from the energy[0] to energy[0]
132 c[j] = new point[noPointsConv];
134 c[j][0].x = 0.0; //these are always 0.0
135 c[j][0].y = 0.0; //because vacum attenuates nothing at all energy levels
137 for (i = 1; i < noPointsConv; i++)
139 c[j][i].x = InFile.getnum(TRUE, &eol1);
140 c[j][i].y = InFile.getnum(FALSE, &eol2);
141 if (eol1 == TRUE || eol2 == TRUE) //too few values provided
143 fprintf(output, "\n **** ERROR: too few point for "
144 "energy level %i \n **** program aborted\n\n", j);
154 ///////////////////////////////////////////////////////////////////////////////
156 int bh_correction::correction(SnarkPrjFile * BHC_PrjFile)
158 INTEGER i, j, k; //loop counters
160 REAL x, y, u, den; //store elementary object descriptions
161 char command[256]; //store command line to be used in system() call
162 char filename[256]; //store file name created on the run
163 char line[256]; //store temporarly a line of a file
167 sprintf(command, "cp %s/q_p0 prjfil", tempDir);
168 if (system(command) != 0)
170 fprintf(output, "\n **** ERROR: could not execute command"
171 "%s\n **** program aborted\n\n", command);
177 //an input file for reconstruction part of iteration,
178 //uses existing prjfil to generate recfil
179 ofstream bhcInputFromPrjfil;
180 //an array of input files to be created based on the recfil created by
181 //reconstruction step during an iterative run
182 //[0] will be used for monochromatic input file
183 //[1] will be used for polychromatic input file
184 ofstream bhcInputFiles[2];
186 //a copy of the original input file for the main snark14 run
187 //it is needed here for information about projection generation, etc.
190 //a recfil from the itermediate reconstruction during the iterative process
191 DIGFileSnarkRec bhcRecfil;
193 unsigned int NoOfPixels, NoOfRecSets;
195 //a reconstruction image for the intermediate reconstruction during
196 //the iterative process
199 /* make sure that the nergy prjfil's are present,
200 * this projfil's contain the origila projection data corrected by multi-
201 * plying by the correction coefficient provided by the user in the input
205 /* create nergy input files that use the same structures as the original
206 * input file, use PROJECTION OLD and run convolution algorithm for
210 //create and open input files
211 //copy the original input file to the temp directory and open it for reading
212 sprintf(command, " cp file11 %s/inputFile.in ", tempDir);
213 if (system(command) != 0)
216 "\n **** ERROR: could not copy input file to inputFile.in\n "
217 " **** program aborted\n\n");
221 //open a new input file for writing
223 sprintf(filename, "%s/bhc_in.in", tempDir);
224 bhcInputFromPrjfil.open(filename);
225 if (!bhcInputFromPrjfil.is_open())
227 fprintf(output, "\n **** ERROR: could not open file %s from "
228 "temp directory\n **** program aborted\n\n", filename);
232 //create input files that use the existing prjfil's to perform the reconstruction
233 bhcInputFromPrjfil << "PICTURE RECONSTUCTION " << GeoPar.nelem << " "
234 << GeoPar.pixsiz << "\n";
235 bhcInputFromPrjfil << "PROJECTION OLD \n";
236 bhcInputFromPrjfil << "EXECUTE CONV \n";
237 bhcInputFromPrjfil << "monochromatic intermediate reconstruction\n";
238 bhcInputFromPrjfil << "hamming 0.54 2 \n";
239 bhcInputFromPrjfil << "END\n";
241 //close the input file
242 bhcInputFromPrjfil.close();
244 //allocalte needed memory
245 //bhcRecfils = new DIGFileSnarkRec[Spctrm.nergy];
246 //bhcRec_buffer = new REAL * [Spctrm.nergy];
247 NoOfPixels = GeoPar.nelem * GeoPar.nelem;
249 /* repeat noIter times
251 for (r = 0; r < noIter; r++)
253 sprintf(filename, "%s/inputFile.in", tempDir);
254 inputFile.open(filename);
255 if (!inputFile.is_open())
258 "\n **** ERROR: could not read the inputFile.in from "
259 "temp directory\n **** program aborted\n\n");
263 /* run snark14 and create new input files based on reconstruction images,
264 * then create monochromatic projections (nergy of them) and polychromatic
265 * projections in order to perform the next step of iterative data refinement
267 * correct the data in projfil's and prepare the input files for the next
268 * round of iterative data refinement step
271 //run snark14 to obtain reconstruction in recfil
272 //for (i = 0; i < Spctrm.nergy; i++) {
273 //rename the existing q_p%i to prjfil
274 sprintf(command, "cd %s\n cp q_p%i prjfil", tempDir, r);
275 if (system(command) != 0)
277 fprintf(output, "\n **** ERROR: could not execute command "
278 "%s\n **** program aborted\n\n", command);
281 //run snark14 with input file in"i".in
282 sprintf(command, "cd %s\n snark14 bhc_in.in >out.out", tempDir);
283 if (system(command) != 0)
285 fprintf(output, "\n **** ERROR: could not execute command "
286 "%s\n **** program aborted\n\n", command);
289 //rename the recfil, so it is not lost at the next run of snark14
290 sprintf(command, "cd %s\n mv recfil recfil%i", tempDir, r);
291 if (system(command) != 0)
293 fprintf(output, "\n **** ERROR: could not execute command "
294 "%s\n **** program aborted\n\n", command);
299 //open just created recfil's and prepare data for reading
301 sprintf(filename, "%s/recfil%i", tempDir, r);
302 if (bhcRecfil.Open(filename) != 0)
304 fprintf(output, "\n **** ERROR: could not open file "
305 "%s (r = %i)\n **** program aborted\n\n", filename, r);
309 // verify that there is only one reconstruction set (convolution
310 // produces only one)
312 if (bhcRecfil.GetNoOfRecSets(&NoOfRecSets) != 0)
315 "\n **** ERROR: could not read the number of reconstruction"
316 " sets from file %s\n **** program aborted\n\n",
320 if (NoOfRecSets != 1)
322 fprintf(output, "\n **** ERROR: the number (= %i) of reconstruction"
323 " sets from file %s is incorrect, should be =1\n "
324 "**** program aborted\n\n", NoOfRecSets, filename);
328 //allocate space for reconstuction data
329 bhcRecBuffer = new REAL[NoOfPixels];
331 //select the reconstruction
332 if (bhcRecfil.SelectRec(0) != 0)
334 fprintf(output, "\n **** ERROR: could not select reconstruction #0 "
336 "**** program aborted\n\n", filename);
340 //read in the data from the reconstruction pixel by pixel
341 bhcRecfil.GetRecData(bhcRecBuffer);
348 // create new set of input files for monochromatic projections (overwrite the old ones)
349 //with structures being pixel by pixel
350 //description based on the reconstruction images from recfil's above;
351 // create an input file for polychromatic projection with density
352 //values based on densities obtained in the reconstruction recfil's
353 //and distibution based on the original input file
355 //open the monochromatic input file for writing
356 //and the polychromatic input file for writing
357 for (i = 0; i < 2; i++)
359 sprintf(filename, "%s/%iin_rec%i.in", tempDir, r, i);
360 bhcInputFiles[i].open(filename);
361 if (!bhcInputFiles[i].is_open())
363 fprintf(output, "\n **** ERROR: could not open file %s from "
364 "temp directory\n **** program aborted\n\n", filename);
369 //start copying initial lines from the original input file to
370 //the new input files
372 // SPECTRUM POLYCHROMATIC imax
373 // d1 p1 d2 p2 ... d1imax pimax
375 for (i = 0; i < 2; i++)
377 bhcInputFiles[i] << "CREATE\n";
380 while (!inputFile.eof())
382 inputFile.getline(line, 255);
383 if (!strncasecmp(line, "SPEC", 4))
384 { //if line starts with the word "SPECTRUM"
385 break; //jump out of the loop
387 //otherwise write the line to all the newly created input files
388 for (i = 0; i < 2; i++)
390 bhcInputFiles[i] << line << "\n";
394 // SPECTRUM MONOCHROMATIC "Spctrm.energy[i]
395 // to the monochromatic input files
397 bhcInputFiles[0] << "SPECTRUM MONOCHROMATIC " << Spctrm.energy[0]
401 // SPECTRUM POLYCHROMATIC "Spctrm.nergy"
402 // Spctrm.energy[0] Spctrm.engwt[0] ... Spctrm.energy[nergy-1] Spctrm.engwt[nergy-1]
403 // to the polychromatic input file
404 bhcInputFiles[1] << "SPECTRUM POLYCHROMATIC " << Spctrm.nergy << "\n";
405 for (i = 0; i < Spctrm.nergy; i++)
407 bhcInputFiles[1] << Spctrm.energy[i] << " "
408 << 100 * Spctrm.engwt[i] << " ";
410 bhcInputFiles[1] << "\n";
413 // to the input files
414 for (i = 0; i < 2; i++)
416 bhcInputFiles[i] << "OBJECTS" << "\n";
419 //write the pixel by pixel description of the reconstruction image
420 // to the input files
421 INTEGER sizeHalf = 0.5 * GeoPar.nelem;
423 for (k = 0; k < NoOfPixels; k++)
425 x = (k % GeoPar.nelem - sizeHalf) * GeoPar.pixsiz;
426 y = (k / GeoPar.nelem - sizeHalf) * GeoPar.pixsiz;
427 u = GeoPar.pixsiz / 2;
428 den = bhcRecBuffer[k];
430 bhcInputFiles[0] << "RECT\t" << x << "\t" << y << "\t" << u << "\t"
431 << u << "\t" << 0 << "\t" << den << "\n";
433 bhcInputFiles[1] << "RECT\t" << x << "\t" << y << "\t" << u << "\t"
434 << u << "\t" << 0 << "\t";
436 bhcInputFiles[1] << den << "\nDENS";
438 //compute the densities using piecewise linear conversion function
439 for (i = 1; i < Spctrm.nergy; i++)
443 REAL m, mm; //slope and y intercept of the line segment
445 for (j = 1; j < noPointsConv; j++)
448 || (j == noPointsConv - 1 && den > c[i][j].x))
450 m = (c[i][j].y - c[i][j - 1].y)
451 / (c[i][j].x - c[i][j - 1].x);
452 mm = c[i][j].y - m * c[i][j].x;
453 nextden = m * den + mm;
456 if (den == c[i][j].x)
463 bhcInputFiles[1] << "\t" << nextden;
466 bhcInputFiles[1] << "\n";
469 //write the lines in input files up to line
471 for (i = 0; i < 2; i++)
473 bhcInputFiles[i] << "LAST 1 \n" << "PHANTOM AVERAGE "
474 << Creacm.nave1 << "\n" << GeoPar.nelem << " "
475 << GeoPar.pixsiz << "\n" << "RAYSUM \n" << "PICTURE TEST \n"
476 << "PROJECTION BHPSEUDO \n"
477 << "psudo run to create prjfil for correction\n";
480 //skip through lines of the original input file until line
483 while (!inputFile.eof())
485 inputFile.getline(line, 256);
486 if (!strncasecmp(line, "GEOM", 4))
487 { //if line starts with the word "GEOM"
488 break; //jump out of the loop
492 //copy the lines from the original input file to all the new
493 // input files until the line
496 for (i = 0; i < 2; i++)
498 bhcInputFiles[i] << line << "\n";
500 while (!inputFile.eof())
502 inputFile.getline(line, 255);
503 //write the line to all the newly created input files
504 for (i = 0; i < 2; i++)
506 bhcInputFiles[i] << line << "\n";
508 if (!strncasecmp(line, "BACK", 4))
509 { //if line starts with the word "BACKGROUND"
510 break; //jump out of the loop
514 //finish the files with END line and close the files
515 for (i = 0; i < 2; i++)
517 bhcInputFiles[i] << "END \n";
518 bhcInputFiles[i].close();
520 //close the original input file
523 //obtain projection data from the input files created above
525 for (i = 0; i < 2; i++)
528 //run snark14 with input file %rin"i".in
529 sprintf(command, "cd %s\n snark14 %iin_rec%i.in >out%i.out",
531 if (system(command) != 0)
533 fprintf(output, "\n **** ERROR: could not execute command"
534 "%s\n **** program aborted\n\n", command);
537 //rename the prjfil, so it is not lost at the next run of snark14
539 sprintf(command, "cd %s\n mv prjfil p%i", tempDir, r);
541 sprintf(command, "cd %s\n mv prjfil m%i", tempDir, r);
542 if (system(command) != 0)
544 fprintf(output, "\n **** ERROR: could not execute command"
545 "%s\n **** program aborted\n\n", command);
550 //correct the projection data following the formula
551 //m[i] = miu q[p][i] + (m'[i] - miu q(p')[i])
554 } //end of repeat noIter times
556 fprintf(output, "\n **** Finished the iterations\n\n");
558 /* write the final value of m[0] to prjfil that will be passed to the remainder
559 * of the run of snark14
562 sprintf(command, "cp %s/q_p%i prjfil", tempDir, noIter);
563 if (system(command) != 0)
565 fprintf(output, "\n **** ERROR: could not execute command"
566 "%s\n **** program aborted\n\n", command);
570 //if (bhcRecfil != NULL) delete [] bhcRecfil;
571 //if (bhc_InputFiles != NULL) delete [] bhc_InputFiles;
572 if (bhcRecBuffer != NULL)
573 delete[] bhcRecBuffer;
577 ///////////////////////////////////////////////////////////////////////////////
578 ///////////////////////////////////////////////////////////////////////////////
579 ///////////////////////////////////////////////////////////////////////////////
581 int bh_correction::idr_correction( INTEGER r)
583 //correct the projection data following the formula
584 //m[i] = miu q[p][i] + (m'[i] - miu q(p')[i])
585 REAL miu = 1.0; //the relaxation parameter, for now permanently set to 1.0
586 unsigned INTEGER i, j, k; //loop counters
587 DIGFileSnarkProj Org, NewM, NewP, Corrected;
590 REAL** Buffer = new REAL*[4]; //buffer to store a projection data
592 unsigned int NoOfRays;
593 unsigned int NoOfProjections;
594 ProjFileMH MainHeader;
595 ProjFileAH ApplicationHeader;
596 ProjFilePrjH ProjectionHeader[3];
598 //open original prjfil corrected by coefficient at energy level 0
600 //read Main Header and Application Header of the file
601 sprintf(filename, "%s/q_p0", tempDir);
602 if (Org.Open(filename) != 0)
604 fprintf(output, "\n **** ERROR: could not open file "
605 "%s (Org)\n **** program aborted\n\n", filename);
609 if (Org.GetNoOfRays(&NoOfRays) != 0)
612 "\n **** ERROR: could not read the number of rays from file "
613 "%s\n **** program aborted\n\n", filename);
617 if (Org.GetNoOfProjs(&NoOfProjections) != 0)
620 "\n **** ERROR: could not read the number of projections from file "
621 "%s\n **** program aborted\n\n", filename);
625 if (Org.GetMainHeader(&MainHeader) != 0)
627 fprintf(output, "\n **** ERROR: could not read Main Header of file "
628 "%s\n **** program aborted\n\n", filename);
632 if (Org.GetAppHeader(&ApplicationHeader) != 0)
635 "\n **** ERROR: could not read Application Header of file "
636 "%s\n **** program aborted\n\n", filename);
640 //open just created monochromatic prjfil at energy level i
642 sprintf(filename, "%s/m%i", tempDir, r);
643 if (NewM.Open(filename) != 0)
645 fprintf(output, "\n **** ERROR: could not open file "
646 "%s\n **** program aborted\n\n", filename);
650 //open just created polychromatic prjfil
652 sprintf(filename, "%s/p%i", tempDir, r);
653 if (NewP.Open(filename) != 0)
655 fprintf(output, "\n **** ERROR: could not open file "
656 "%s\n **** program aborted\n\n", filename);
660 //open prjfil that will contain corrected data,
662 //set its Main Header and Application Header
663 sprintf(filename, "%s/q_p%i", tempDir, r + 1);
664 if (Corrected.Open(filename, &MainHeader) != 0)
666 fprintf(output, "\n **** ERROR: could not open file "
667 "%s\n **** program aborted\n\n", filename);
671 if (Corrected.SetAppHeader(&ApplicationHeader) != 0)
674 "\n **** ERROR: could not set Application Header of file "
675 "%s\n **** program aborted\n\n", filename);
679 //allocate memory to store a single projection data
680 for (j = 0; j < 4; j++)
682 (Buffer[j]) = new REAL[NoOfRays];
685 for (j = 0; j < NoOfProjections; j++)
688 // array set selection
689 if (Org.SelectProj(j) != 0)
691 fprintf(output, "\n **** ERROR: could not select projection of file"
692 "Org\n **** program aborted\n\n");
695 //get projection header
696 if (Org.GetProjHeader(&ProjectionHeader[0]) != 0)
699 "\n **** ERROR: could not get projection header from file"
700 "Org\n **** program aborted\n\n");
703 //get projection data
704 if (Org.GetProjData(Buffer[0]) != 0)
707 "\n **** ERROR: could not get projection data from file"
708 "Org\n **** program aborted\n\n");
711 // array set selection
712 if (NewM.SelectProj(j) != 0)
714 fprintf(output, "\n **** ERROR: could not select projection of file"
715 "NewM\n **** program aborted\n\n");
719 //get projection data
720 if (NewM.GetProjData(Buffer[1]) != 0)
723 "\n **** ERROR: could not get projection data from file"
724 "NewM\n **** program aborted\n\n");
727 // array set selection
728 if (NewP.SelectProj(j) != 0)
730 fprintf(output, "\n **** ERROR: could not select projection of file"
731 "NewP\n **** program aborted\n\n");
735 //get projection data
736 if (NewP.GetProjData(Buffer[2]) != 0)
739 "\n **** ERROR: could not get projection data from file"
740 "NewP\n **** program aborted\n\n");
743 //calculate corrected value for each ray in the projection
745 for (k = 0; k < NoOfRays; k++)
748 REAL den = Buffer[2][k];
751 for (i = 0; i <= noDegree; i++)
753 newden += denden * a[i];
754 denden = denden * den;
757 Buffer[3][k] = miu * Buffer[0][k] + Buffer[1][k] - miu * newden;
760 //add corrected projection to the Corrected file
761 if (Corrected.AppendProj(&ProjectionHeader[0], Buffer[3]) != 0)
763 fprintf(output, "\n **** ERROR: could not append projection to file"
764 "Corrected\n **** program aborted\n\n");
774 //delete allocated memory
775 if (Buffer[0] != NULL)
777 if (Buffer[1] != NULL)
779 if (Buffer[2] != NULL)
781 if (Buffer[3] != NULL)