src/snark/snark.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/snark.cpp $
   5  $LastChangedRevision: 177 $
   6  $Date: 2015-08-03 16:44:07 -0400 (Mon, 03 Aug 2015) $
   7  $Author: zye $
   8  ***********************************************************
   9
  10  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  11  *                                                               *
  12  *                         S N A R K   9 3                       *
  13  *                                                               *
  14  *                A PICTURE RECONSTRUCTION PROGRAM               *
  15  *                                                               *
  16  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  17
  18              WRITTEN- DECEMBER 1974
  19              STUART W. ROWLAND
  20              ROGER W. PERETTI
  21              DEPT. OF COMPUTER SCIENCE
  22              S.U.N.Y. AT BUFFALO
  23              4226 RIDGE LEA ROAD
  24              AMHERST, NEW YORK     14226
  25              TEL. (716) 831-1351
  26
  27 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  28
  29              REVISED BY-  G. T. HERMAN
  30                           K.J. CHEN
  31              SUMMER 1975
  32
  33 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  34
  35    SNARK76 - PRELIMINARY VERSION - SPRING 1976
  36 THIS VERSION OF SNARK INCORPORATES A NUMBER OF IMPROVEMENTS TO
  37 THE SNARK75 PROGRAM.  THE MOST SIGNIFICANT OF WHICH ARE-
  38    1) THE ABILITY TO GENERATE PROJECTION DATA SIMULATING A
  39       POLYCHROMATIC X-RAY SOURCE.
  40    2) THE INTRODUCTION OF AN OVERLAY STRUCTURE GREATLY REDUCING
  41       THE MEMORY REQUIREMENTS OF SNARK.  ALSO, CREATR AND SKUNK
  42       HAVE BEEN MODIFIED SO THAT THEY ARE CONSIDERABLY SMALLER.
  43    3) COMMENTS AND MESSAGES HAVE BEEN MODIFIED SO THAT THEY ARE
  44       MORE UNDERSTANDABLE.
  45    4) NON-STANDARD SUBSCRIPTING HAS BEEN REMOVED.
  46    5) ONLY THE RECONSTRUCTION FILE (RECFIL) NEED BE PRESENT DURING
  47       POST-PROCESSING FOR MOST OPERATIONS.  THE EXCEPTION IS THE
  48       DISTANCE COMMAND WHICH REQUIRES THE PROJECTION FILE (PRJFIL)
  49       TO COMPUTE THE RESIDUAL.
  50    6) THE ABILITY TO SPECIFY THE RANDOM NUMBER GENERATOR SEED HAS
  51       BEEN ADDED TO CREATE.
  52    7) A SIMPLE PHOTON SCATTERING MODEL HAS BEEN ADDED TO CREATE.
  53 THESE CHANGES HAVE BEEN IMPLEMENTED BY-
  54 S. W. ROWLAND AND G. T. HERMAN.
  55
  56 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  57
  58    SNARK77 - AUTUMN 1977
  59 THIS MODIFICATION OF SNARK76 PROVIDES A SIMPLER AND MORE POWERFUL
  60 COMMAND LANGUAGE STRUCTURE.  THE MAJOR MODIFICATIONS ARE:
  61    1) THE ADDITION OF THE ROUTINES GETINT AND GETNUM TO EXPAND THE
  62       CAPABILITIES OF THE COMMAND LANGUAGE.  GETINT REPLACES THE
  63       INTEGER OPTION WHICH WAS PART OF THE OLD GETWRD.  GETNUM
  64       PROVIDES A FREE FORM FLOATING POINT INPUT FACILITY TO
  65       SNARK77.
  66    2) THE CREATE (AND THEREFORE PICTURE AND PROJECTION) INPUT
  67       SEQUENCE HAS BEEN MODIFIED SO THAT IT IS MORE UNDERSTANDABLE
  68       WITHOUT REFERENCING THE SNARK MANUAL.  ALL OF THE INFORMATIO
  69       USED TO CREATE A PHANTOM OR PROJECTION DATA IS RETAINED IN
  70       THE GENERATED DATA SET HEADERS SO THAT THIS INFORMATION WILL
  71       NOT BE LOST TO FUTURE USERS OF THE DATA SET.
  72    3) A NAME NOW MUST BE ASSIGNED TO EACH ALGORITHM EXECUTION.
  73       THIS NAME IS PRINTED ON ALL POST-PROCESSING OUTPUT TO IDENTI
  74       THE PARTICULAR RECONSTRUCTION BEING USED.  IT IS SUGGESTED
  75       THAT THIS CARD INCLUDE AN INDICATION OF THE PARAMETERS USED
  76       BY THE ALGORITHM.
  77    4) ALMOST ALL OF THE SNARK CONTROL CARD FORMATS HAVE CHANGED.
  78    5) A CAPABILITY TO GENERATE A PARALLEL PROJECTION DATA SET FROM
  79       A DIVERGENT PROJECTION DATA SET VIA THE REBINNING FACILITY
  80       OF RDPROJ HAS BEEN ADDED.
  81    6) THE LINES COMMAND HAS BEEN EXTENDED TO GENERATE THE OUTPUT
  82       ON A DIGITAL PLOTTER.
  83
  84    THE OVERLAY STATEMENTS FROM THE CDC3500 VERSION OF SNARK76 HAVE
  85 BEEN REMOVED FOR THIS BURROUGHS B1700 VERSION.  SINCE THE SIZE OF
  86 THE PROGRAM HAS INCREASED TREMENDOUSLY, OVERLAYING WILL HAVE TO BE
  87 INSTALLED ON NON-PAGING COMPUTERS.  THE FOLLOWING TWO-LEVEL
  88 STRUCTURE IS SUGGESTED:
  89    (0,0)- THE MAIN PROGRAM AND COMMON SUBROUTINES.
  90    (1,0)- THE CREATER SUBSYSTEM.
  91    (2,0)- RDPICT AND RDPROJ.
  92    (3,0)- EXALG.
  93    (3,N)- EACH OF THE RECONSTRUCTION ALGORITHMS.
  94    (4,0)- THE POST-PROCESSING ROUTINES.
  95 EVEN WITH THIS STRUCTURE, IT MAY BE THAT SOME OF THE LARGER
  96 ALGORITHMS (QUAD, FOR EXAMPLE) MAY NOT FIT IN YOUR LIMITED
  97 MEMORY AND SHOULD THEREFORE BE DELETED FROM THE SYSTEM.
  98
  99 THESE CHANGES HAVE BEEN IMPLEMENTED BY-
 100 S. W. ROWLAND AND THE MEDICAL IMAGE PROCESSING GROUP (MIPG).
 101
 102 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 103
 104    SNARK89 Version 1.01 - July 1990
 105 The calculation of totden and totlen in the STRIP case in rdproj.f
 106 has been changed to agree with the manual.  Checks of nelem and
 107 pixsiz have been inserted in rdpict.f.
 108 These changes were made by G.T. Herman and D. Odhner.
 109
 110 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 111
 112 SNARK05
 113
 114 The major advances that are incorporated into SNARK05 are:
 115
 116 1. SNARK05 is implemented in C++;
 117 2. the file structures for holding projection data, phantoms and
 118    reconstructions have been redesigned to include XML headers;
 119 3. all itertive algorithms are now capable of handling image
 120    representations that use "blobs" as well as "pixels";
 121 4. the efficient data access ordering is a standard feature;
 122 5. all artificial restrictions on sizes have been removed - the only
 123    remaining restrictions are those imposed by the hardware, compiler
 124    and operating system.
 125
 126 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 127
 128 MAIN PROGRAM FOR SNARK SYSTEM
 129 THIS ROUTINE READS THE SNARK CONTROL CARDS, INTERPRETS THE MAJOR
 130 COMMAND ON THE CARD, AND INVOKES THAT ROUTINE.
 131
 132 */
 133
 134 #include <cstdio>
 135 #include <string>
 136 #include <sys/types.h>
 137 #include <sys/stat.h>
 138 #include <fcntl.h>
 139 #include <unistd.h>
 140
 141 #include "blkdta.h"
 142 #include "uiod.h"
 143
 144 #include "basis.h"
 145 #include "second.h"
 146 #include "mode.h"
 147 #include "settra.h"
 148 #include "creatr.h"
 149 #include "rdpict.h"
 150 #include "rdproj.h"
 151 #include "select.h"
 152 #include "stopex.h"
 153 #include "exalg.h"
 154 #include "eval.h"
 155 #include "disply.h"
 156 #include "pnch.h"
 157 #include "doline.h"
 158 #include "halft.h"
 159 #include "settra.h"
 160 #include "infile.h"
 161 #include "recfile.h"
 162 #include "file11.h"
 163 #include "int2str.h"
 164 #include "objects.h"
 165 #include "superior.h"
 166
 167 /*
 168   "JUST THE PLACE FOR A SNARK!' THE BELLMAN CRIED,
 169        AS HE LANDED HIS CREW WITH CARE;
 170   SUPPORTING EACH MAN ON THE TOP OF THE TIDE
 171        BY A FINGER ENTWINED IN HIS HAIR.
 172   "JUST THE PLACE FOR A SNARK!  I HAVE SAID IT TWICE:
 173        THAT ALONE SHOULD ENCOURAGE THE CREW.
 174   JUST THE PLACE FOR A SNARK!  I HAVE SAID IT THRICE:
 175        WHAT I TELL YOU THREE TIMES IS TRUE."
 176   FIT THE FIRST - THE LANDING
 177   THE HUNTING OF THE SNARK (AN AGONY IN EIGHT FITS)
 178   BY LEWIS CARROLL
 179 */
 180
 181 #define NUMBER_OF_COMMANDS 16
 182
 183 enum SnarkCommands
 184 {
 185   SNARK_COMMAND_MODE,
 186   SNARK_COMMAND_TRACE,
 187   SNARK_COMMAND_END,
 188   SNARK_COMMAND_CREATE,
 189   SNARK_COMMAND_PICTURE,
 190   SNARK_COMMAND_PROJECTION,
 191   SNARK_COMMAND_BASSIS,
 192   SNARK_COMMAND_SELECT,
 193   SNARK_COMMAND_STOP,
 194   SNARK_COMMAND_EXECUTE,
 195   SNARK_COMMAND_EVALUATE,
 196   SNARK_COMMAND_DISPLAY,
 197   SNARK_COMMAND_PUNCH,
 198   SNARK_COMMAND_LINES,
 199   SNARK_COMMAND_SKUNK,
 200   SNARK_COMMAND_SUPERIORIZE
 201 };
 202
 203 static const INTEGER main_codes[NUMBER_OF_COMMANDS] =
 204 {
 205   CHAR2INT('m','o','d','e'),
 206   CHAR2INT('t','r','a','c'),
 207   CHAR2INT('e','n','d',' '),
 208   CHAR2INT('c','r','e','a'),
 209   CHAR2INT('p','i','c','t'),
 210   CHAR2INT('p','r','o','j'),
 211   CHAR2INT('b','a','s','i'),
 212   CHAR2INT('s','e','l','e'),
 213   CHAR2INT('s','t','o','p'),
 214   CHAR2INT('e','x','e','c'),
 215   CHAR2INT('e','v','a','l'),
 216   CHAR2INT('d','i','s','p'),
 217   CHAR2INT('p','u','n','c'),
 218   CHAR2INT('l','i','n','e'),
 219   CHAR2INT('s','k','u','n'),
 220   CHAR2INT('s','u','p','e')
 221 };
 222
 223 static const char* kLockFilename = "snark.lock";
 224 static int gLockFd = -1;
 225
 226 int FindKeywordCode(int word)
 227 {
 228   for(int i = 0; i < NUMBER_OF_COMMANDS; i++) {       // find keyword code
 229     if(word == main_codes[i]) {
 230       return i;
 231     }
 232   }
 233   return -1;
 234 }
 235
 236 INTEGER GetCommandPhase(INTEGER code)
 237 {
 238   static const INTEGER Phase[NUMBER_OF_COMMANDS] =
 239   {
 240 //    0, // *     Comment
 241     0, // mode  Mode
 242     0, // trac  Trace
 243     0, // end   End
 244     1, // crea  Create
 245     2, // pict  Picture
 246     2, // proj  Projection
 247     3, // basi  Basis
 248     3, // sele  Select
 249     3, // stop  Stop
 250     3, // exec  Execue
 251     4, // eval  Evaluate
 252     4, // disp  Display
 253     4, // punc  Punch
 254     4, // line  Lines
 255     4, // skun  Skunk
 256     3  // supe  Superiorize
 257   };
 258
 259   if(code >= NUMBER_OF_COMMANDS) {
 260     return 4;
 261   }
 262
 263   return Phase[code];
 264 }
 265
 266 void atExitCleanup(void) {
 267 #ifndef __CYGWIN__
 268   if (gLockFd >= 0) {
 269     close(gLockFd);
 270     if (unlink(kLockFilename) != 0) {
 271       fprintf(stderr, "Error unlinking lock file\n");
 272     }
 273   }
 274 #endif
 275 }
 276
 277 int snark(int argc, char *argv[])
 278 {
 279   ///static const INTEGER syst = CHAR2INT('a','y','s','t');
 280   //static const INTEGER snar = CHAR2INT('s','n','a','r');
 281   //static const INTEGER rand_1 = CHAR2INT('r','a','n','d');
 282
 283   INTEGER word;
 284   BOOLEAN eol;
 285   //BOOLEAN flag;
 286   INTEGER lphase;  // last phase
 287   INTEGER phase;   // current phase
 288   BOOLEAN iptflg;   // picture command called
 289   BOOLEAN iprflg;   // projection command called
 290
 291   BOOLEAN islflg;   // select command called, by wei 1/05
 292
 293   REAL begin;
 294   REAL curr_time;
 295   REAL dur;
 296   int CommandCode;
 297
 298   // bug 179 - swr - 10/30/05
 299 #ifndef __CYGWIN__
 300   // bug 167 - swr - 9/24/05
 301   gLockFd = open(kLockFilename,O_RDWR|O_CREAT,0640);
 302   if (gLockFd < 0) {
 303     fprintf(stderr, "can't create lock file in current directory - check directory permissions\n");
 304     exit(1); /* can not open */
 305   }
 306   if (lockf(gLockFd,F_TLOCK,0)<0) {
 307     fprintf(stderr, "multiple snark executions in the same directory are not allowed!\n");
 308     exit(0); /* can not lock */
 309   }
 310   /* only first instance continues */
 311   char str[10];
 312   sprintf(str,"%d\n",getpid());
 313   write(gLockFd,str,strlen(str)); /* record pid to lockfile */
 314
 315   atexit(atExitCleanup);
 316
 317 #endif
 318
 319   // bug 187 - allow input filename on command line - swr - 11/11/05
 320   if (argc > 1) {
 321     close(fileno(stdin));
 322     if (open(argv[1], O_RDONLY) < 0) {
 323       perror(argv[1]);
 324       exit(1);
 325     }
 326   }
 327
 328   InFile.Open();
 329
 330   iptflg = FALSE;
 331   iprflg = FALSE;
 332   islflg = FALSE;     //initialization islflg, by wei 1/05
 333
 334   std::string progName = argv[0];    // bug89 - fix header - swr - 2/25/05
 335
 336   progName = progName.substr(progName.find_last_of('/')+1);
 337
 338   fprintf(output, "   %s.s170710 - A PICTURE RECONSTRUCTION PROGRAM\n", progName.c_str());
 339
 340   lphase = 0;
 341
 342   second(&begin);
 343
 344   for(;;) {
 345
 346     InFile.getnxt(TRUE);        // get next line of input  // changed "FALSE" to "TRUE. Lajos, Jan 25, 2005
 347
 348     word = InFile.getwrd(FALSE, &eol, main_codes, NUMBER_OF_COMMANDS);   // check for the keyword
 349
 350     // added echoing. Lajos, Jan 25, 2005
 351     if(eol == TRUE) {
 352       InFile.echoline(TRUE);
 353       continue;   // got to next line if end of line
 354     } else InFile.echoline(FALSE);
 355
 356     // find keyword code
 357     if((CommandCode = FindKeywordCode(word)) < 0) {
 358       // VALID COMMAND NOT FOUND
 359
 360       fprintf(output, " **** command ignored");   // invalid command
 361       continue;                          // goto next line
 362
 363     }
 364
 365     // VALID COMMAND FOUND
 366
 367     phase = GetCommandPhase(CommandCode);
 368
 369     if(phase != 0) {                // check phase
 370       if(phase < lphase) {
 371
 372         // INCORRECT CONTROL CARD SEQUENCE
 373         fprintf(output, "\n **** command sequence error");
 374         fprintf(output, "\n **** program aborted\n");
 375         return -1;
 376       }
 377
 378       // if end of execs close recfile
 379       // if changing phase to 4 or end encoutnered before phase 4
 380       if((lphase < 4) && (phase == 4)) {
 381         RecFile.Close();
 382       }
 383
 384       lphase = phase;
 385     }
 386     else {
 387       if((CommandCode == SNARK_COMMAND_END) && (lphase < 4)) {
 388         RecFile.Close();
 389       }
 390     }
 391
 392     // INSERT OVERLAY CALL BASED ON LPHASE HERE
 393     // PASS I-LPHASE*5 TO THE OVERLAY TO CALL THE RIGHT FUNCTION
 394     // DELETE THE COMPUTED GO TO AND ITS ASSOCIATED CALLS
 395
 396     switch(CommandCode) {
 397
 398 //    case SNARK_COMMAND_COMENT: // comment
 399 //      break;
 400
 401     case SNARK_COMMAND_MODE: // ESTABLISH CONSTRAINT MODE
 402       mode();
 403       break;
 404
 405     case SNARK_COMMAND_TRACE: // CHANGE THE TRACE SWITCH
 406       settra();
 407       break;
 408
 409     case SNARK_COMMAND_END:
 410       ///fprintf(output, "\n         %10i storage units needed", maxsp);
 411       //fprintf(output, "\n         %10i storage units needed\n", 0);   // Obsolete
 412
 413       // free the dinamically allocated global array objects
 414       if(MAX_NUMBER_OF_OBJECTS>0) delete[] objects;
 415
 416       fprintf(output, "\n");
 417       return 0;
 418
 419     case SNARK_COMMAND_CREATE: // CREATE TEST PATTERN AND/OR PROJECTION DATA
 420       if(File11.Open(TRUE) != 0) {
 421         fprintf(output, "\n **** unable to open file11 for wrting");
 422         fprintf(output, "\n **** program aborted\n");
 423               return -1;
 424       }
 425       creatr();
 426       File11.Close();
 427       break;
 428
 429     case SNARK_COMMAND_PICTURE: // SET PICT CALLED FLAG
 430       File11.isOpen = FALSE;
 431       rdpict(); //  READ TEST PATTERN
 432       iptflg = TRUE;
 433       break;
 434
 435     case SNARK_COMMAND_PROJECTION: // BE SURE THAT PICTURE CMD HAS BEEN ENCOUNTERED
 436       // IF NOT  ERROR MSG  AND STOP
 437       if(iptflg) {
 438         rdproj(); // READ PROJECTION DATA
 439
 440         if (File11.isOpen) File11.Close();
 441         iprflg = TRUE;
 442       }
 443       else {
 444         ///write (output,1060)
 445         fprintf(output, "\n **** command picture is required before projection");
 446         fprintf(output, "\n **** program aborted\n");
 447         return -1;
 448       }
 449       break;
 450
 451     case SNARK_COMMAND_BASSIS: // ESTABLISH THE BASIS
 452       basis();
 453       break;
 454
 455     case SNARK_COMMAND_SELECT: // RAY SELECTION CONDITIONS FOR PICK
 456       select(FALSE);
 457        islflg = TRUE;                             //user specified select, by wei 1/05
 458       break;
 459
 460     case SNARK_COMMAND_STOP: // ESTABLISH ALGORITHM TERMINATING CONDITIONS
 461       stopex();
 462       break;
 463
 464     case SNARK_COMMAND_EXECUTE: // EXECUTE THE RECONSTRUCTION ALGORITHM
 465       if(iprflg) {
 466         if(!islflg) select(TRUE);
 467         islflg = TRUE;          //no select command, select default efficient , by wei 1/05
 468         exalg();
 469       }
 470       else {
 471         fprintf(output, "\n **** command projection is required before execute");
 472         fprintf(output, "\n **** program aborted\n");
 473         return -1;
 474       }
 475       break;
 476
 477     case SNARK_COMMAND_EVALUATE: // EVALUATION OF THE RECONSTRUCTIONS
 478       Eval.eval_1();
 479       break;
 480
 481     case SNARK_COMMAND_DISPLAY: // DIGITAL DISPLAY OF THE RECONSTRUCTIONS
 482       disply();
 483       break;
 484
 485     case SNARK_COMMAND_PUNCH: // PUNCH THE RECONSTRUCTIONS
 486       pnch();
 487       break;
 488
 489     case SNARK_COMMAND_LINES: // LINES DISPLAY OF THE DIFFERENCE
 490       doline();
 491       break;
 492
 493     case SNARK_COMMAND_SKUNK: // PSEUDO HALF TONE DISPLSY OF THE RECONSTRUCTIONS
 494       halft();
 495       break;
 496
 497     case SNARK_COMMAND_SUPERIORIZE: // ESTABLISH SUPERIORIZATION OPTIONS
 498       initSuperiorization();
 499       break;
 500     }
 501
 502     second(&curr_time);
 503     dur = curr_time - begin;
 504     begin = curr_time;
 505     fprintf(output, "\n         %10.3f seconds used for processing command %s\n", dur, int2str(word)); // changed precision to three digits - swr 1/21/06
 506   };
 507
 508   fprintf(output, "\n");
 509
 510   return 0;
 511 }
 512