src/snark/halft.cpp

   1 /*
   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/halft.cpp $
   5  $LastChangedRevision: 118 $
   6  $Date: 2014-07-09 14:22:29 -0400 (Wed, 09 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  POST-PROCESSING ROUTINE TO GENERATE PSEUDO-HALF-TONE DISPLAYS OF
  11  RECONSTRUCTIONS ON THE LINE PRINTER.
  12  */
  13
  14 #include <cstdlib>
  15 #include <cstdio>
  16 #include <cstring>
  17 #include <cmath>
  18
  19 #include "blkdta.h"
  20
  21 #include "uiod.h"
  22 #include "int2str.h"
  23
  24 #include "recfile.h"
  25 #include "infile.h"
  26
  27 #include "modefl.h"
  28 #include "halft.h"
  29
  30 void strtrim(char* out, const char* in)
  31 {
  32         unsigned int first = 0;
  33         //find first non white from the left
  34         while ((in[first] == ' ') || (in[first] == '\t'))
  35         {  // removed third term of the condition. lajos, Nov 18, 2004
  36                 first++;
  37         }
  38
  39         unsigned int last = strlen(in) - 1;
  40         //find first non white from the right
  41         while (last > first)
  42         {
  43                 if ((in[last] != ' ') && (in[last] != '\t'))
  44                 {
  45                         break;
  46                 }
  47                 last--;
  48         }
  49
  50         unsigned int i;
  51
  52         for (i = 0; i < last - first + 1; i++)
  53         {
  54                 out[i] = in[first + i];
  55         }
  56
  57         out[i] = 0;
  58 }
  59
  60 int SavePGM(char* filename, unsigned char* array, int xdim, int ydim,
  61                 char* coment)
  62 {
  63         FILE* fpout;
  64
  65         if ((fpout = fopen(filename, "w")) == NULL)
  66         {
  67                 return -1;
  68         }
  69
  70
  71         fprintf(fpout, "P5\n");
  72         fprintf(fpout, "#%s\n", coment);
  73         fprintf(fpout, "%d %d\n", xdim, ydim);
  74         fprintf(fpout, "255\n");
  75         fwrite(array, xdim * ydim, sizeof(unsigned char), fpout);
  76
  77         fclose(fpout);
  78
  79         return 0;
  80 }
  81
  82 void SaveImage(
  83 REAL* x,
  84 REAL amax,
  85 REAL amin,
  86 INTEGER nrow,
  87 INTEGER mcol,
  88 INTEGER nran,
  89 INTEGER mp,
  90 CHAR* filename,
  91 CHAR* coment)
  92 {
  93         int i;
  94
  95         INTEGER np;
  96         REAL range;
  97         REAL ot;
  98
  99         // allocate space for image buffer
 100         unsigned char* array = new unsigned char[nrow * mcol];
 101         // changed "CHAR" to "char". lajos Nov 17, 2004
 102
 103         if (array == NULL)
 104         {
 105                 fprintf(output, "\n **** memory allocation error");
 106                 fprintf(output, "\n **** program aborted\n");
 107                 exit(-1); // Error
 108         }
 109
 110         // FORM INTENSITY VALUES MEAN AND DIFFERENCE ARRAYS
 111
 112         np = mp;
 113
 114         range = amax - amin;
 115
 116         if (np == 2)
 117         {
 118                 fprintf(output, "\n          intensity levels printed\n");
 119         }
 120         else
 121         {
 122                 np = 1;
 123                 fprintf(output, "\n          amplitude levels printed\n");
 124         }
 125
 126         range = (REAL) 1.0 / range;
 127
 128         fprintf(output, "\n          columns %5i to %5i being printed\n", 1, mcol);
 129
 130         for (i = 0; i < nrow * mcol; i++)
 131         {
 132
 133                 // AMPLITUDE OR INTENSITY
 134
 135                 ot = MAX0((REAL) 0.0, MIN0((REAL) 1.0, ((x[i] - amin) * range)));
 136
 137                 if (np == 2)
 138                         ot = (REAL) sqrt(ot);
 139
 140                 array[i] = (unsigned char) (255 * ot);
 141                 // changed "(CHAR)" to "(unsigned char)". lajos Nov 18, 2004
 142
 143         }
 144
 145         SavePGM(filename, array, mcol, nrow, coment);
 146
 147         delete[] array;  // bug 92 - Lajos - 03/02/2005
 148
 149         return;
 150 }
 151
 152 void halft()
 153 {
 154         // was in post but it is realy local
 155         CHAR prjnam[81];
 156         CHAR phnnam[81];
 157         CHAR recnam[81];
 158
 159         // was in geo but it is realy local
 160         INTEGER nelem;
 161         INTEGER area;
 162         //REAL    pixsiz;
 163
 164         // bug 124 - added minimum and maximum - swr - 7/6/05
 165         static const INTEGER half_codes[] =
 166         {
 167                         CHAR2INT('p', 'h', 'a', 'n'), CHAR2INT('d', 'i', 'f', 'f'),
 168                         CHAR2INT('m', 'i', 'n', 'i'), CHAR2INT('m', 'a', 'x', 'i'),
 169                         CHAR2INT('a', 'm', 'p', 'l'), CHAR2INT('i', 'n', 't', 'e') };
 170
 171         int i;
 172
 173         INTEGER mp;
 174         INTEGER nran;
 175         REAL amin;
 176         REAL amax;
 177
 178         INTEGER word;
 179         INTEGER type;
 180         REAL* Phantom;
 181         REAL* Recon;
 182         INTEGER flag[51];
 183         CHAR algn[5];
 184         unsigned INTEGER count;
 185         unsigned INTEGER iter;
 186         BOOLEAN eol;
 187         BOOLEAN phan;
 188         BOOLEAN minfl;
 189         BOOLEAN maxfl;
 190
 191         CHAR coment[200];
 192         CHAR FileName[100];
 193
 194         type = 0;
 195         phan = FALSE;
 196         minfl = FALSE;
 197         maxfl = FALSE;
 198         mp = 1;
 199         nran = -1;
 200
 201         // read input line
 202         // bug 124 - rewrote input line processing - swr - 7/6/05
 203         int pos = 0;
 204         int length = sizeof(half_codes) / sizeof(half_codes[0]);
 205         word = InFile.getwrd(FALSE, &eol, &half_codes[pos], length - pos);
 206         while (!eol)
 207         {
 208                 switch (word)
 209                 {
 210                 case CHAR2INT('p', 'h', 'a', 'n'):
 211                         phan = TRUE;
 212                         pos = 2;
 213                         break;
 214                 case CHAR2INT('d', 'i', 'f', 'f'):
 215                         type = 1;
 216                         pos = 2;
 217                         break;
 218                 case CHAR2INT('m', 'i', 'n', 'i'):
 219                         minfl = TRUE;
 220                         amin = InFile.getnum(FALSE, &eol);
 221                         pos = 3;
 222                         break;
 223                 case CHAR2INT('m', 'a', 'x', 'i'):
 224                         maxfl = TRUE;
 225                         amax = InFile.getnum(FALSE, &eol);
 226                         pos = 4;
 227                         break;
 228                 case CHAR2INT('a', 'm', 'p', 'l'):
 229                         eol = TRUE;
 230                         break;
 231                 case CHAR2INT('i', 'n', 't', 'e'):
 232                         mp = 2;
 233                         eol = TRUE;
 234                         break;
 235                 }
 236                 if (!eol)
 237                         word = InFile.getwrd(FALSE, &eol, &half_codes[pos], length - pos);
 238         };
 239
 240         if (RecFile.Open("recfil") != 0)
 241         {
 242                 fprintf(output, "\n **** unable to open recfil");
 243                 fprintf(output, "\n **** SKUNK execution aborted\n");
 244                 return;
 245         }
 246
 247         if (RecFile.GetNelem(&nelem) != 0)
 248         {
 249                 ;
 250         }
 251
 252         if (RecFile.GetProjName(prjnam) != 0)
 253         {
 254                 ;
 255         }
 256
 257         area = nelem * nelem;
 258
 259         InFile.listit(flag);
 260
 261         // if phantom or differential read phantom
 262         if (phan || (type != 0))
 263         {
 264
 265                 Phantom = new REAL[area];
 266
 267                 if (RecFile.ReadPhan(phnnam, Phantom) != 0)
 268                 {
 269                         fprintf(output, "\n **** phantom not present");
 270                         fprintf(output, "\n **** SKUNK execution aborted\n");
 271                         delete[] Phantom;  // bug 92 - Lajos - 03/02/2005
 272                         return;
 273                 }
 274         }
 275
 276         if (phan)
 277         {
 278                 fprintf(output, "\n          skunk display of phantom\n");
 279
 280                 fprintf(output, "\n          phantom name:    %s", phnnam);
 281
 282                 // bug 124 - restored flag checks - swr - 7/5/05
 283                 if (!minfl)
 284                 {
 285                         minfl = TRUE;
 286                         amin = Phantom[0];
 287                         for (i = 1; i < area; i++)
 288                         {
 289                                 if (Phantom[i] < amin)
 290                                         amin = Phantom[i];
 291                         }
 292                 }
 293
 294                 if (!maxfl)
 295                 {
 296                         maxfl = TRUE;
 297                         amax = Phantom[0];
 298                         for (i = 1; i < area; i++)
 299                         {
 300                                 if (Phantom[i] > amax)
 301                                         amax = Phantom[i];
 302                         }
 303                 }
 304
 305                 if (amax <= amin)
 306                 {
 307                         fprintf(output,
 308                                         "\n **** maximal density is less or equal to minimal density");
 309                         fprintf(output, "\n **** SKUNK execution aborted\n");
 310                         delete[] Phantom;  // bug 92 - Lajos - 03/02/2005
 311                         return;
 312                 }
 313
 314                 sprintf(coment,
 315                                 "This is phantom was created by snark14: Phantom Name %s",
 316                                 phnnam);
 317
 318                 char PhantomName[41];
 319
 320                 strncpy(PhantomName, phnnam, 40);
 321                 PhantomName[40] = 0; // terminate string
 322
 323                 char PhantomName2[41];
 324
 325                 // trim whitespaces
 326                 strtrim(PhantomName2, PhantomName);
 327
 328                 sprintf(FileName, "%s.pgm", PhantomName2);
 329
 330                 // bug 169 - swr - 9/25/05
 331                 for (i = 0; i < 100; ++i)
 332                 {
 333                         if (FileName[i] == ' ')
 334                                 FileName[i] = '_';
 335                 }
 336
 337                 SaveImage(Phantom, amax, amin, nelem, nelem, nran, mp, FileName,
 338                                 coment);
 339         }
 340
 341         Recon = new REAL[area];
 342
 343         while (RecFile.ReadRec(recnam, algn, &count, &iter, Recon) == 0)
 344         {
 345
 346                 // continue if reconstruction not selected
 347                 if (flag[iter - 1] == 0)
 348                 {
 349                         continue;
 350                 }
 351
 352                 if (type == 0)
 353                 {
 354                         fprintf(output,
 355                                         "\n          skunk display of reconstruction using %s at iter %4i\n",
 356                                         algn, count);
 357                 }
 358                 else
 359                 {
 360                         fprintf(output,
 361                                         "\n          skunk difference picture using %s at iter %4i\n",
 362                                         algn, count);
 363
 364                         for (i = 0; i < area; i++)
 365                         {
 366                                 Recon[i] = (REAL) fabs(Phantom[i] - Recon[i]);
 367                         }
 368                 }
 369
 370                 fprintf(output, "\n          projection data: %s", prjnam);
 371
 372                 fprintf(output, "\n          execution name:  %s", recnam);
 373
 374                 // bug 124 - restored flag checks - swr - 7/5/05
 375                 if (!minfl)
 376                 {
 377                         minfl = TRUE;
 378                         amin = Recon[0];
 379                         for (i = 1; i < area; i++)
 380                         {
 381                                 if (Recon[i] < amin)
 382                                         amin = Recon[i];
 383                         }
 384                 }
 385
 386                 if (!maxfl)
 387                 {
 388                         maxfl = TRUE;
 389                         amax = Recon[0];
 390                         for (i = 1; i < area; i++)
 391                         {
 392                                 if (Recon[i] > amax)
 393                                         amax = Recon[i];
 394                         }
 395                 }
 396
 397                 if (amax <= amin)
 398                 {
 399                         fprintf(output,
 400                                         "\n **** maximal density is less or equal to minimal density");
 401                         fprintf(output, "\n **** SKUNK execution aborted\n");
 402                         break;
 403                 }
 404
 405                 char ProjectionName[13];
 406
 407                 strncpy(ProjectionName, prjnam, 12);
 408                 ProjectionName[12] = 0; // terminate string
 409
 410                 char ProjectionName2[13];
 411
 412                 // trim whitespaces
 413                 strtrim(ProjectionName2, ProjectionName);
 414
 415                 char ReconstructionName[13];
 416
 417                 strncpy(ReconstructionName, recnam, 12);
 418                 ReconstructionName[12] = 0; // terminate string
 419
 420                 char ReconstructionName2[13];
 421
 422                 // trim whitespaces
 423                 strtrim(ReconstructionName2, ReconstructionName);
 424
 425                 sprintf(FileName, "%s_%s_%4s_%04u_%c_%c.pgm", ProjectionName2,
 426                                 ReconstructionName2, algn, count, (type != 0) ? 'd' : 'r',
 427                                 (mp == 2) ? 'i' : 'a');
 428                 sprintf(coment,
 429                                 "This is reconstruction was created by snark14: Algorithm Name %4s Iteration %d",
 430                                 algn, count);
 431
 432                 // bug 169 - swr - 9/25/05
 433                 for (i = 0; i < 100; ++i)
 434                 {
 435                         if (FileName[i] == ' ')
 436                                 FileName[i] = '_';
 437                 }
 438
 439                 SaveImage(Recon, amax, amin, nelem, nelem, nran, mp, FileName, coment);
 440         }
 441
 442         delete[] Recon;  // bug 92 - Lajos - 03/02/2005
 443
 444         if (phan || (type != 0))
 445         {
 446                 delete[] Phantom;  // bug 92 - Lajos - 03/02/2005
 447         }
 448
 449         RecFile.Close();
 450
 451         fprintf(output, "\n");
 452         return;
 453 }