tools/phm2helix.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **   Name:          phm2helix.cpp
   5 **   Purpose:       Take projections of a phantom object
   6 **   Programmer:    Ian Kay
   7 **   Date Started:  Aug 2001
   8 **
   9 **  This is part of the CTSim program
  10 **  Copyright (C) 1983-2000 Kevin Rosenberg
  11 **
  12 **  $Id: phm2helix.cpp,v 1.1 2001/09/24 09:40:42 kevin Exp $
  13 **
  14 **  This program is free software; you can redistribute it and/or modify
  15 **  it under the terms of the GNU General Public License (version 2) as
  16 **  published by the Free Software Foundation.
  17 **
  18 **  This program is distributed in the hope that it will be useful,
  19 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  20 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21 **  GNU General Public License for more details.
  22 **
  23 **  You should have received a copy of the GNU General Public License
  24 **  along with this program; if not, write to the Free Software
  25 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  26 ******************************************************************************/
  27
  28 #include "ct.h"
  29 #include "timer.h"
  30
  31
  32 enum { O_PHANTOMPROG, O_PHMFILE,  O_DESC, O_NRAY, O_ROTANGLE, O_GEOMETRY, O_FOCAL_LENGTH, O_CENTER_DETECTOR_LENGTH,
  33   O_VIEW_RATIO, O_SCAN_RATIO, O_OFFSETVIEW,  O_TRACE, O_VERBOSE, O_HELP, O_DEBUG, O_VERSION };
  34
  35 static struct option phm2helix_options[] =
  36 {
  37   {"phantom", 1, 0, O_PHANTOMPROG},
  38   {"desc", 1, 0, O_DESC},
  39   {"nray", 1, 0, O_NRAY},
  40   {"rotangle", 1, 0, O_ROTANGLE},
  41   {"geometry", 1, 0, O_GEOMETRY},
  42   {"focal-length", 1, 0, O_FOCAL_LENGTH},
  43   {"center-detector-length", 1, 0, O_CENTER_DETECTOR_LENGTH},
  44   {"view-ratio", 1, 0, O_VIEW_RATIO},
  45   {"scan-ratio", 1, 0, O_SCAN_RATIO},
  46   {"offsetview", 1, 0, O_OFFSETVIEW},
  47   {"trace", 1, 0, O_TRACE},
  48   {"verbose", 0, 0, O_VERBOSE},
  49   {"help", 0, 0, O_HELP},
  50   {"debug", 0, 0, O_DEBUG},
  51   {"version", 0, 0, O_VERSION},
  52   {"phmfile", 1, 0, O_PHMFILE},
  53   {0, 0, 0, 0}
  54 };
  55
  56 static const char* g_szIdStr = "$Id: phm2helix.cpp,v 1.1 2001/09/24 09:40:42 kevin Exp $";
  57
  58
  59 void
  60 phm2helix_usage (const char *program)
  61 {
  62           std::cout << "usage: " << fileBasename(program) << " outfile ndet nview phmprog [OPTIONS]\n";
  63           std::cout << "Calculate (projections) through time varying phantom object  \n\n";
  64           std::cout << "     outfile          Name of output file for projectsions\n";
  65           std::cout << "     ndet             Number of detectors\n";
  66           std::cout << "     nview            Number of rotated views\n";
  67           std::cout << "     phmprog          Name of phm generation executable\n";
  68           std::cout << "     --phmfile        Temp phantom file name \n";
  69           std::cout << "     --desc           Description of raysum\n";
  70           std::cout << "     --nray           Number of rays per detector (default = 1)\n";
  71           std::cout << "     --rotangle       Angle to rotate view through (fraction of a circle)\n";
  72           std::cout << "                      (default = select appropriate for geometry)\n";
  73           std::cout << "     --geometry       Geometry of scanning\n";
  74           std::cout << "        parallel      Parallel scan beams (default)\n";
  75           std::cout << "        equilinear    Equilinear divergent scan beams\n";
  76           std::cout << "        equiangular   Equiangular divergent scan beams\n";
  77           std::cout << "     --focal-length   Focal length ratio (ratio to radius of phantom)\n";
  78           std::cout << "                      (default = 1)\n";
  79           std::cout << "     --view-ratio     Length to view (view diameter to phantom diameter)\n";
  80           std::cout << "                      (default = 1)\n";
  81           std::cout << "     --scan-ratio     Length to scan (scan diameter to view diameter)\n";
  82           std::cout << "                      (default = 1)\n";
  83           std::cout << "     --offsetview      Intial gantry offset in  'views' (default = 0)\n";
  84           std::cout << "     --trace          Trace level to use\n";
  85           std::cout << "        none          No tracing (default)\n";
  86           std::cout << "        console       Trace text level\n";
  87           std::cout << "        phantom       Trace phantom image\n";
  88           std::cout << "        proj          Trace projections\n";
  89           std::cout << "        plot          Trace plot\n";
  90           std::cout << "        clipping      Trace clipping\n";
  91           std::cout << "     --verbose        Verbose mode\n";
  92           std::cout << "     --debug          Debug mode\n";
  93           std::cout << "     --version        Print version\n";
  94           std::cout << "     --help           Print this help message\n";
  95 }
  96
  97
  98 int
  99 phm2helix_main (int argc, char* const argv[])
 100 {
 101         Phantom phm;
 102         std::string optGeometryName = Scanner::convertGeometryIDToName(Scanner::GEOMETRY_PARALLEL);
 103         char *opt_outfile = NULL;
 104         std::string opt_desc;
 105         std::string opt_PhmProg;
 106         std::string opt_PhmFileName = "tmpphmfile";
 107         int opt_ndet;
 108         int opt_nview;
 109         int opt_offsetview = 0;
 110         int opt_nray = 1;
 111         double dOptFocalLength = 2.;
 112         double dOptCenterDetectorLength = 2;
 113         double dOptViewRatio = 1.;
 114         double dOptScanRatio = 1.;
 115         int opt_trace = Trace::TRACE_NONE;
 116         int opt_verbose = 0;
 117         int opt_debug = 0;
 118         double opt_rotangle = -1;
 119         char* endptr = NULL;
 120         char* endstr;
 121
 122         Timer timerProgram;
 123
 124         while (1) {
 125                 int c = getopt_long(argc, argv, "", phm2helix_options, NULL);
 126
 127                 if (c == -1)
 128                 break;
 129
 130         switch (c) {
 131         case O_VERBOSE:
 132                 opt_verbose = 1;
 133                 break;
 134         case O_DEBUG:
 135                 opt_debug = 1;
 136                 break;
 137         case O_TRACE:
 138                 if ((opt_trace = Trace::convertTraceNameToID(optarg))
 139                                                                 == Trace::TRACE_INVALID) {
 140                         phm2helix_usage(argv[0]);
 141                         return (1);
 142                 }
 143                 break;
 144                           case O_PHMFILE:
 145                                 opt_PhmFileName = optarg;
 146                                 break;
 147                           case O_DESC:
 148                                 opt_desc = optarg;
 149                                 break;
 150                           case O_ROTANGLE:
 151                                 opt_rotangle = strtod(optarg, &endptr);
 152                                 endstr = optarg + strlen(optarg);
 153                                 if (endptr != endstr) {
 154                                         std::cerr << "Error setting --rotangle to " << optarg << std::endl;
 155                                         phm2helix_usage(argv[0]);
 156                                         return (1);
 157                                 }
 158                                 break;
 159                           case O_GEOMETRY:
 160                                 optGeometryName = optarg;
 161                                 break;
 162                           case O_FOCAL_LENGTH:
 163                                 dOptFocalLength = strtod(optarg, &endptr);
 164                                 endstr = optarg + strlen(optarg);
 165                                 if (endptr != endstr) {
 166                                         std::cerr << "Error setting --focal-length to " << optarg << std::endl;
 167                                         phm2helix_usage(argv[0]);
 168                                         return (1);
 169                                 }
 170                                 break;
 171                           case  O_CENTER_DETECTOR_LENGTH:
 172                                 dOptCenterDetectorLength = strtod(optarg, &endptr);
 173                                 endstr = optarg + strlen(optarg);
 174                                 if (endptr != endstr) {
 175                                         std::cerr << "Error setting --center-detector-length to " << optarg << std::endl;
 176                                         phm2helix_usage(argv[0]);
 177                                         return (1);
 178                                 }
 179                           break;
 180                           case O_VIEW_RATIO:
 181                                 dOptViewRatio = strtod(optarg, &endptr);
 182                                 endstr = optarg + strlen(optarg);
 183                                 if (endptr != endstr) {
 184                                         std::cerr << "Error setting --view-ratio to " << optarg << std::endl;
 185                                         phm2helix_usage(argv[0]);
 186                                         return (1);
 187                                 }
 188                                 break;
 189                           case O_SCAN_RATIO:
 190                                 dOptScanRatio = strtod(optarg, &endptr);
 191                                 endstr = optarg + strlen(optarg);
 192                                 if (endptr != endstr) {
 193                                         std::cerr << "Error setting --scan-ratio to " << optarg << std::endl;
 194                                         phm2helix_usage(argv[0]);
 195                                         return (1);
 196                                 }
 197                                 break;
 198                           case O_NRAY:
 199                                 opt_nray = strtol(optarg, &endptr, 10);
 200                                 endstr = optarg + strlen(optarg);
 201                                 if (endptr != endstr) {
 202                                   std::cerr << "Error setting --nray to %s" << optarg << std::endl;
 203                                   phm2helix_usage(argv[0]);
 204                                   return (1);
 205                                 }
 206                                 break;
 207                           case O_OFFSETVIEW:
 208                                 opt_offsetview = strtol(optarg, &endptr, 10);
 209                                 endstr = optarg + strlen(optarg);
 210                                 if (endptr != endstr) {
 211                                   std::cerr << "Error setting --offsetview to %s" << optarg << std::endl;
 212                                   phm2helix_usage(argv[0]);
 213                                   return (1);
 214                                 }
 215                                 break;
 216                           case O_VERSION:
 217 #ifdef VERSION
 218                                   std::cout << "Version: " << VERSION << std::endl << g_szIdStr << std::endl;
 219 #else
 220                                           std::cout << "Unknown version number\n";
 221 #endif
 222                                   return (0);
 223                           case O_HELP:
 224                           case '?':
 225                                 phm2helix_usage(argv[0]);
 226                                 return (0);
 227                           default:
 228                                 phm2helix_usage(argv[0]);
 229                                 return (1);
 230                         } // end of switch
 231           } // end of while loop
 232
 233           if (optind + 4 != argc) {
 234                 phm2helix_usage(argv[0]);
 235                 return (1);
 236           }
 237
 238           opt_outfile = argv[optind];
 239           opt_ndet = strtol(argv[optind+1], &endptr, 10);
 240           endstr = argv[optind+1] + strlen(argv[optind+1]);
 241           if (endptr != endstr) {
 242                 std::cerr << "Error setting --ndet to " << argv[optind+1] << std::endl;
 243                 phm2helix_usage(argv[0]);
 244                 return (1);
 245           }
 246           opt_nview = strtol(argv[optind+2], &endptr, 10);
 247           endstr = argv[optind+2] + strlen(argv[optind+2]);
 248           if (endptr != endstr) {
 249                 std::cerr << "Error setting --nview to " << argv[optind+2] << std::endl;
 250                 phm2helix_usage(argv[0]);
 251                 return (1);
 252           }
 253           opt_PhmProg = argv[optind+3];
 254
 255           if (opt_rotangle < 0) {
 256                 if (optGeometryName.compare ("parallel") == 0)
 257                   opt_rotangle = 0.5;
 258                 else
 259                   opt_rotangle = 1.0;
 260           }
 261
 262           std::ostringstream desc;
 263           desc << "phm2helix: NDet=" << opt_ndet
 264                    << ", Nview=" << opt_nview
 265                    << ", NRay=" << opt_nray
 266                    << ", RotAngle=" << opt_rotangle
 267                    << ", OffsetView =" << opt_offsetview
 268                    << ", Geometry=" << optGeometryName
 269                    << ", PhantomProg=" << opt_PhmProg
 270                    << ", PhmFileName=" << opt_PhmFileName;
 271           if (opt_desc.length()) {
 272                 desc << ": " << opt_desc;
 273           }
 274           opt_desc = desc.str();
 275
 276           opt_rotangle *= TWOPI;
 277
 278           int stat;
 279           char extcommand[100];
 280           if(opt_debug != 0)
 281                         std::cout  <<  opt_PhmProg  <<  " " << 0 << " " <<  opt_nview << " " << opt_PhmFileName  << std::endl;
 282            //extcommand <<  opt_PhmProg  <<  " " << 0 << " " <<  opt_nview << " " << opt_PhmFileName ;
 283
 284           sprintf(extcommand, "%s %d %d %s",    opt_PhmProg.c_str(), 0, opt_nview, opt_PhmFileName.c_str() );
 285
 286            stat = system( extcommand );
 287           if (stat != 0 )
 288                         std::cerr << "Error executing external phantom program " << opt_PhmProg << " with command " << extcommand << std::endl;
 289
 290           phm.createFromFile (opt_PhmFileName.c_str());
 291           remove(opt_PhmFileName.c_str());
 292
 293           Scanner scanner (phm, optGeometryName.c_str(), opt_ndet, opt_nview,
 294                                 opt_offsetview, opt_nray, opt_rotangle, dOptFocalLength,
 295                                 dOptCenterDetectorLength, dOptViewRatio, dOptScanRatio);
 296           if (scanner.fail()) {
 297                  std::cout << "Scanner Creation Error: " << scanner.failMessage()
 298                                                         << std::endl;
 299                  return (1);
 300           }
 301
 302           Projections pjGlobal(scanner);
 303
 304
 305           for( int iView = 0; iView < opt_nview; iView++ ){
 306                 if(opt_debug != 0)
 307                         std::cout  <<  opt_PhmProg  <<  " " << iView << " " <<  opt_nview << " " << opt_PhmFileName  << std::endl;
 308            //extcommand <<  opt_PhmProg  <<  " " << iView << " " <<  opt_nview << " " << opt_PhmFileName ;
 309
 310                         sprintf(extcommand, "%s %d %d %s",      opt_PhmProg.c_str(), iView, opt_nview, opt_PhmFileName.c_str() );
 311                 stat = system( extcommand );
 312
 313                 if (stat != 0 )
 314                         std::cerr << "Error executing external phantom program " << opt_PhmProg << " with command " << extcommand << std::endl;
 315                 Phantom phmtmp;
 316                 phmtmp.createFromFile (opt_PhmFileName.c_str());
 317
 318                 scanner.collectProjections (pjGlobal, phmtmp, iView,
 319                           1, scanner.offsetView(), true, opt_trace);
 320                 remove(opt_PhmFileName.c_str());
 321           }
 322
 323
 324           pjGlobal.setCalcTime (timerProgram.timerEnd());
 325                 pjGlobal.setRemark (opt_desc);
 326                 pjGlobal.write (opt_outfile);
 327                 if (opt_verbose) {
 328                   phm.print (std::cout);
 329                   std::cout << std::endl;
 330                 std::ostringstream os;
 331                 pjGlobal.printScanInfo (os);
 332                 std::cout << os.str() << std::endl;
 333                 std::cout << "  Remark: " << pjGlobal.remark() << std::endl;
 334                 std::cout << "Run time: " << pjGlobal.calcTime() << " seconds\n";
 335           }
 336
 337           return (0);
 338 }
 339
 340
 341 #ifndef NO_MAIN
 342 int
 343 main (int argc, char* argv[])
 344 {
 345   int retval = 1;
 346
 347   try {
 348     retval = phm2helix_main(argc, argv);
 349 #if HAVE_DMALLOC
 350     //    dmalloc_shutdown();
 351 #endif
 352   } catch (exception e) {
 353     std::cerr << "Exception: " << e.what() << std::endl;
 354   } catch (...) {
 355     std::cerr << "Unknown exception\n";
 356   }
 357
 358   return (retval);
 359 }
 360 #endif
 361