X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;ds=sidebyside;f=src%2Fphm2pj.cpp;fp=src%2Fphm2pj.cpp;h=fe831548437ad183353f8882e7d3f3479df79443;hb=f4a23943110823118f35756dd41fbd6707f04511;hp=0000000000000000000000000000000000000000;hpb=b5857e74e5735455b5ef11cbea5044ae7b2e8a0d;p=ctsim.git

diff --git a/src/phm2pj.cpp b/src/phm2pj.cpp
new file mode 100644
index 0000000..fe83154
--- /dev/null
+++ b/src/phm2pj.cpp
@@ -0,0 +1,342 @@
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+**   Name:          phm2pj.cpp
+**   Purpose:       Take projections of a phantom object
+**   Programmer:    Kevin Rosenberg
+**   Date Started:  1984
+**
+**  This is part of the CTSim program
+**  Copyright (C) 1983-2000 Kevin Rosenberg
+**
+**  $Id: phm2pj.cpp,v 1.1 2000/06/17 20:12:15 kevin Exp $
+**
+**  This program is free software; you can redistribute it and/or modify
+**  it under the terms of the GNU General Public License (version 2) as
+**  published by the Free Software Foundation.
+**
+**  This program is distributed in the hope that it will be useful,
+**  but WITHOUT ANY WARRANTY; without even the implied warranty of
+**  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+**  GNU General Public License for more details.
+**
+**  You should have received a copy of the GNU General Public License
+**  along with this program; if not, write to the Free Software
+**  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+******************************************************************************/
+
+#include "ct.h"
+#include "timer.h"
+
+
+enum { O_PHANTOM, O_DESC, O_NRAY, O_ROTANGLE, O_PHMFILE,
+       O_TRACE, O_VERBOSE, O_HELP, O_DEBUG, O_VERSION };
+
+static struct option phm2pj_options[] = 
+{
+  {"phantom", 1, 0, O_PHANTOM},
+  {"phmfile", 1, 0, O_PHMFILE},
+  {"desc", 1, 0, O_DESC},
+  {"nray", 1, 0, O_NRAY},
+  {"rotangle", 1, 0, O_ROTANGLE},
+  {"trace", 1, 0, O_TRACE},
+  {"verbose", 0, 0, O_VERBOSE},
+  {"help", 0, 0, O_HELP},
+  {"debug", 0, 0, O_DEBUG},
+  {"version", 0, 0, O_VERSION},
+  {0, 0, 0, 0}
+};
+
+
+void 
+phm2pj_usage (const char *program)
+{
+  cout << "usage: " << fileBasename(program) << " outfile ndet nview [--phantom phantom-name] [--phmfile filename] [OPTIONS]" << endl;
+  cout << "Calculate (projections) through phantom object, either a predefined --phantom or a --phmfile" << endl;
+  cout << "" << endl;
+  cout << "     outfile      Name of output file for raysums" << endl;
+  cout << "     ndet         Number of detectors" << endl;
+  cout << "     nview        Number of rotated views" << endl;
+  cout << "     --phantom    Phantom to use for projection" << endl;
+  cout << "        herman    Herman head phantom" << endl;
+  cout << "        rowland   Rowland head phantom" << endl;
+  cout << "        browland  Bordered Rowland head phantom" << endl;
+  cout << "        unitpulse Unit pulse phantom" << endl;
+  cout << "     --phmfile    Get Phantom from phantom file" << endl;
+  cout << "     --desc       Description of raysum" << endl;
+  cout << "     --nray       Number of rays per detector (default = 1)" << endl;
+  cout << "     --rotangle   Degrees to rotate view through, multiple of PI (default = 1)" << endl;
+  cout << "     --trace      Trace level to use" << endl;
+  cout << "        none      No tracing (default)" << endl;
+  cout << "        text      Trace text level" << endl;
+  cout << "        phm       Trace phantom image" << endl;
+  cout << "        rays      Trace rays" << endl;
+  cout << "        plot      Trace plot" << endl;
+  cout << "        clipping  Trace clipping" << endl;
+  cout << "     --verbose    Verbose mode" << endl;
+  cout << "     --debug      Debug mode" << endl;
+  cout << "     --version    Print version" << endl;
+  cout << "     --help       Print this help message" << endl;
+}
+
+#ifdef HAVE_MPI
+void GatherProjectionsMPI (MPIWorld& mpiWorld, Projections& pjGlobal, Projections& pjLocal, const int opt_debug);
+#endif
+
+int 
+phm2pj_main (int argc, char* argv[])
+{
+  Phantom phm;
+  ScannerGeometry opt_geometry = DETECTOR_PARALLEL;
+  char *opt_outfile = NULL;
+  string opt_desc;
+  string opt_phmfilename;
+  int opt_ndet, opt_nview;
+  int opt_nray = 1;
+  int opt_trace = 0;
+  int opt_phmnum = -1;
+  int opt_verbose = 0;
+  int opt_debug = 0;
+  double opt_rotangle = 1;
+  char* endptr = NULL;
+  char* endstr;
+
+#ifdef HAVE_MPI
+  MPIWorld mpiWorld (argc, argv);
+#endif
+
+  Timer timerProgram;
+
+#ifdef HAVE_MPI
+  if (mpiWorld.getRank() == 0) {
+#endif
+    while (1) {
+      int c = getopt_long(argc, argv, "", phm2pj_options, NULL);
+
+      if (c == -1)
+	break;
+      
+      switch (c) {
+      case O_PHANTOM:
+	if ((opt_phmnum = opt_set_phantom (optarg)) < 0) {
+	  phm2pj_usage(argv[0]);
+	  return (1);
+	}
+	phm.create (opt_phmnum);
+	break;
+      case O_PHMFILE:
+#ifdef HAVE_MPI
+	if (mpiWorld.getRank() == 0)
+	  cerr << "Can not read phantom from file in MPI mode" << endl;
+	return (1);
+#endif
+	opt_phmfilename = optarg;
+	phm.createFromFile (opt_phmfilename.c_str());
+	break;
+      case O_VERBOSE:
+	opt_verbose = 1;
+	break;
+      case O_DEBUG:
+	opt_debug = 1;
+	break;
+	break;
+      case O_TRACE:
+	if ((opt_trace = opt_set_trace(optarg)) < 0) {
+	  phm2pj_usage(argv[0]);
+	  return (1);
+	}
+	break;
+      case O_DESC:
+	opt_desc = optarg;
+	break;
+      case O_ROTANGLE:
+	opt_rotangle = strtod(optarg, &endptr);
+	endstr = optarg + strlen(optarg);
+	if (endptr != endstr) {
+	  cerr << "Error setting --rotangle to " << optarg << endl;
+	  phm2pj_usage(argv[0]);
+	  return (1);
+	}
+	break;
+      case O_NRAY:
+	opt_nray = strtol(optarg, &endptr, 10);
+	endstr = optarg + strlen(optarg);
+	if (endptr != endstr) {
+	  cerr << "Error setting --nray to %s" << optarg << endl;
+	  phm2pj_usage(argv[0]);
+	  return (1);
+	}
+	break;
+        case O_VERSION:
+#ifdef VERSION
+	cout << "Version: " << VERSION << endl;
+#else
+        cout << "Unknown version number" << endl;
+#endif
+	  exit(0);
+      case O_HELP:
+      case '?':
+	phm2pj_usage(argv[0]);
+	return (0);
+      default:
+	phm2pj_usage(argv[0]);
+	return (1);
+      }
+    }
+  
+    if (phm.nPElem() == 0) {
+      cerr << "No phantom defined" << endl;
+      phm2pj_usage(argv[0]);
+      return (1);
+    }
+    if (optind + 3 != argc) {
+      phm2pj_usage(argv[0]);
+      return (1);
+    }
+
+    opt_outfile = argv[optind];
+    opt_ndet = strtol(argv[optind+1], &endptr, 10);
+    endstr = argv[optind+1] + strlen(argv[optind+1]);
+    if (endptr != endstr) {
+      cerr << "Error setting --ndet to " << argv[optind+1] << endl;
+      phm2pj_usage(argv[0]);
+      return (1);
+    }
+    opt_nview = strtol(argv[optind+2], &endptr, 10);
+    endstr = argv[optind+2] + strlen(argv[optind+2]);
+    if (endptr != endstr) {
+      cerr << "Error setting --nview to " << argv[optind+2] << endl;
+      phm2pj_usage(argv[0]);
+      return (1);
+    }
+
+    ostringstream desc;
+    desc << "Raysum_Collect: NDet=" << opt_ndet << ", Nview=" << opt_nview << ", NRay=" << opt_nray << ", RotAngle=" << opt_rotangle << ", ";
+    if (opt_phmfilename.length()) {
+      desc << "PhantomFile=" << opt_phmfilename;
+    } else if (opt_phmnum != -1) {
+      desc << "Phantom=" << name_of_phantom(opt_phmnum);
+    }
+    if (opt_desc.length()) {
+      desc << ": " << opt_desc;
+    }
+    opt_desc = desc.str();
+#ifdef HAVE_MPI
+  }
+#endif
+
+#ifdef HAVE_MPI
+  mpiWorld.getComm().Bcast (&opt_rotangle, 1, MPI::DOUBLE, 0);
+  mpiWorld.getComm().Bcast (&opt_nview, 1, MPI::INT, 0);
+  mpiWorld.getComm().Bcast (&opt_ndet, 1, MPI::INT, 0);
+  mpiWorld.getComm().Bcast (&opt_nray, 1, MPI::INT, 0);
+  mpiWorld.getComm().Bcast (&opt_phmnum, 1, MPI::INT, 0);
+  mpiWorld.getComm().Bcast (&opt_verbose, 1, MPI::INT, 0);
+  mpiWorld.getComm().Bcast (&opt_debug, 1, MPI::INT, 0);
+  mpiWorld.getComm().Bcast (&opt_trace, 1, MPI::INT, 0);
+
+  if (mpiWorld.getRank() > 0 && opt_phmnum >= 0)
+    phm.create (opt_phmnum);
+#endif
+
+  opt_rotangle *= PI;
+  Scanner scanner (phm, opt_geometry, opt_ndet, opt_nview, opt_nray, opt_rotangle);
+
+#ifdef HAVE_MPI
+  mpiWorld.setTotalWorkUnits (opt_nview);
+
+  Projections pjGlobal;
+  if (mpiWorld.getRank() == 0) {
+    pjGlobal = Projections (scanner);
+  }
+  
+  Scanner localScanner (phm, opt_geometry, opt_ndet, mpiWorld.getMyLocalWorkUnits(), opt_nray, opt_rotangle);
+  Projections pjLocal (localScanner);
+  if (opt_debug)
+    cout << "pjLocal->nview = " << pjLocal.nView() << " (process " << mpiWorld.getRank() << ")" << endl;;
+
+  TimerMPI timerProject (mpiWorld.getComm());
+  localScanner.collectProjections (pjLocal, phm, mpiWorld.getMyStartWorkUnit(), opt_trace);
+  if (opt_verbose)
+    timerProject.timerEndAndReport ("Time to collect projections");
+
+  TimerMPI timerGather (mpiWorld.getComm());
+  GatherProjectionsMPI (mpiWorld, pjGlobal, pjLocal, opt_debug);
+  if (opt_verbose)
+    timerGather.timerEndAndReport ("Time to gather projections");
+#else
+  Projections pjGlobal (scanner);
+  scanner.collectProjections (pjGlobal, phm, 0, opt_trace);
+#endif
+  
+#ifdef HAVE_MPI
+  if (mpiWorld.getRank() == 0) 
+#endif
+    {
+      pjGlobal.setCalcTime (timerProgram.timerEnd());
+      pjGlobal.setRemark (opt_desc);
+      pjGlobal.write (opt_outfile);
+      if (opt_verbose) {
+	phm.print();
+	cout << endl;
+	pjGlobal.printScanInfo();
+	cout << endl;
+	cout << "Remark: " << pjGlobal.remark() << endl;
+	cout << "Run time: " << pjGlobal.calcTime() << " seconds" << endl;
+      }
+    }
+
+  return (0);
+}
+
+
+/* FUNCTION
+ *    GatherProjectionsMPI
+ *
+ * SYNOPSIS
+ *    Gather's raysums from all processes in pjLocal in pjGlobal in process 0
+ */
+
+#ifdef HAVE_MPI
+void GatherProjectionsMPI (MPIWorld& mpiWorld, Projections& pjGlobal, Projections& pjLocal, const int opt_debug)
+{
+  for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++) {
+    DetectorArray& detArray = pjLocal.getDetectorArray(iw);
+    double viewAngle = detArray.viewAngle();
+    int nDet = detArray.nDet();
+    DetectorValue* detval = detArray.detValues();
+
+    mpiWorld.getComm().Send(&viewAngle, 1, MPI::DOUBLE, 0, 0);
+    mpiWorld.getComm().Send(&nDet, 1, MPI::INT, 0, 0);
+    mpiWorld.getComm().Send(detval, nDet, MPI::FLOAT, 0, 0);
+  }
+
+  if (mpiWorld.getRank() == 0) {
+    for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) {
+      for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) {
+	MPI::Status status;
+	double viewAngle;
+	int nDet;
+	DetectorArray detArray = pjGlobal.getDetectorArray(iw);
+	DetectorValue* detval = detArray.detValues();
+
+	mpiWorld.getComm().Recv(&viewAngle, 1, MPI::DOUBLE, iProc, 0, status);
+	mpiWorld.getComm().Recv(&nDet, 1, MPI::INT, iProc, 0, status);
+	mpiWorld.getComm().Recv(detval, nDet, MPI::FLOAT, iProc, 0, status);
+	detArray.setViewAngle (viewAngle);
+      }
+    }
+  }
+
+}
+#endif
+
+
+#ifndef NO_MAIN
+int 
+main (int argc, char* argv[])
+{
+  return (phm2pj_main(argc, argv));
+}
+#endif
+