-#SUBDIRS = ctsim
-bin_PROGRAMS = pjrec phm2pj pj2if @lamprograms@ ifinfo phm2if if-1 if-2 if2img
-bin_SCRIPTS = sample-ctsim.sh
-EXTRA_PROGRAMS = pjrec-lam phm2pj-lam phm2if-lam
-INCLUDES=@my_includes@
-EXTRA_DIST=Makefile.nt mpiworld.cpp
-SOURCE_DEPEND=../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a ../libctgraphics/libctgraphics.a
+bin_PROGRAMS=ctsim
-pjrec_SOURCES = pjrec.cpp
-pjrec_LDADD=@ctlibs@
-pjrec_DEPENDENCIES=$(SOURCE_DEPEND)
-phm2pj_SOURCES=phm2pj.cpp
-phm2pj_LDADD=@ctlibs@
-phm2pj_DEPENDENCIES=$(SOURCE_DEPEND)
-phm2if_SOURCES = phm2if.cpp
-phm2if_LDADD=@ctlibs@
-phm2if_DEPENDENCIES=$(SOURCE_DEPEND)
-if2img_SOURCES = if2img.cpp
-if2img_LDADD=@ctlibs@
-if2img_DEPENDENCIES=$(SOURCE_DEPEND)
-pj2if_SOURCES = pj2if.cpp
-pj2if_LDADD=@ctlibs@
-pj2if_DEPENDENCIES=$(SOURCE_DEPEND)
-if_1_SOURCES=if-1.cpp
-if_1_LDADD=@ctlibs@
-if_1_DEPENDENCIES=$(SOURCE_DEPEND)
-if_2_SOURCES=if-2.cpp
-if_2_LDADD=@ctlibs@
-if_2_DEPENDENCIES=$(SOURCE_DEPEND)
-ifinfo_SOURCES = ifinfo.cpp
-ifinfo_LDADD=@ctlibs@
-ifinfo_DEPENDENCIES=$(SOURCE_DEPEND)
-
-pjrec_lam_SOURCES=pjrec.cpp
-pjrec_lam_LDADD=@ctlamlibs@
-phm2if_lam_SOURCES=phm2if.cpp
-phm2if_lam_LDADD=@ctlamlibs@
-phm2pj_lam_SOURCES=phm2pj.cpp
-phm2pj_lam_LDADD=@ctlamlibs@
-
-realclean:
- rm -f *.pgm *.if *~ *.pj
-
-if USE_LAM
-CC_LAM = $(lamdir)/bin/balky
-LAM_EXTRA_SRC = mpiworld.cpp
-
-pjrec-lam: pjrec.cpp mpiworld.cpp ../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a
- $(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI pjrec.cpp -o pjrec-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@
-
-phm2pj-lam: phm2pj.cpp mpiworld.cpp ../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a
- $(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI phm2pj.cpp -o phm2pj-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@
-
-phm2if-lam: phm2if.cpp mpiworld.cpp ../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a
- $(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI phm2if.cpp -o phm2if-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@
-
-endif
-
-shared: pjrec.cpp phm2pj.cpp
- $(CC) @DEFS@ $(CFLAGS) $(INCLUDES) -DNO_MAIN -shared phm2pj.cpp pjrec.cpp @ctlibs@ $(LDFLAGS) -o ir.so
+ctsim_SOURCES=ctsim.cpp docs.cpp views.cpp
+ctsim_DEPENDENCIES=../libctgraphics/libctgraphics.a ../libctsupport/libctsupport.a ../libctsim/libctsim.a ../include/ct.h
+ctsim_LDADD=-L../libctgraphics -L../libctsupport -L../libctsim @wxlibs@
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: ctsim.h
+** Purpose: Header file for top-level routines of CTSim program
+** Programmer: Kevin Rosenberg
+** Date Started: July 2000
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: ctsim.cpp,v 1.4 2000/07/13 07:01:59 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
+******************************************************************************/
+
+// For compilers that support precompilation, includes "wx/wx.h".
+#include "wx/wxprec.h"
+
+#ifdef __BORLANDC__
+#pragma hdrstop
+#endif
+
+#ifndef WX_PRECOMP
+#include "wx/wx.h"
+#endif
+
+#if !wxUSE_DOC_VIEW_ARCHITECTURE
+#error You must set wxUSE_DOC_VIEW_ARCHITECTURE to 1 in setup.h!
+#endif
+
+#include "wx/docview.h"
+
+#include "ctsim.h"
+#include "docs.h"
+#include "views.h"
+
+IMPLEMENT_APP(CTSimApp)
+
+CTSimApp::CTSimApp(void)
+ : m_docManager(NULL), m_pFrame(NULL)
+{
+}
+
+bool
+CTSimApp::OnInit(void)
+{
+ m_docManager = new wxDocManager;
+
+ (void) new wxDocTemplate (m_docManager, "ImageFile", "*.if", "", "if", "ImageFile doc", "ImageFile View", CLASSINFO(ImageFileDocument), CLASSINFO(ImageFileView));
+
+ (void) new wxDocTemplate (m_docManager, "ProjectionFile", "*.pj", "", "pj", "ProjectionFile doc", "ProjectionFile View", CLASSINFO(ProjectionFileDocument), CLASSINFO(ProjectionFileView));
+
+ //// Create the main frame window
+ m_pFrame = new MainFrame(m_docManager, (wxFrame *) NULL, -1, "CTSim", wxPoint(0, 0), wxSize(500, 400), wxDEFAULT_FRAME_STYLE);
+
+ //// Make a menubar
+ wxMenu *file_menu = new wxMenu;
+
+ // file_menu->Append(wxID_NEW, "&New...");
+ file_menu->Append(wxID_OPEN, "&Open...");
+
+ file_menu->AppendSeparator();
+ file_menu->Append(MAINMENU_FILE_EXIT, "E&xit");
+
+ // A nice touch: a history of files visited. Use this menu.
+ m_docManager->FileHistoryUseMenu(file_menu);
+
+ wxMenu *help_menu = new wxMenu;
+ help_menu->Append(MAINMENU_HELP_ABOUT, "&About");
+
+ wxMenuBar *menu_bar = new wxMenuBar;
+
+ menu_bar->Append(file_menu, "&File");
+ menu_bar->Append(help_menu, "&Help");
+
+ m_pFrame->SetMenuBar(menu_bar);
+
+ m_pFrame->Centre(wxBOTH);
+ m_pFrame->Show(true);
+
+ SetTopWindow (m_pFrame);
+ return true;
+}
+
+int
+CTSimApp::OnExit(void)
+{
+ delete m_docManager;
+ return 0;
+}
+
+
+// Top-level window for CTSim
+
+IMPLEMENT_CLASS(MainFrame, wxDocParentFrame)
+BEGIN_EVENT_TABLE(MainFrame, wxDocParentFrame)
+ EVT_MENU(MAINMENU_HELP_ABOUT, MainFrame::OnAbout)
+ EVT_MENU(MAINMENU_FILE_EXIT, MainFrame::OnExit)
+END_EVENT_TABLE()
+
+MainFrame::MainFrame(wxDocManager *manager, wxFrame *frame, wxWindowID id, const wxString& title, const wxPoint& pos, const wxSize& size, const long type)
+ : wxDocParentFrame(manager, frame, id, title, pos, size, type)
+{
+ CreateStatusBar();
+ SetStatusText ("Welcome to CTSim");
+}
+
+void
+MainFrame::OnAbout(wxCommandEvent& WXUNUSED(event) )
+{
+ wxMessageBox("CTSim\nAuthor: Kevin Rosenberg <kevin@rosenberg.net>\nUsage: ctsim", "About CTSim", wxOK | wxICON_INFORMATION, this);
+}
+
+void
+MainFrame::OnExit (wxCommandEvent& WXUNUSED(event) )
+{
+ Close(true);
+}
+
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: ctsim.cpp
+** Purpose: Top-level routines for CTSim program
+** Programmer: Kevin Rosenberg
+** Date Started: July 2000
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: ctsim.h,v 1.1 2000/07/13 07:01:59 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
+******************************************************************************/
+
+#ifndef __CTSIMH__
+#define __CTSIMH__
+
+#include "wx/docview.h"
+
+class wxDocManager;
+class MainFrame;
+class CTSimApp: public wxApp
+{
+public:
+ CTSimApp(void);
+ bool OnInit(void);
+ int OnExit(void);
+ MainFrame* getMainFrame(void) const
+ { return m_pFrame; }
+
+protected:
+ wxDocManager* m_docManager;
+ MainFrame* m_pFrame;
+};
+
+DECLARE_APP(CTSimApp)
+
+// Define a new frame
+class MainFrame: public wxDocParentFrame
+{
+ DECLARE_CLASS(MainFrame)
+
+public:
+ MainFrame(wxDocManager *manager, wxFrame *frame, wxWindowID id, const wxString& title, const wxPoint& pos, const wxSize& size, const long type);
+
+ void OnAbout (wxCommandEvent& event);
+ void OnExit (wxCommandEvent& event);
+
+ DECLARE_EVENT_TABLE()
+};
+
+extern MainFrame *GetMainFrame(void);
+
+
+enum {
+ MAINMENU_HELP_ABOUT = 500,
+ MAINMENU_FILE_EXIT,
+ IFMENU_FILE_PROPERTIES,
+ PJMENU_FILE_PROPERTIES,
+ PJMENU_FILE_RECONSTRUCT,
+};
+
+#endif
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: doc.cpp
+** Purpose: Document routines for CTSim program
+** Programmer: Kevin Rosenberg
+** Date Started: July 2000
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: docs.cpp,v 1.1 2000/07/13 07:01:59 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
+******************************************************************************/
+
+#ifdef __GNUG__
+// #pragma implementation
+#endif
+
+// For compilers that support precompilation, includes "wx/wx.h".
+#include "wx/wxprec.h"
+
+#ifdef __BORLANDC__
+#pragma hdrstop
+#endif
+
+#ifndef WX_PRECOMP
+#include "wx/wx.h"
+#endif
+#include "wx/txtstrm.h"
+
+#if !wxUSE_DOC_VIEW_ARCHITECTURE
+#error You must set wxUSE_DOC_VIEW_ARCHITECTURE to 1 in setup.h!
+#endif
+
+#include "docs.h"
+#include "views.h"
+
+
+
+// ImageFileDocument
+
+IMPLEMENT_DYNAMIC_CLASS(ImageFileDocument, wxDocument)
+
+bool ImageFileDocument::OnSaveDocument(const wxString& filename)
+{
+ if (!m_imageFile.fileWrite (filename))
+ return false;
+ Modify(false);
+ return true;
+}
+
+bool ImageFileDocument::OnOpenDocument(const wxString& filename)
+{
+ if (! m_imageFile.fileRead (filename))
+ return false;
+
+ SetFilename(filename, true);
+ Modify(false);
+ UpdateAllViews();
+ return true;
+}
+
+bool ImageFileDocument::IsModified(void) const
+{
+ return wxDocument::IsModified();
+}
+
+void ImageFileDocument::Modify(bool mod)
+{
+ wxDocument::Modify(mod);
+}
+
+// ProjectionFileDocument
+
+IMPLEMENT_DYNAMIC_CLASS(ProjectionFileDocument, wxDocument)
+
+bool ProjectionFileDocument::OnSaveDocument(const wxString& filename)
+{
+ if (!m_projectionFile.write (filename))
+ return false;
+ Modify(false);
+ return true;
+}
+
+bool ProjectionFileDocument::OnOpenDocument(const wxString& filename)
+{
+ if (! m_projectionFile.read (filename))
+ return false;
+
+ SetFilename(filename, true);
+ Modify(false);
+ UpdateAllViews();
+ return true;
+}
+
+bool ProjectionFileDocument::IsModified(void) const
+{
+ return wxDocument::IsModified();
+}
+
+void ProjectionFileDocument::Modify(bool mod)
+{
+ wxDocument::Modify(mod);
+}
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: doc.h
+** Purpose: Header file for Document routines of CTSim program
+** Programmer: Kevin Rosenberg
+** Date Started: July 2000
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: docs.h,v 1.1 2000/07/13 07:01:59 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
+******************************************************************************/
+
+#ifdef __GNUG__
+// #pragma interface
+#endif
+
+#ifndef __DOCSH__
+#define __DOCSH__
+
+#include "wx/docview.h"
+#include "imagefile.h"
+#include "projections.h"
+
+
+class ImageFileDocument: public wxDocument
+{
+ DECLARE_DYNAMIC_CLASS(ImageFileDocument)
+private:
+ ImageFile m_imageFile;
+
+public:
+ virtual bool OnSaveDocument(const wxString& filename);
+ virtual bool OnOpenDocument(const wxString& filename);
+ virtual bool IsModified(void) const;
+ virtual void Modify(bool mod);
+
+ ImageFileDocument(void) {}
+ ~ImageFileDocument(void) {}
+
+ const ImageFile& getImageFile(void) const
+ { return m_imageFile; }
+
+ ImageFile& getImageFile(void)
+ { return m_imageFile; }
+};
+
+
+class ProjectionFileDocument: public wxDocument
+{
+ DECLARE_DYNAMIC_CLASS(ProjectionFileDocument)
+private:
+ Projections m_projectionFile;
+
+public:
+ virtual bool OnSaveDocument(const wxString& filename);
+ virtual bool OnOpenDocument(const wxString& filename);
+ virtual bool IsModified(void) const;
+ virtual void Modify(bool mod);
+
+ ProjectionFileDocument(void) {}
+ ~ProjectionFileDocument(void) {}
+
+ const Projections& getProjections(void) const
+ { return m_projectionFile; }
+
+ Projections& getProjections(void)
+ { return m_projectionFile; }
+};
+
+class TextEditDocument: public wxDocument
+{
+ DECLARE_DYNAMIC_CLASS(TextEditDocument)
+private:
+public:
+ virtual bool OnSaveDocument(const wxString& filename);
+ virtual bool OnOpenDocument(const wxString& filename);
+ virtual bool IsModified(void) const;
+ virtual void Modify(bool mod);
+
+ TextEditDocument(void) {}
+ ~TextEditDocument(void) {}
+};
+
+
+#endif
+++ /dev/null
-/*****************************************************************************
-** FILE IDENTIFICATION
-**
-** Name: if-1.cpp
-** Purpose: Manipulate a single image file
-** Programmer: Kevin Rosenberg
-** Date Started: Aug 1984
-**
-** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
-**
-** $Id: if-1.cpp,v 1.9 2000/06/26 21:15:24 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
-******************************************************************************/
-
-/* FILE
- * if-1.c Filter a single IF file
- */
-
-#include "ct.h"
-
-enum {O_LOG, O_EXP, O_SQRT, O_SQR, O_INVERT, O_VERBOSE, O_HELP, O_VERSION};
-
-static struct option my_options[] =
-{
- {"invert", 0, 0, O_INVERT},
- {"verbose", 0, 0, O_VERBOSE},
- {"log", 0, 0, O_LOG},
- {"exp", 0, 0, O_EXP},
- {"sqr", 0, 0, O_SQR},
- {"sqrt", 0, 0, O_SQRT},
- {"help", 0, 0, O_HELP},
- {"version", 0, 0, O_VERSION},
- {0, 0, 0, 0}
-};
-
-void
-if1_usage (const char *program)
-{
- cout << "usage: " << fileBasename(program) << " infile outfile [OPTIONS]" << endl;
- cout << "Generate a IF file from a IF file" << endl;
- cout << endl;
- cout << " --invert Invert image" << endl;
- cout << " --log Natural logrithm of image" << endl;
- cout << " --exp Natural exponential of image" << endl;
- cout << " --sqr Square of image" << endl;
- cout << " --sqrt Square root of image" << endl;
- cout << " --verbose Verbose modem" << endl;
- cout << " --version Print version" << endl;
- cout << " --help Print this help message" << endl;
-}
-
-int
-if1_main (int argc, char *const argv[])
-{
- ImageFile *im_in;
- ImageFile *im_out;
- char *in_file;
- char *out_file;
- int opt_verbose = 0;
- int opt_invert = 0;
- int opt_log = 0;
- int opt_exp = 0;
- int opt_sqr = 0;
- int opt_sqrt = 0;
-
- while (1)
- {
- int c = getopt_long (argc, argv, "", my_options, NULL);
-
- if (c == -1)
- break;
-
- switch (c)
- {
- case O_INVERT:
- opt_invert = 1;
- break;
- case O_LOG:
- opt_log = 1;
- break;
- case O_SQR:
- opt_sqr = 1;
- break;
- case O_SQRT:
- opt_sqrt = 1;
- break;
- case O_EXP:
- opt_exp = 1;
- break;
- case O_VERBOSE:
- opt_verbose = 1;
- break;
- case O_VERSION:
-#ifdef VERSION
- cout << "Version " << VERSION << endl;
-#else
- cout << "Unknown version number" << endl;
-#endif
- return (0);
- case O_HELP:
- case '?':
- if1_usage(argv[0]);
- return (0);
- default:
- if1_usage(argv[0]);
- return (1);
- }
- }
-
- if (optind + 2 != argc)
- {
- if1_usage(argv[0]);
- return (1);
- }
-
- in_file = argv[optind];
- out_file = argv[optind + 1];
-
-
- string histString;
-
- if (opt_invert || opt_log || opt_exp || opt_sqr || opt_sqrt) {
- int ix, iy;
-
- im_in = new ImageFile ();
- im_in->fileRead (in_file);
- int nx = im_in->nx();
- int ny = im_in->ny();
- im_out = new ImageFile (nx, ny);
-
- ImageFileArray vIn = im_in->getArray();
- ImageFileArray vOut = im_out->getArray();
-
- if (opt_invert) {
- for (ix = 0; ix < nx; ix++)
- for (iy = 0; iy < ny; iy++)
- vOut[ix][iy] = - vIn[ix][iy];
- histString = "Invert transformation";
- }
- if (opt_log) {
- for (ix = 0; ix < nx; ix++)
- for (iy = 0; iy < ny; iy++)
- vOut[ix][iy] = log (vIn[ix][iy]);
- histString = "Log transformation";
- }
- if (opt_exp) {
- for (ix = 0; ix < nx; ix++)
- for (iy = 0; iy < ny; iy++)
- vOut[ix][iy] = exp (vIn[ix][iy]);
- histString = "Exp transformation";
- }
- if (opt_sqr) {
- for (ix = 0; ix < nx; ix++)
- for (iy = 0; iy < ny; iy++)
- vOut[ix][iy] = vIn[ix][iy] * vIn[ix][iy];
- histString = "Sqr transformation";
- }
- if (opt_sqrt) {
- for (ix = 0; ix < nx; ix++)
- for (iy = 0; iy < ny; iy++)
- vOut[ix][iy] = sqrt (vIn[ix][iy]);
- histString = "Sqrt transformation";
- }
-
- im_out->labelsCopy (*im_in);
- im_out->labelAdd (Array2dFileLabel::L_HISTORY, histString.c_str());
- im_out->fileWrite (out_file);
- }
-
- return (0);
-}
-
-#ifndef NO_MAIN
-int
-main (int argc, char *const argv[])
-{
- int retval = 1;
-
- try {
- retval = if1_main(argc, argv);
- } catch (exception e) {
- cerr << "Exception: " << e.what() << endl;
- } catch (...) {
- cerr << "Unknown exception" << endl;
- }
-
- return (retval);
-}
-#endif
-
+++ /dev/null
-/*****************************************************************************
-** FILE IDENTIFICATION
-**
-** Name: if-2.cpp
-** Purpose: Manipulate two image files
-** Programmer: Kevin Rosenberg
-** Date Started: May 2000
-**
-** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
-**
-** $Id: if-2.cpp,v 1.10 2000/07/09 08:16:18 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
-******************************************************************************/
-
-/* FILE
- * if-2.c Generate a IF file from two input IF files
- */
-
-#include "ct.h"
-#include "timer.h"
-
-enum {O_ADD, O_SUB, O_MUL, O_COMP, O_ROW_PLOT, O_COLUMN_PLOT, O_VERBOSE, O_HELP, O_VERSION};
-
-static struct option my_options[] =
-{
- {"add", 0, 0, O_ADD},
- {"sub", 0, 0, O_SUB},
- {"mul", 0, 0, O_MUL},
- {"comp", 0, 0, O_COMP},
- {"column-plot", 1, 0, O_COLUMN_PLOT},
- {"row-plot", 1, 0, O_ROW_PLOT},
- {"verbose", 0, 0, O_VERBOSE},
- {"help", 0, 0, O_HELP},
- {"version", 0, 0, O_VERSION},
- {0, 0, 0, 0}
-};
-
-void
-if2_usage (const char *program)
-{
- cout << "usage: " << fileBasename(program) << " infile1 infile2 outfile [OPTIONS]" << endl;
- cout << "Perform functions on two input image files" << endl;
- cout << endl;
- cout << " infile1 Name of first input IF file" << endl;
- cout << " infile2 Name of second input IF file" << endl;
- cout << " outfile Name of output IF file" << endl;
- cout << " --add Add images" << endl;
- cout << " --sub Subtract image 2 from image 1" << endl;
- cout << " --mul Multiply images" << endl;
- cout << " --comp Compare images" << endl;
- cout << " --column-plot n Plot column\n";
- cout << " --row-plot n Plot row\n";
- cout << " --verbose Verbose modem" << endl;
- cout << " --version Print version" << endl;
- cout << " --help Print this help message" << endl;
-}
-
-int
-if2_main (int argc, char *const argv[])
-{
- ImageFile* pim_in1;
- ImageFile* pim_in2;
- ImageFile* pim_out;
- string in_file1;
- string in_file2;
- string out_file;
- int opt_verbose = 0;
- int opt_add = 0;
- int opt_sub = 0;
- int opt_mul = 0;
- int opt_comp = 0;
- int opt_outputFile = 0;
- int opt_rowPlot = -1;
- int opt_columnPlot = -1;
- Timer timerProgram;
-
- while (1) {
- char* endptr;
- int c = getopt_long (argc, argv, "", my_options, NULL);
-
- if (c == -1)
- break;
-
- switch (c) {
- case O_ADD:
- opt_add = 1;
- opt_outputFile = 1;
- break;
- case O_SUB :
- opt_sub = 1;
- opt_outputFile = 1;
- break;
- case O_MUL:
- opt_mul = 1;
- opt_outputFile = 1;
- break;
- case O_ROW_PLOT:
- opt_rowPlot = strtol(optarg, &endptr, 10);
- if (endptr != optarg + strlen(optarg)) {
- if2_usage(argv[0]);
- }
- case O_COLUMN_PLOT:
- opt_columnPlot = strtol(optarg, &endptr, 10);
- if (endptr != optarg + strlen(optarg)) {
- if2_usage(argv[0]);
- }
- case O_COMP:
- opt_comp = 1;
- break;
- case O_VERBOSE:
- opt_verbose = 1;
- break;
- case O_VERSION:
-#ifdef VERSION
- cout << "Version " << VERSION << endl;
-#else
- cout << "Unknown version number" << endl;
-#endif
- return (0);
- case O_HELP:
- case '?':
- if2_usage(argv[0]);
- return (0);
- default:
- if2_usage(argv[0]);
- return (1);
- }
- }
-
- if (opt_outputFile && (optind + 3 != argc)) {
- if2_usage(argv[0]);
- return (1);
- }
- else if (! opt_outputFile && optind + 2 != argc) {
- if2_usage(argv[0]);
- return (1);
- }
-
- in_file1 = argv[optind];
- in_file2 = argv[optind + 1];
- if (opt_outputFile)
- out_file = argv[optind + 2];
-
- pim_in1 = new ImageFile ();
- pim_in2 = new ImageFile ();
- ImageFile& im_in1 = *pim_in1;
- ImageFile& im_in2 = *pim_in2;
-
- if (! im_in1.fileRead(in_file1) || ! im_in2.fileRead(in_file2)) {
- sys_error (ERR_WARNING, "Error reading an image");
- return (1);
- }
-
- if (im_in1.nx() != im_in2.nx() || im_in1.ny() != im_in2.ny()) {
- sys_error (ERR_SEVERE, "Error: Size of image 1 (%d,%d) and image 2 (%d,%d) do not match",
- im_in1.nx(), im_in1.ny(), im_in2.nx(), im_in2.ny());
- return(1);
- }
- if (im_in1.nx() < 0 || im_in1.ny() < 0) {
- sys_error (ERR_SEVERE, "Error: Size of image < 0");
- return(1);
- }
-
- ImageFileArray v1 = im_in1.getArray();
- ImageFileArray v2 = im_in2.getArray();
- ImageFileArray vout = NULL;
-
- if (opt_outputFile) {
- pim_out = new ImageFile (im_in1.nx(), im_in1.ny());
- vout = pim_out->getArray();
- }
-
- string strOperation;
- if (opt_add) {
- strOperation = "Add Images";
- for (int ix = 0; ix < im_in1.nx(); ix++) {
- ImageFileColumn in1 = v1[ix];
- ImageFileColumn in2 = v2[ix];
- ImageFileColumn out = vout[ix];
- for (int iy = 0; iy < im_in1.ny(); iy++)
- *out++ = *in1++ + *in2++;
- }
- } else if (opt_sub) {
- strOperation = "Subtract Images";
- for (int ix = 0; ix < im_in1.nx(); ix++) {
- ImageFileColumn in1 = v1[ix];
- ImageFileColumn in2 = v2[ix];
- ImageFileColumn out = vout[ix];
- for (int iy = 0; iy < im_in1.ny(); iy++)
- *out++ = *in1++ - *in2++;
- }
- } else if (opt_mul) {
- strOperation = "Multiply Images";
- for (int ix = 0; ix < im_in1.nx(); ix++) {
- ImageFileColumn in1 = v1[ix];
- ImageFileColumn in2 = v2[ix];
- ImageFileColumn out = vout[ix];
- for (int iy = 0; iy < im_in1.ny(); iy++)
- *out++ = *in1++ * *in2++;
- }
- }
- if (opt_comp) {
- double d, r, e;
- im_in1.comparativeStatistics (im_in2, d, r, e);
- cout << "d=" << d << ", r=" << r << ", e=" << e << endl;
- }
- if (opt_columnPlot > 0) {
- if (opt_columnPlot >= im_in1.nx() || opt_columnPlot >= im_in2.nx()) {
- sys_error (ERR_SEVERE, "column-plot > nx");
- return (1);
- }
- double plot_xaxis [im_in1.nx()];
- for (int i = 0; i < im_in1.nx(); i++)
- plot_xaxis[i] = i;
-#if HAVE_SGP
-#if 0
-#else
- ezset ("clear.");
- ezset ("xticks major 5.");
- ezset ("xlabel Column");
- ezset ("ylabel Pixel");
- ezset ("curves 2");
- ezset ("box.");
- ezset ("grid.");
- ezplot (v1[opt_columnPlot], plot_xaxis, im_in1.ny());
- ezplot (v2[opt_columnPlot], plot_xaxis, im_in2.ny());
-#endif
- char str[256];
- cout << "Press enter to continue" << flush;
- fgets(str, sizeof(str), stdin);
- sgp2_close (sgp2_get_active_win());
-#endif
- }
-
- if (opt_rowPlot > 0) {
- if (opt_rowPlot >= im_in1.ny() || opt_rowPlot >= im_in2.ny()) {
- sys_error (ERR_SEVERE, "row_plot > ny");
- return (1);
- }
- double plot_xaxis [im_in1.ny()];
- double v1Row[im_in1.nx()], v2Row[im_in2.nx()];
-
- for (int i = 0; i < im_in1.ny(); i++)
- plot_xaxis[i] = i;
- for (int i = 0; i < im_in1.nx(); i++)
- v1Row[i] = v1[opt_rowPlot][i];
- for (int i = 0; i < im_in2.nx(); i++)
- v2Row[i] = v2[opt_rowPlot][i];
-
-#if HAVE_SGP
-#if 0
-#else
- ezset ("clear.");
- ezset ("xticks major 5.");
- ezset ("xlabel Column");
- ezset ("ylabel Pixel");
- ezset ("curves 2");
- ezset ("box.");
- ezset ("grid.");
- ezplot (v1Row, plot_xaxis, im_in1.nx());
- ezplot (v2Row, plot_xaxis, im_in2.nx());
-#endif
- char str[256];
- cout << "Press enter to continue" << flush;
- fgets(str, sizeof(str), stdin);
- sgp2_close (sgp2_get_active_win());
-#endif
- }
-
- if (opt_outputFile) {
- pim_out->labelsCopy (im_in1, "if-2 file 1: ");
- pim_out->labelsCopy (im_in2, "if-2 file 2: ");
- pim_out->labelAdd (Array2dFileLabel::L_HISTORY, strOperation.c_str(), timerProgram.timerEnd());
- pim_out->fileWrite(out_file);
- }
-
- return (0);
-}
-
-#ifndef NO_MAIN
-int
-main (int argc, char *const argv[])
-{
- int retval = 1;
-
- try {
- retval = if2_main(argc, argv);
- } catch (exception e) {
- cerr << "Exception: " << e.what() << endl;
- } catch (...) {
- cerr << "Unknown exception" << endl;
- }
-
- return (retval);
-}
-#endif
-
+++ /dev/null
-/*****************************************************************************
-** FILE IDENTIFICATION
-**
-** Name: if2img.cpp
-** Purpose: Convert an image file to a viewable image
-** Programmer: Kevin Rosenberg
-** Date Started: April 2000
-**
-** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
-**
-** $Id: if2img.cpp,v 1.13 2000/07/11 10:32:44 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"
-
-
-enum { O_SCALE, O_MIN, O_MAX, O_AUTO, O_CENTER, O_STATS, O_FORMAT, O_LABELS,
- O_HELP, O_VERBOSE, O_VERSION, O_DEBUG };
-
-static struct option my_options[] =
-{
- {"scale", 1, 0, O_SCALE},
- {"min", 1, 0, O_MIN},
- {"max", 1, 0, O_MAX},
- {"auto", 1, 0, O_AUTO},
- {"center", 1, 0, O_CENTER},
- {"format", 1, 0, O_FORMAT},
- {"labels", 0, 0, O_LABELS},
- {"stats", 0, 0, O_STATS},
- {"verbose", 0, 0, O_VERBOSE},
- {"debug", 0, 0, O_DEBUG},
- {"help", 0, 0, O_HELP},
- {"version", 0, 0, O_VERSION},
- {0, 0, 0, 0}
-};
-
-enum { O_AUTO_FULL, O_AUTO_STD0_1, O_AUTO_STD0_5, O_AUTO_STD1, O_AUTO_STD2, O_AUTO_STD3 };
-static const char O_AUTO_FULL_STR[]="full";
-static const char O_AUTO_STD0_1_STR[]="std0.1";
-static const char O_AUTO_STD0_5_STR[]="std0.5";
-static const char O_AUTO_STD1_STR[]="std1";
-static const char O_AUTO_STD2_STR[]="std2";
-static const char O_AUTO_STD3_STR[]="std3";
-
-enum { O_CENTER_MEAN, O_CENTER_MODE };
-static const char O_CENTER_MEAN_STR[]="mean";
-static const char O_CENTER_MODE_STR[]="mode";
-
-enum { O_FORMAT_GIF, O_FORMAT_PNG, O_FORMAT_PNG16, O_FORMAT_PGM, O_FORMAT_PGMASC, O_FORMAT_DISP };
-static const char O_FORMAT_GIF_STR[]="gif";
-static const char O_FORMAT_PNG_STR[]="png" ;
-static const char O_FORMAT_PNG16_STR[]="png16" ;
-static const char O_FORMAT_PGM_STR[]="pgm";
-static const char O_FORMAT_PGMASC_STR[]="pgmasc";
-static const char O_FORMAT_DISP_STR[]="disp";
-
-void
-if2img_usage (const char *program)
-{
- cout << "usage: " << fileBasename(program) << " ifname outfile [OPTIONS]" << endl;
- cout << "Convert IF file to an image file" << endl;
- cout << endl;
- cout << " ifname Name of input file" << endl;
- cout << " outfile Name of output file" << endl;
- cout << " --format Output format" << endl;
- cout << " pgm PGM (portable graymap) format (default)" << endl;
- cout << " pgmasc PGM (portable graymap) ASCII format" << endl;
-#ifdef HAVE_PNG
- cout << " png PNG (8-bit) format" << endl;
- cout << " png16 PNG (16-bit) format" << endl;
-#endif
-#if HAVE_G2
- cout << " gif GIF format" << endl;
-#endif
- cout << " disp Display on screen" << endl;
- cout << " --center Center of window" << endl;
- cout << " mode Mode is center of window (default)" << endl;
- cout << " mean Mean is center of window" << endl;
- cout << " --auto Set auto window" << endl;
- cout << " full Use full window (default)" << endl;
- cout << " std0.1 Use 0.1 standard deviation about center" << endl;
- cout << " std0.5 Use 0.5 standard deviation about center" << endl;
- cout << " std1 Use one standard deviation about center" << endl;
- cout << " std2 Use two standard deviations about center" << endl;
- cout << " std3 Use three standard deviations about center" << endl;
- cout << " --scale Scaling factor for output size" << endl;
- cout << " --min Set minimum intensity" << endl;
- cout << " --max Set maximum intensity" << endl;
- cout << " --stats Print image statistics" << endl;
- cout << " --labels Print image labels" << endl;
- cout << " --debug Set debug mode" << endl;
- cout << " --verbose Set verbose mode" << endl;
- cout << " --version Print version" << endl;
- cout << " --help Print this help message" << endl;
-}
-
-
-int
-if2img_main (int argc, char *const argv[])
-{
- ImageFile* pim = NULL;
- double densmin = HUGE_VAL, densmax = -HUGE_VAL;
- char *in_file, *out_file;
- int opt_verbose = 0;
- int opt_scale = 1;
- int opt_auto = O_AUTO_FULL;
- int opt_set_max = 0;
- int opt_set_min = 0;
- int opt_stats = 0;
- int opt_debug = 0;
- int opt_center = O_CENTER_MODE;
- int opt_format = O_FORMAT_PGM;
- int opt_labels = 0;
-
- while (1)
- {
- int c = getopt_long (argc, argv, "", my_options, NULL);
- char *endptr = NULL;
- char *endstr;
-
- if (c == -1)
- break;
-
- switch (c)
- {
- case O_MIN:
- opt_set_min = 1;
- densmin = strtod(optarg, &endptr);
- endstr = optarg + strlen(optarg);
- if (endptr != endstr)
- {
- sys_error (ERR_SEVERE, "Error setting --min to %s", optarg);
- if2img_usage(argv[0]);
- return (1);
- }
- break;
- case O_MAX:
- opt_set_max = 1;
- densmax = strtod(optarg, &endptr);
- endstr = optarg + strlen(optarg);
- if (endptr != endstr)
- {
- sys_error (ERR_SEVERE, "Error setting --max to %s", optarg);
- if2img_usage(argv[0]);
- return (1);
- }
- break;
- case O_SCALE:
- opt_scale = strtol(optarg, &endptr, 10);
- endstr = optarg + strlen(optarg);
- if (endptr != endstr)
- {
- sys_error (ERR_SEVERE, "Error setting --scale to %s", optarg);
- if2img_usage(argv[0]);
- return (1);
- }
- break;
- case O_AUTO:
- if (strcmp(optarg, O_AUTO_FULL_STR) == 0)
- opt_auto = O_AUTO_FULL;
- else if (strcmp(optarg, O_AUTO_STD1_STR) == 0)
- opt_auto = O_AUTO_STD1;
- else if (strcmp(optarg, O_AUTO_STD0_5_STR) == 0)
- opt_auto = O_AUTO_STD0_5;
- else if (strcmp(optarg, O_AUTO_STD0_1_STR) == 0)
- opt_auto = O_AUTO_STD0_1;
- else if (strcmp(optarg, O_AUTO_STD2_STR) == 0)
- opt_auto = O_AUTO_STD2;
- else if (strcmp(optarg, O_AUTO_STD3_STR) == 0)
- opt_auto = O_AUTO_STD3;
- else
- {
- sys_error (ERR_SEVERE, "Invalid auto mode %s", optarg);
- if2img_usage(argv[0]);
- return (1);
- }
- break;
- case O_CENTER:
- if (strcmp(optarg, O_CENTER_MEAN_STR) == 0)
- opt_center = O_CENTER_MEAN;
- else if (strcmp(optarg, O_CENTER_MODE_STR) == 0)
- opt_center = O_CENTER_MODE;
- else
- {
- sys_error (ERR_SEVERE, "Invalid center mode %s", optarg);
- if2img_usage(argv[0]);
- return (1);
- }
- break;
- case O_FORMAT:
- if (strcmp(optarg, O_FORMAT_PGM_STR) == 0)
- opt_format = O_FORMAT_PGM;
- else if (strcmp(optarg, O_FORMAT_PGMASC_STR) == 0)
- opt_format = O_FORMAT_PGMASC;
-#if HAVE_PNG
- else if (strcmp(optarg, O_FORMAT_PNG_STR) == 0)
- opt_format = O_FORMAT_PNG;
- else if (strcmp(optarg, O_FORMAT_PNG16_STR) == 0)
- opt_format = O_FORMAT_PNG16;
-#endif
-#if HAVE_GIF
- else if (strcmp(optarg, O_FORMAT_GIF_STR) == 0)
- opt_format = O_FORMAT_GIF;
-#endif
- else if (strcmp(optarg, O_FORMAT_DISP_STR) == 0)
- opt_format = O_FORMAT_DISP;
- else {
- sys_error (ERR_SEVERE, "Invalid format mode %s", optarg);
- if2img_usage(argv[0]);
- return (1);
- }
- break;
- case O_LABELS:
- opt_labels = 1;
- break;
- case O_VERBOSE:
- opt_verbose = 1;
- break;
- case O_DEBUG:
- opt_debug = 1;
- break;
- case O_STATS:
- opt_stats = 1;
- break;
- case O_VERSION:
-#ifdef VERSION
- cout << "Version " << VERSION;
-#else
- cout << "Unknown version number" << endl;
-#endif
- return (0);
- case O_HELP:
- case '?':
- if2img_usage(argv[0]);
- return (0);
- default:
- if2img_usage(argv[0]);
- return (1);
- }
- }
-
- if ((opt_format == O_FORMAT_DISP && optind + 1 != argc)
- || (opt_format != O_FORMAT_DISP && optind + 2 != argc))
- {
- if2img_usage(argv[0]);
- return (1);
- }
-
- in_file = argv[optind];
-
- if (opt_format != O_FORMAT_DISP)
- out_file = argv[optind+1];
- else out_file = NULL;
-
- pim = new ImageFile ();
- ImageFile& im = *pim;
- if (! im.fileRead(in_file)) {
- sys_error (ERR_FATAL, "File %s does not exist", in_file);
- return (1);
- }
-
- if (opt_labels)
- im.printLabels(cout);
-
- if (opt_stats || (! (opt_set_max && opt_set_min))) {
- double min, max, mean, mode, median, stddev;
- double window;
- im.statistics(min, max, mean, mode, median, stddev);
-
- if (opt_auto == O_AUTO_FULL) {
- if (! opt_set_max)
- densmax = max;
- if (! opt_set_min)
- densmin = min;
- }
- if (opt_stats || opt_auto != O_AUTO_FULL) {
-
- if (opt_auto == O_AUTO_FULL)
- ;
- else if (opt_auto == O_AUTO_STD1)
- window = stddev;
- else if (opt_auto == O_AUTO_STD0_1)
- window = stddev * 0.1;
- else if (opt_auto == O_AUTO_STD0_5)
- window = stddev * 0.5;
- else if (opt_auto == O_AUTO_STD2)
- window = stddev * 2;
- else if (opt_auto == O_AUTO_STD3)
- window = stddev * 3;
- else {
- sys_error (ERR_SEVERE, "Internal Error: Invalid auto mode %d", opt_auto);
- return (1);
- }
- }
- if (opt_stats) {
- cout <<"nx: " << im.nx() << endl;
- cout <<"ny: " << im.ny() << endl;
- cout <<"min: " << min << endl;
- cout <<"max: " << max << endl;
- cout <<"mean: " << mean << endl;
- cout <<"mode: " << mode << endl;
- cout <<"stddev: " << stddev << endl;
- }
- if (opt_auto != O_AUTO_FULL) {
- double center;
-
- if (opt_center == O_CENTER_MODE)
- center = mode;
- else if (opt_center == O_CENTER_MEAN)
- center = mean;
- else {
- sys_error (ERR_SEVERE, "Internal Error: Invalid center mode %d", opt_center);
- return (1);
- }
- if (! opt_set_max)
- densmax = center + window;
- if (! opt_set_min)
- densmin = center - window;
- }
- }
-
- if (opt_stats) {
- cout << "min display: " << densmin << endl;
- cout << "max display: " << densmax << endl;
- }
-
- if (opt_format == O_FORMAT_PGM)
- im.writeImagePGM (out_file, opt_scale, opt_scale, densmin, densmax);
- else if (opt_format == O_FORMAT_PGMASC)
- im.writeImagePGMASCII (out_file, opt_scale, opt_scale, densmin, densmax);
-#if HAVE_PNG
- else if (opt_format == O_FORMAT_PNG)
- im.writeImagePNG (out_file, 8, opt_scale, opt_scale, densmin, densmax);
- else if (opt_format == O_FORMAT_PNG16)
- im.writeImagePNG (out_file, 16, opt_scale, opt_scale, densmin, densmax);
-#endif
-#if HAVE_GIF
- else if (opt_format == O_FORMAT_GIF)
- im.writeImageGIF (out_file, opt_scale, opt_scale, densmin, densmax);
-#endif
- else if (opt_format == O_FORMAT_DISP) {
-#if HAVE_SGP
- im.displayScaling (opt_scale, densmin, densmax);
- cout << "Press enter to continue\n";
- cio_kb_getc();
- sgp2_close(sgp2_get_active_win());
-#endif
- }
- else
- {
- sys_error (ERR_SEVERE, "Internal Error: Invalid format mode %d", opt_format);
- return (1);
- }
- return (0);
-}
-
-
-#ifndef NO_MAIN
-int
-main (int argc, char *const argv[])
-{
- int retval = 1;
-
- try {
- retval = if2img_main(argc, argv);
- } catch (exception e) {
- cerr << "Exception: " << e.what() << endl;
- } catch (...) {
- cerr << "Unknown exception" << endl;
- }
-
- return (retval);
-}
-#endif
+++ /dev/null
-/*****************************************************************************
-** FILE IDENTIFICATION
-**
-** Name: ifinfo.cpp
-** Purpose: Display information about an image file
-** Programmer: Kevin Rosenberg
-** Date Started: April 2000
-**
-** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
-**
-** $Id: ifinfo.cpp,v 1.13 2000/07/09 08:16:18 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
-******************************************************************************/
-
-/* FILE
- * ifinfo.c Display info on sdf files
- */
-
-#include "ct.h"
-
-enum { O_LABELS, O_STATS, O_NO_STATS, O_NO_LABELS, O_VERBOSE, O_HELP, O_VERSION, O_DEBUG };
-
-static struct option my_options[] =
-{
- {"labels", 0, 0, O_LABELS},
- {"no-labels", 0, 0, O_NO_LABELS},
- {"stats", 0, 0, O_STATS},
- {"no-stats", 0, 0, O_NO_STATS},
- {"debug", 0, 0, O_DEBUG},
- {"verbose", 0, 0, O_VERBOSE},
- {"help", 0, 0, O_HELP},
- {"version", 0, 0, O_VERSION},
- {0, 0, 0, 0}
-};
-
-void
-ifinfo_usage (const char *program)
-{
- cout << "usage: " << fileBasename(program) << " image-filename [OPTIONS]" << endl;
- cout << "Imagefile information" << endl;
- cout << endl;
- cout << " infile Name of input IF file" << endl;
- cout << " --display Display image" << endl;
- cout << " --labels Print image labels (default)" << endl;
- cout << " --no-labels Do not print image labels" << endl;
- cout << " --stats Print image statistics (default)" << endl;
- cout << " --no-stats Do not print image statistics" << endl;
- cout << " --debug Debug mode" << endl;
- cout << " --verbose Verbose mode" << endl;
- cout << " --version Print version" << endl;
- cout << " --help Print this help message" << endl;
-}
-
-int
-ifinfo_main (int argc, char *const argv[])
-{
- ImageFile *im = NULL;
- string in_file;
- int opt_verbose = 0;
- int opt_stats = 1;
- int opt_labels = 1;
- int opt_debug = 0;
-
- while (1)
- {
- int c = getopt_long (argc, argv, "", my_options, NULL);
-
- if (c == -1)
- break;
-
- switch (c)
- {
- case O_LABELS:
- opt_labels = 1;
- break;
- case O_STATS:
- opt_stats = 1;
- break;
- case O_NO_LABELS:
- opt_labels = 0;
- break;
- case O_NO_STATS:
- opt_stats = 0;
- break;
- case O_VERBOSE:
- opt_verbose = 1;
- break;
- case O_DEBUG:
- opt_debug = 0;
- break;
- case O_VERSION:
-#ifdef VERSION
- cout << "Version " << VERSION << endl;
-#else
- cout << "Unknown version number" << endl;
-#endif
- return (0);
- case O_HELP:
- case '?':
- ifinfo_usage(argv[0]);
- return (0);
- default:
- ifinfo_usage(argv[0]);
- return (1);
- }
- }
-
- if (optind + 1 != argc) {
- ifinfo_usage (argv[0]);
- return (1);
- }
-
- in_file = argv[optind];
-
- im = new ImageFile ();
- if (! im->fileRead (in_file)) {
- sys_error (ERR_WARNING, "Unable to read file %s", in_file.c_str());
- return (1);
- }
-
- if (opt_labels)
- im->printLabels (cout);
-
- if (opt_stats) {
- cout << "Size: (" << im->nx() << "," << im->ny() << ")" << endl;
- im->printStatistics (cout);
- }
-
- return (0);
-}
-
-#ifndef NO_MAIN
-int
-main (int argc, char *const argv[])
-{
- int retval = 1;
-
- try {
- retval = ifinfo_main(argc, argv);
- } catch (exception e) {
- cerr << "Exception: " << e.what() << endl;
- } catch (...) {
- cerr << "Unknown exception" << endl;
- }
-
- return (retval);
-}
-#endif
+++ /dev/null
-/*****************************************************************************
-** FILE IDENTIFICATION
-**
-** Name: phm2if.cpp
-** Purpose: Convert an phantom object to an image file
-** Programmer: Kevin Rosenberg
-** Date Started: 1984
-**
-** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
-**
-** $Id: phm2if.cpp,v 1.16 2000/07/04 22:21:01 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_NSAMPLE, O_FILTER, O_TRACE, O_VERBOSE, O_HELP,
- O_PHMFILE, O_FILTER_DOMAIN, O_FILTER_BW, O_FILTER_PARAM, O_DEBUG, O_VERSION };
-
-static struct option my_options[] =
-{
- {"phantom", 1, 0, O_PHANTOM},
- {"phmfile", 1, 0, O_PHMFILE},
- {"desc", 1, 0, O_DESC},
- {"nsample", 1, 0, O_NSAMPLE},
- {"filter", 1, 0, O_FILTER},
- {"filter-domain", 1, 0, O_FILTER_DOMAIN},
- {"filter-bw", 1, 0, O_FILTER_BW},
- {"filter-param", 1, 0, O_FILTER_PARAM},
- {"trace", 1, 0, O_TRACE},
- {"verbose", 0, 0, O_VERBOSE},
- {"debug", 0, 0, O_DEBUG},
- {"help", 0, 0, O_HELP},
- {"version", 0, 0, O_VERSION},
- {0, 0, 0, 0}
-};
-
-void
-phm2if_usage (const char *program)
-{
- cout << "phm2if_usage: " << fileBasename(program) << " outfile nx ny [--phantom phantom-name] [--phmfile filename] [--filter filter-name] [OPTIONS]" << endl;
- cout << "Generate phantom image from a predefined --phantom or a --phmfile" << endl;
- cout << endl;
- cout << " outfile Name of output file for image" << endl;
- cout << " nx Number of pixels X-axis" << endl;
- cout << " ny Number of pixels Y-axis" << endl;
- cout << " --phantom Phantom to use for projection" << endl;
- cout << " herman Herman head phantom" << endl;
- cout << " bherman Bordered 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 Generate Phantom from phantom file" << endl;
- cout << " --filter Generate Phantom from a filter function" << endl;
- cout << " abs_bandlimit Abs * Bandlimiting" << endl;
- cout << " abs_sinc Abs * Sinc" << endl;
- cout << " abs_cos Abs * Cosine" << endl;
- cout << " abs_hamming Abs * Hamming" << endl;
- cout << " shepp Shepp-Logan" << endl;
- cout << " bandlimit Bandlimiting" << endl;
- cout << " sinc Sinc" << endl;
- cout << " cos Cosine" << endl;
- cout << " triangle Triangle" << endl;
- cout << " hamming Hamming" << endl;
- cout << " --filter-param Alpha level for Hamming filter" << endl;
- cout << " --filter-domain Set domain of filter" << endl;
- cout << " spatial Spatial domain (default)" << endl;
- cout << " freq Frequency domain" << endl;
- cout << " --filter-bw Filter bandwidth (default = 1)" << endl;
- cout << " --desc Description of raysum" << endl;
- cout << " --nsample Number of samples per axis per pixel (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" << endl;
- cout << " rays Trace rays" << endl;
- cout << " plot Trace plot" << endl;
- cout << " clipping Trace clipping" << endl;
- cout << " --debug Debug mode" << endl;
- cout << " --verbose Verbose mode" << endl;
- cout << " --version Print version" << endl;
- cout << " --help Print this help message" << endl;
-}
-
-#ifdef HAVE_MPI
-void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* imGlobal, ImageFile* imLocal, const int opt_debug);
-#endif
-
-int
-phm2if_main (int argc, char* argv[])
-{
- ImageFile* imGlobal = NULL;
- Phantom phm;
- int opt_nx = 0, opt_ny = 0;
- int opt_nsample = 1;
- string optPhmName = Phantom::PHM_HERMAN_STR;
- string optFilterName;
- string optDomainName = SignalFilter::DOMAIN_SPATIAL_STR;
- char *opt_outfile = NULL;
- int opt_debug = 0;
- string opt_desc;
- string opt_phmFileName;
- char *endstr, *endptr;
- double opt_filter_param = 1;
- double opt_filter_bw = 1.;
- int opt_trace = TRACE_NONE;
- bool opt_verbose = false;
-#ifdef HAVE_MPI
- ImageFile* imLocal = NULL;
- MPIWorld mpiWorld (argc, argv);
-#endif
-
- Timer timerProgram;
-
-#ifdef HAVE_MPI
- if (mpiWorld.getRank() == 0) {
-#endif
- while (1) {
- int c = getopt_long(argc, argv, "", my_options, NULL);
- char *endptr = NULL;
- char *endstr;
-
- if (c == -1)
- break;
-
- switch (c) {
- case O_PHANTOM:
- optPhmName = optarg;
- break;
- case O_PHMFILE:
- opt_phmFileName = optarg;
- break;
- case O_VERBOSE:
- opt_verbose = true;
- break;
- case O_DEBUG:
- opt_debug = 1;
- break;
- case O_TRACE:
- if ((opt_trace = convertTraceNameToID(optarg)) == TRACE_INVALID) {
- phm2if_usage(argv[0]);
- return (1);
- }
- break;
- case O_FILTER:
- optFilterName = optarg;
- break;
- case O_FILTER_DOMAIN:
- optDomainName = optarg;
- break;
- case O_DESC:
- opt_desc = optarg;
- break;
- case O_FILTER_BW:
- opt_filter_bw = strtod(optarg, &endptr);
- endstr = optarg + strlen(optarg);
- if (endptr != endstr) {
- sys_error(ERR_SEVERE,"Error setting --filter-bw to %s\n", optarg);
- phm2if_usage(argv[0]);
- return (1);
- }
- break;
- case O_FILTER_PARAM:
- opt_filter_param = strtod(optarg, &endptr);
- endstr = optarg + strlen(optarg);
- if (endptr != endstr) {
- sys_error(ERR_SEVERE,"Error setting --filter-param to %s\n", optarg);
- phm2if_usage(argv[0]);
- return (1);
- }
- break;
- case O_NSAMPLE:
- opt_nsample = strtol(optarg, &endptr, 10);
- endstr = optarg + strlen(optarg);
- if (endptr != endstr) {
- sys_error(ERR_SEVERE,"Error setting --nsample to %s\n", optarg);
- phm2if_usage(argv[0]);
- return (1);
- }
- break;
- case O_VERSION:
-#ifdef VERSION
- cout << "Version " << VERSION << endl;
-#else
- cerr << "Unknown version number" << endl;
-#endif
- case O_HELP:
- case '?':
- phm2if_usage(argv[0]);
- return (0);
- default:
- phm2if_usage(argv[0]);
- return (1);
- }
- }
-
- if (optPhmName == "" && optFilterName == "" && opt_phmFileName == "") {
- cerr << "No phantom defined" << endl << endl;
- phm2if_usage(argv[0]);
- return (1);
- }
-
- if (optind + 3 != argc) {
- phm2if_usage(argv[0]);
- return (1);
- }
- opt_outfile = argv[optind];
- opt_nx = strtol(argv[optind+1], &endptr, 10);
- endstr = argv[optind+1] + strlen(argv[optind+1]);
- if (endptr != endstr) {
- sys_error(ERR_SEVERE,"Error setting nx to %s\n", argv[optind+1]);
- phm2if_usage(argv[0]);
- return (1);
- }
- opt_ny = strtol(argv[optind+2], &endptr, 10);
- endstr = argv[optind+2] + strlen(argv[optind+2]);
- if (endptr != endstr) {
- sys_error(ERR_SEVERE,"Error setting ny to %s\n", argv[optind+2]);
- phm2if_usage(argv[0]);
- return (1);
- }
-
- ostringstream oss;
- oss << "phm2if: nx=" << opt_nx << ", ny=" << opt_ny << ", nsample=" << opt_nsample << ", ";
- if (opt_phmFileName != "")
- oss << "phantom=" << opt_phmFileName;
- else if (optPhmName != "")
- oss << "phantom=" << optPhmName;
- else if (optFilterName != "") {
- oss << "filter=" << optFilterName << " - " << optDomainName;
- }
- if (opt_desc != "")
- oss << ": " << opt_desc;
- opt_desc = oss.str();
-
- if (optPhmName != "") {
- phm.createFromPhantom (optPhmName.c_str());
- if (phm.fail()) {
- cout << phm.failMessage() << endl << endl;
- phm2if_usage(argv[0]);
- return (1);
- }
- }
-
- if (opt_phmFileName != "") {
- phm.createFromFile(opt_phmFileName.c_str());
-#ifdef HAVE_MPI
- if (mpiWorld.getRank() == 0)
- cerr << "Can't use phantom from file in MPI mode" << endl;
- return (1);
-#endif
- }
-
- if (opt_verbose)
- cout << "Rasterize Phantom to Image" << endl << endl;
-#ifdef HAVE_MPI
- }
-#endif
-
-#ifdef HAVE_MPI
- TimerCollectiveMPI timerBcast (mpiWorld.getComm());
- mpiWorld.BcastString (optPhmName);
- 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);
- mpiWorld.getComm().Bcast (&opt_nx, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&opt_ny, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&opt_nsample, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&opt_filter_param, 1, MPI::DOUBLE, 0);
- mpiWorld.getComm().Bcast (&opt_filter_bw, 1, MPI::DOUBLE, 0);
-
- mpiWorld.BcastString (optFilterName);
- mpiWorld.BcastString (optDomainName);
-
- if (opt_verbose)
- timerBcast.timerEndAndReport ("Time to broadcast variables");
-
- mpiWorld.setTotalWorkUnits (opt_nx);
-
- if (mpiWorld.getRank() > 0 && optPhmName != "")
- phm.createFromPhantom (optPhmName.c_str());
-
- if (mpiWorld.getRank() == 0) {
- imGlobal = new ImageFile (opt_nx, opt_ny);
- }
- imLocal = new ImageFile (opt_nx, opt_ny);
-#else
- imGlobal = new ImageFile (opt_nx, opt_ny);
-#endif
-
- ImageFileArray v;
-#ifdef HAVE_MPI
- if (mpiWorld.getRank() == 0)
- v = imGlobal->getArray ();
-
- if (phm.getComposition() == P_UNIT_PULSE) {
- if (mpiWorld.getRank() == 0) {
- v[opt_nx/2][opt_ny/2] = 1.;
- }
- } else if (optFilterName != "") {
- if (mpiWorld.getRank() == 0) {
- imGlobal->filterResponse (optDomainName.c_str(), opt_filter_bw, optFilterName.c_str(), opt_filter_param);
- }
- } else {
- TimerCollectiveMPI timerRasterize (mpiWorld.getComm());
- phm.convertToImagefile (*imLocal, opt_nsample, opt_trace, mpiWorld.getMyStartWorkUnit(), mpiWorld.getMyLocalWorkUnits());
- if (opt_verbose)
- timerRasterize.timerEndAndReport ("Time to rasterize phantom");
-
- TimerCollectiveMPI timerGather (mpiWorld.getComm());
- mpi_gather_image (mpiWorld, imGlobal, imLocal, opt_debug);
- if (opt_verbose)
- timerGather.timerEndAndReport ("Time to gather image");
- }
-#else
- v = imGlobal->getArray ();
- if (phm.getComposition() == P_UNIT_PULSE) {
- v[opt_nx/2][opt_ny/2] = 1.;
- } else if (optFilterName != "") {
- imGlobal->filterResponse (optDomainName.c_str(), opt_filter_bw, optFilterName.c_str(), opt_filter_param);
- } else {
-#if HAVE_SGP
- if (opt_trace >= TRACE_PHM)
- phm.show();
-#endif
- phm.convertToImagefile (*imGlobal, opt_nsample, opt_trace);
- }
-#endif
-
-#ifdef HAVE_MPI
- if (mpiWorld.getRank() == 0)
-#endif
- {
- double calctime = timerProgram.timerEnd ();
- imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, opt_desc.c_str(), calctime);
- imGlobal->fileWrite (opt_outfile);
- if (opt_verbose)
- cout << "Time to rasterized phantom: " << calctime << " seconds" << endl;
-
- if (opt_trace >= TRACE_PHM) {
- double dmin, dmax;
- int nscale;
-
- printf ("Enter display size scale (nominal = 1): ");
- scanf ("%d", &nscale);
- printf ("Enter minimum and maximum densities (min, max): ");
- scanf ("%lf %lf", &dmin, &dmax);
- imGlobal->displayScaling (nscale, dmin, dmax);
- }
- }
-
- return (0);
-}
-
-
-
-#ifdef HAVE_MPI
-void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* imGlobal, ImageFile* imLocal, const int opt_debug)
-{
- ImageFileArray vLocal = imLocal->getArray();
- ImageFileArray vGlobal = NULL;
- int nyLocal = imLocal->ny();
-
- if (mpiWorld.getRank() == 0)
- vGlobal = imGlobal->getArray();
-
- for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++)
- mpiWorld.getComm().Send(vLocal[iw], nyLocal, imLocal->getMPIDataType(), 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;
- mpiWorld.getComm().Recv(vGlobal[iw], nyLocal, imLocal->getMPIDataType(), iProc, 0, status);
- }
- }
- }
-
-}
-#endif
-
-#ifndef NO_MAIN
-int
-main (int argc, char* argv[])
-{
- int retval = 1;
-
- try {
- retval = phm2if_main(argc, argv);
- } catch (exception e) {
- cerr << "Exception: " << e.what() << endl;
- } catch (...) {
- cerr << "Unknown exception" << endl;
- }
-
- return (retval);
-}
-#endif
+++ /dev/null
-/*****************************************************************************
-** 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.6 2000/07/04 22:21:01 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 << " bherman Bordered 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;
- string optGeometryName = Scanner::GEOMETRY_PARALLEL_STR;
- char *opt_outfile = NULL;
- string opt_desc;
- string optPhmFileName;
- int opt_ndet;
- int opt_nview;
- int opt_nray = 1;
- int opt_trace = 0;
- string optPhmName = Phantom::PHM_HERMAN_STR;
- 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:
- optPhmName = optarg;
- break;
- case O_PHMFILE:
- optPhmFileName = optarg;
- break;
- case O_VERBOSE:
- opt_verbose = 1;
- break;
- case O_DEBUG:
- opt_debug = 1;
- break;
- break;
- case O_TRACE:
- if ((opt_trace = convertTraceNameToID(optarg)) == TRACE_INVALID) {
- 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
- return (0);
- case O_HELP:
- case '?':
- phm2pj_usage(argv[0]);
- return (0);
- default:
- phm2pj_usage(argv[0]);
- return (1);
- }
- }
-
- if (optPhmName == "" && optPhmFileName == "") {
- cerr << "No phantom defined" << endl << 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 << "phm2pj: NDet=" << opt_ndet << ", Nview=" << opt_nview << ", NRay=" << opt_nray << ", RotAngle=" << opt_rotangle << ", ";
- if (optPhmFileName.length()) {
- desc << "PhantomFile=" << optPhmFileName;
- } else if (optPhmName != "") {
- desc << "Phantom=" << optPhmName;
- }
- if (opt_desc.length()) {
- desc << ": " << opt_desc;
- }
- opt_desc = desc.str();
-
- if (optPhmName != "") {
- phm.createFromPhantom (optPhmName.c_str());
- if (phm.fail()) {
- cout << phm.failMessage() << endl << endl;
- phm2pj_usage(argv[0]);
- return (1);
- }
- }
-
- if (optPhmFileName != "") {
-#ifdef HAVE_MPI
- cerr << "Can not read phantom from file in MPI mode" << endl;
- return (1);
-#endif
- phm.createFromFile (optPhmFileName.c_str());
- }
-
-#ifdef HAVE_MPI
- }
-#endif
-
-#ifdef HAVE_MPI
- TimerCollectiveMPI timerBcast(mpiWorld.getComm());
- mpiWorld.BcastString (optPhmName);
- 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_verbose, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&opt_debug, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&opt_trace, 1, MPI::INT, 0);
- timerBcast.timerEndAndReport ("Time to broadcast variables");
-
- if (mpiWorld.getRank() > 0 && optPhmName != "")
- phm.createFromPhantom (optPhmName.c_str());
-#endif
-
- opt_rotangle *= PI;
- Scanner scanner (phm, optGeometryName.c_str(), opt_ndet, opt_nview, opt_nray, opt_rotangle);
- if (scanner.fail()) {
- cout << "Scanner Creation Error: " << scanner.failMessage() << endl;
- return (1);
- }
-#ifdef HAVE_MPI
- mpiWorld.setTotalWorkUnits (opt_nview);
-
- Projections pjGlobal;
- if (mpiWorld.getRank() == 0)
- pjGlobal.initFromScanner (scanner);
-
- if (opt_verbose)
- pjGlobal.printScanInfo();
-
- Projections pjLocal (scanner);
- pjLocal.setNView (mpiWorld.getMyLocalWorkUnits());
-
- if (opt_debug)
- cout << "pjLocal->nview = " << pjLocal.nView() << " (process " << mpiWorld.getRank() << ")" << endl;;
-
- TimerCollectiveMPI timerProject (mpiWorld.getComm());
- scanner.collectProjections (pjLocal, phm, mpiWorld.getMyStartWorkUnit(), opt_trace);
- if (opt_verbose)
- timerProject.timerEndAndReport ("Time to collect projections");
-
- TimerCollectiveMPI 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[])
-{
- int retval = 1;
-
- try {
- retval = phm2pj_main(argc, argv);
- } catch (exception e) {
- cerr << "Exception: " << e.what() << endl;
- } catch (...) {
- cerr << "Unknown exception" << endl;
- }
-
- return (retval);
-}
-#endif
-
+++ /dev/null
-/*****************************************************************************
-** FILE IDENTIFICATION
-**
-** Name: pj2if.cpp
-** Purpose: Convert an projection data file to an image file
-** Programmer: Kevin Rosenberg
-** Date Started: April 2000
-**
-** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
-**
-** $Id: pj2if.cpp,v 1.6 2000/06/28 15:25:34 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
-******************************************************************************/
-
-/* FILE
- * pj2if.c Convert Raysum to image
- *
- * DATE
- * Apr 1999
- */
-
-#include "ct.h"
-#include "timer.h"
-
-
-enum { O_VERBOSE, O_HELP, O_VERSION };
-
-static struct option my_options[] =
-{
- {"verbose", 0, 0, O_VERBOSE},
- {"help", 0, 0, O_HELP},
- {"version", 0, 0, O_VERSION},
- {0, 0, 0, 0}
-};
-
-
-void
-pj2if_usage (const char *program)
-{
- cout << "usage: " << fileBasename(program) << " in-proj-file out-if-file [OPTIONS]" << endl;
- cout << "Converts a projection file to a IF file" << endl;
- cout << endl;
- cout << " --verbose Verbose mode" << endl;
- cout << " --version Print version" << endl;
- cout << " --help Print this help message" << endl;
-}
-
-
-
-int
-pj2if_main (const int argc, char *const argv[])
-{
- char *pj_name, *im_name;
- bool optVerbose = false;
- extern int optind;
- Timer timerProgram;
-
- while (1)
- {
- int c = getopt_long (argc, argv, "", my_options, NULL);
- if (c == -1)
- break;
-
- switch (c)
- {
- case O_VERBOSE:
- optVerbose = true;
- break;
- case O_VERSION:
-#ifdef VERSION
- cout << "Version " << VERSION << endl;
-#else
- cout << "Unknown version number" << endl;
-#endif
- return (0);
- case O_HELP:
- case '?':
- pj2if_usage(argv[0]);
- return (0);
- default:
- pj2if_usage(argv[0]);
- return (1);
- }
- }
-
- if (argc - optind != 2) {
- pj2if_usage(argv[0]);
- return (1);
- }
-
- pj_name = argv[optind];
- im_name = argv[optind + 1];
-
- Projections pj;
- if (! pj.read (pj_name)) {
- sys_error (ERR_SEVERE, "Can not open projection file %s", pj_name);
- return (1);
- }
-
- if (optVerbose)
- pj.printScanInfo();
-
- ImageFile im (pj.nDet(), pj.nView());
-
- ImageFileArray v = im.getArray();
-
- for (int iy = 0; iy < pj.nView(); iy++)
- {
- DetectorArray& detarray = pj.getDetectorArray (iy);
- const DetectorValue* detval = detarray.detValues();
- for (int ix = 0; ix < pj.nDet(); ix++)
- {
- v[ix][iy] = detval[ix];
- }
- }
-
- im.labelAdd (pj.getLabel());
- im.labelAdd (Array2dFileLabel::L_HISTORY, "Conversion from .pj to .if", timerProgram.timerEnd());
- im.fileWrite (im_name);
-
- return(0);
-}
-
-
-#ifndef NO_MAIN
-int
-main (const int argc, char *const argv[])
-{
- int retval = 1;
-
- try {
- retval = pj2if_main(argc, argv);
- } catch (exception e) {
- cerr << "Exception: " << e.what() << endl;
- } catch (...) {
- cerr << "Unknown exception" << endl;
- }
-
- return (retval);
-}
-#endif
+++ /dev/null
-/*****************************************************************************
-** FILE IDENTIFICATION
-**
-** Name: pjrec.cpp
-** Purpose: Reconstruct an image from projections
-** Programmer: Kevin Rosenberg
-** Date Started: Aug 1984
-**
-** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
-**
-** $Id: pjrec.cpp,v 1.10 2000/07/11 10:32:44 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_INTERP, O_FILTER, O_FILTER_METHOD, O_ZEROPAD, O_FILTER_PARAM, O_BACKPROJ, O_PREINTERPOLATION_FACTOR, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION};
-
-static struct option my_options[] =
-{
- {"interp", 1, 0, O_INTERP},
- {"preinterpolation-factor", 1, 0, O_PREINTERPOLATION_FACTOR},
- {"filter", 1, 0, O_FILTER},
- {"filter-method", 1, 0, O_FILTER_METHOD},
- {"zeropad", 1, 0, O_ZEROPAD},
- {"filter-param", 1, 0, O_FILTER_PARAM},
- {"backproj", 1, 0, O_BACKPROJ},
- {"trace", 1, 0, O_TRACE},
- {"debug", 0, 0, O_DEBUG},
- {"verbose", 0, 0, O_VERBOSE},
- {"help", 0, 0, O_HELP},
- {"version", 0, 0, O_VERSION},
- {0, 0, 0, 0}
-};
-
-
-void
-pjrec_usage (const char *program)
-{
- cout << "usage: " << fileBasename(program) << " raysum-file image-file nx-image ny-image [OPTIONS]" << endl;
- cout << "Image reconstruction from raysum projections" << endl;
- cout << endl;
- cout << " raysum-file Input raysum file" << endl;
- cout << " image-file Output image file in SDF2D format" << endl;
- cout << " nx-image Number of columns in output image" << endl;
- cout << " ny-image Number of rows in output image" << endl;
- cout << " --interp Interpolation method during backprojection" << endl;
- cout << " nearest Nearest neighbor interpolation" << endl;
- cout << " linear Linear interpolation" << endl;
-#if HAVE_BSPLINE_INTERP
- cout << " bspline B-spline interpolation" << endl;
-#endif
- cout << " --preinterpolate Preinterpolation factor (default = 1)\n";
- cout << " Used only with frequency-based filtering\n";
- cout << " --filter Filter name" << endl;
- cout << " abs_bandlimit Abs * Bandlimiting (default)" << endl;
- cout << " abs_sinc Abs * Sinc" << endl;
- cout << " abs_cos Abs * Cosine" << endl;
- cout << " abs_hamming Abs * Hamming" << endl;
- cout << " shepp Shepp-Logan" << endl;
- cout << " bandlimit Bandlimiting" << endl;
- cout << " sinc Sinc" << endl;
- cout << " cos Cosine" << endl;
- cout << " triangle Triangle" << endl;
- cout << " hamming Hamming" << endl;
- cout << " --filter-method Filter method before backprojections\n";;
- cout << " convolution Spatial filtering (default)\n";
- cout << " fourier Frequency filtering with discete fourier\n";
- cout << " fourier_table Frequency filtering with table lookup fourier\n";
- cout << " fft Fast Fourier Transform\n";
-#if HAVE_FFTW
- cout << " fftw Fast Fourier Transform West library\n";
- cout << " rfftw Fast Fourier Transform West (real-mode) library\n";
-#endif
- cout << " --zeropad n Set zeropad level (default = 0)\n";
- cout << " set n to number of powers to two to pad\n";
- cout << " --backproj Backprojection Method" << endl;
- cout << " trig Trigometric functions at every point" << endl;
- cout << " table Trigometric functions with precalculated table" << endl;
- cout << " diff Difference method" << endl;
- cout << " diff2 Optimized difference method (default)" << endl;
- cout << " idiff2 Optimized difference method with integer math" << endl;
- cout << " idiff3 Highly-optimized difference method with integer math" << endl;
- cout << " --filter-param Alpha level for Hamming filter" << endl;
- cout << " --trace Set tracing to level" << endl;
- cout << " none No tracing (default)" << endl;
- cout << " text Text level tracing" << endl;
- cout << " phm Trace phantom" << endl;
- cout << " rays Trace allrays" << endl;
- cout << " plot Trace plotting" << endl;
- cout << " clipping Trace clipping" << endl;
- cout << " --verbose Turn on verbose mode" << endl;
- cout << " --debug Turn on debug mode" << endl;
- cout << " --version Print version" << endl;
- cout << " --help Print this help message" << endl;
-}
-
-
-#ifdef HAVE_MPI
-static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int debug);
-static void ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* imGlobal);
-#endif
-
-
-int
-pjrec_main (int argc, char * argv[])
-{
- Projections projGlobal;
- ImageFile* imGlobal = NULL;
- char* filenameProj = NULL;
- char* filenameImage = NULL;
- string remark;
- char *endptr;
- int optVerbose = 0;
- int optDebug = 0;
- int optZeroPad = 0;
- int optTrace = TRACE_NONE;
- double optFilterParam = -1;
- string optFilterName = SignalFilter::FILTER_ABS_BANDLIMIT_STR;
- string optFilterMethodName = SignalFilter::FILTER_METHOD_CONVOLUTION_STR;
- string optInterpName = Backprojector::INTERP_LINEAR_STR;
- string optBackprojName = Backprojector::BPROJ_IDIFF2_STR;
- // string optFilterMethodName = SignalFilter::FILTER_METHOD_CONVOLUTION_STR;
- int optPreinterpolationFactor = 1;
- int nx, ny;
-#ifdef HAVE_MPI
- ImageFile* imLocal;
- int mpi_nview, mpi_ndet;
- double mpi_detinc, mpi_rotinc, mpi_phmlen;
- MPIWorld mpiWorld (argc, argv);
-#endif
-
- Timer timerProgram;
-
-#ifdef HAVE_MPI
- if (mpiWorld.getRank() == 0) {
-#endif
- while (1) {
- int c = getopt_long(argc, argv, "", my_options, NULL);
- char *endptr = NULL;
-
- if (c == -1)
- break;
-
- switch (c)
- {
- case O_INTERP:
- optInterpName = optarg;
- break;
- case O_PREINTERPOLATION_FACTOR:
- optPreinterpolationFactor = strtol(optarg, &endptr, 10);
- if (endptr != optarg + strlen(optarg)) {
- pjrec_usage(argv[0]);
- return(1);
- }
- break;
- case O_FILTER:
- optFilterName = optarg;
- break;
- case O_FILTER_METHOD:
- optFilterMethodName = optarg;
- break;
- case O_BACKPROJ:
- optBackprojName = optarg;
- break;
- case O_FILTER_PARAM:
- optFilterParam = strtod(optarg, &endptr);
- if (endptr != optarg + strlen(optarg)) {
- pjrec_usage(argv[0]);
- return(1);
- }
- break;
- case O_ZEROPAD:
- optZeroPad = strtol(optarg, &endptr, 10);
- if (endptr != optarg + strlen(optarg)) {
- pjrec_usage(argv[0]);
- return(1);
- }
- break;
- case O_VERBOSE:
- optVerbose = 1;
- break;
- case O_DEBUG:
- optDebug = 1;
- break;
- case O_TRACE:
- if ((optTrace = convertTraceNameToID(optarg)) == TRACE_INVALID) {
- pjrec_usage(argv[0]);
- return (1);
- }
- break;
- case O_VERSION:
-#ifdef VERSION
- cout << "Version " << VERSION << endl;
-#else
- cout << "Unknown version number" << endl;
-#endif
- return (0);
- case O_HELP:
- case '?':
- pjrec_usage(argv[0]);
- return (0);
- default:
- pjrec_usage(argv[0]);
- return (1);
- }
- }
-
- if (optind + 4 != argc) {
- pjrec_usage(argv[0]);
- return (1);
- }
-
- filenameProj = argv[optind];
-
- filenameImage = argv[optind + 1];
-
- nx = strtol(argv[optind + 2], &endptr, 10);
- ny = strtol(argv[optind + 3], &endptr, 10);
-
- ostringstream filterDesc;
- if (optFilterParam >= 0)
- filterDesc << optFilterName << ": alpha=" << optFilterParam;
- else
- filterDesc << optFilterName;
-
- ostringstream label;
- label << "pjrec: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << optInterpName << ", preinterpolation=" << optPreinterpolationFactor << ", " << optBackprojName;
- remark = label.str();
-
- if (optVerbose)
- cout << "Remark: " << remark << endl;
-#ifdef HAVE_MPI
- }
-#endif
-
-#ifdef HAVE_MPI
- if (mpiWorld.getRank() == 0) {
- projGlobal.read (filenameProj);
- if (optVerbose)
- projGlobal.printScanInfo();
-
- mpi_ndet = projGlobal.nDet();
- mpi_nview = projGlobal.nView();
- mpi_detinc = projGlobal.detInc();
- mpi_phmlen = projGlobal.phmLen();
- mpi_rotinc = projGlobal.rotInc();
- }
-
- TimerCollectiveMPI timerBcast (mpiWorld.getComm());
- mpiWorld.BcastString (optBackprojName);
- mpiWorld.BcastString (optFilterName);
- mpiWorld.BcastString (optFilterMethodName);
- mpiWorld.BcastString (optInterpName);
- mpiWorld.getComm().Bcast (&optVerbose, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&optDebug, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&optTrace, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&optFilterParam, 1, MPI::DOUBLE, 0);
- mpiWorld.getComm().Bcast (&optZeroPad, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&optPreinterpolationFactor, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&mpi_ndet, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&mpi_nview, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&mpi_detinc, 1, MPI::DOUBLE, 0);
- mpiWorld.getComm().Bcast (&mpi_phmlen, 1, MPI::DOUBLE, 0);
- mpiWorld.getComm().Bcast (&mpi_rotinc, 1, MPI::DOUBLE, 0);
- mpiWorld.getComm().Bcast (&nx, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&ny, 1, MPI::INT, 0);
- if (optVerbose)
- timerBcast.timerEndAndReport ("Time to broadcast variables");
-
- mpiWorld.setTotalWorkUnits (mpi_nview);
-
- Projections projLocal (mpiWorld.getMyLocalWorkUnits(), mpi_ndet);
- projLocal.setDetInc (mpi_detinc);
- projLocal.setPhmLen (mpi_phmlen);
- projLocal.setRotInc (mpi_rotinc);
-
- TimerCollectiveMPI timerScatter (mpiWorld.getComm());
- ScatterProjectionsMPI (mpiWorld, projGlobal, projLocal, optDebug);
- if (optVerbose)
- timerScatter.timerEndAndReport ("Time to scatter projections");
-
- if (mpiWorld.getRank() == 0) {
- imGlobal = new ImageFile (nx, ny);
- }
-
- imLocal = new ImageFile (nx, ny);
-#else
- projGlobal.read (filenameProj);
- if (optVerbose)
- projGlobal.printScanInfo();
-
- imGlobal = new ImageFile (nx, ny);
-#endif
-
-#ifdef HAVE_MPI
- TimerCollectiveMPI timerReconstruct (mpiWorld.getComm());
- projLocal.reconstruct (*imLocal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optPreinterpolationFactor, optBackprojName.c_str(), optTrace);
- if (optVerbose)
- timerReconstruct.timerEndAndReport ("Time to reconstruct");
-
- TimerCollectiveMPI timerReduce (mpiWorld.getComm());
- ReduceImageMPI (mpiWorld, imLocal, imGlobal);
- if (optVerbose)
- timerReduce.timerEndAndReport ("Time to reduce image");
-#else
- projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optPreinterpolationFactor, optBackprojName.c_str(), optTrace);
-#endif
-
-#ifdef HAVE_MPI
- if (mpiWorld.getRank() == 0)
-#endif
- {
- double calcTime = timerProgram.timerEnd();
- imGlobal->labelAdd (projGlobal.getLabel());
- imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, remark.c_str(), calcTime);
- imGlobal->fileWrite (filenameImage);
- if (optVerbose)
- cout << "Run time: " << calcTime << " seconds" << endl;
- }
-#ifdef HAVE_MPI
- MPI::Finalize();
-#endif
-
- return (0);
-}
-
-
-//////////////////////////////////////////////////////////////////////////////////////
-// MPI Support Routines
-//
-//////////////////////////////////////////////////////////////////////////////////////
-
-#ifdef HAVE_MPI
-static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int optDebug)
-{
- if (mpiWorld.getRank() == 0) {
- for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) {
- for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) {
- DetectorArray& detarray = projGlobal.getDetectorArray( iw );
- int nDet = detarray.nDet();
- DetectorValue* detval = detarray.detValues();
-
- double viewAngle = detarray.viewAngle();
- mpiWorld.getComm().Send(&nDet, 1, MPI::INT, iProc, 0);
- mpiWorld.getComm().Send(&viewAngle, 1, MPI::DOUBLE, iProc, 0);
- mpiWorld.getComm().Send(detval, nDet, MPI::FLOAT, iProc, 0);
- }
- }
- }
-
- for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++) {
- MPI::Status status;
- int nDet;
- double viewAngle;
- DetectorValue* detval = projLocal.getDetectorArray(iw).detValues();
-
- mpiWorld.getComm().Recv(&nDet, 1, MPI::INT, 0, 0, status);
- mpiWorld.getComm().Recv(&viewAngle, 1, MPI::DOUBLE, 0, 0, status);
- mpiWorld.getComm().Recv(detval, nDet, MPI::FLOAT, 0, 0, status);
- projLocal.getDetectorArray(iw).setViewAngle( viewAngle );
- }
-}
-
-static void
-ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* imGlobal)
-{
- ImageFileArray vLocal = imLocal->getArray();
-
- for (int ix = 0; ix < imLocal->nx(); ix++) {
- void *recvbuf = NULL;
- if (mpiWorld.getRank() == 0) {
- ImageFileArray vGlobal = imGlobal->getArray();
- recvbuf = vGlobal[ix];
- }
- mpiWorld.getComm().Reduce (vLocal[ix], recvbuf, imLocal->ny(), imLocal->getMPIDataType(), MPI::SUM, 0);
- }
-}
-
-#endif
-
-
-#ifndef NO_MAIN
-int
-main (int argc, char* argv[])
-{
- int retval = 1;
-
- try {
- retval = pjrec_main(argc, argv);
- } catch (exception e) {
- cerr << "Exception: " << e.what() << endl;
- } catch (...) {
- cerr << "Unknown exception" << endl;
- }
-
- return (retval);
-}
-#endif
-
+++ /dev/null
-#!/bin/sh
-
-if test "$1" != "" ; then
- bin=$1
-else
- bin="@prefix@/bin/"
-fi
-
-if test "$1" = "clean" ; then
- rm -f sample-phm.png sample-phm16.png sample-phm.if sample-pj.pj sample-pj.if sample-pj.png sample-pj16.png sample-rec.if sample-rec.png sample-rec16.png
- exit
-fi
-
-# Generate phantom image
-
-${bin}phm2if sample-phm.if 256 256 --nsample 2 --phantom herman
-if [ -f sample-phm.if ] ; then
- ${bin}if2img sample-phm.if sample-phm.png --format png
- ${bin}if2img sample-phm.if sample-phm16.png --format png16
-fi
-
-# Simulate CT data collection and generate raysum sinugram for display
-${bin}phm2pj sample-pj.pj 367 320 --nray 2 --phantom herman
-if [ -f sample-pj.pj ]; then
- ${bin}pj2if sample-pj.pj sample-pj.if
-fi
-if [ -f sample-pj.if ]; then
- ${bin}if2img sample-pj.if sample-pj.png --format png
- ${bin}if2img sample-pj.if sample-pj16.png --format png16
-fi
-
-# Reconstruct raysums and generate image for display
-${bin}pjrec sample-pj.pj sample-rec.if 256 256
-if [ -f sample-rec.if ]; then
- ${bin}if2img sample-rec.if sample-rec.png --format png
- ${bin}if2img sample-rec.if sample-rec16.png --format png16
-
- ${bin}if-2 sample-phm.if sample-rec.if --comp
-fi
-
-# Files sample-phm.png, sample-pj.png, and sample-rec.png are ready for display
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: view.cpp
+** Purpose: View & Canvas routines for CTSim program
+** Programmer: Kevin Rosenberg
+** Date Started: July 2000
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: views.cpp,v 1.1 2000/07/13 07:01:59 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
+******************************************************************************/
+
+#ifdef __GNUG__
+// #pragma implementation
+#endif
+
+// For compilers that support precompilation, includes "wx/wx.h".
+#include "wx/wxprec.h"
+
+#ifdef __BORLANDC__
+#pragma hdrstop
+#endif
+
+#ifndef WX_PRECOMP
+#include "wx/wx.h"
+#endif
+
+#if !wxUSE_DOC_VIEW_ARCHITECTURE
+#error You must set wxUSE_DOC_VIEW_ARCHITECTURE to 1 in setup.h!
+#endif
+
+#include "ctsim.h"
+#include "docs.h"
+#include "views.h"
+#include <sstream>
+#include "ct.h"
+//#include "imagefile.h"
+//#include "phantom.h"
+//#include "scanner.h"
+//#include "projections.h"
+
+// ImageFileCanvas
+
+BEGIN_EVENT_TABLE(ImageFileCanvas, wxScrolledWindow)
+ EVT_MOUSE_EVENTS(ImageFileCanvas::OnMouseEvent)
+END_EVENT_TABLE()
+
+
+ImageFileCanvas::ImageFileCanvas (ImageFileView* v, wxFrame *frame, const wxPoint& pos, const wxSize& size, const long style)
+ : wxScrolledWindow(frame, -1, pos, size, style)
+{
+ m_pView = v;
+}
+
+void
+ImageFileCanvas::OnDraw(wxDC& dc)
+{
+ if (m_pView)
+ m_pView->OnDraw(& dc);
+}
+
+void
+ImageFileCanvas::OnMouseEvent(wxMouseEvent& event)
+{
+ if (! m_pView)
+ return;
+
+ wxClientDC dc(this);
+ PrepareDC(dc);
+
+ dc.SetPen(*wxBLACK_PEN);
+
+ wxPoint pt(event.GetLogicalPosition(dc));
+
+ if (event.LeftUp())
+ cout << pt.x << "x" << pt.y << endl;
+}
+
+
+// ImageFileView
+
+IMPLEMENT_DYNAMIC_CLASS(ImageFileView, wxView)
+
+BEGIN_EVENT_TABLE(ImageFileView, wxView)
+ EVT_MENU(IFMENU_FILE_PROPERTIES, ImageFileView::OnProperties)
+END_EVENT_TABLE()
+
+bool ImageFileView::m_bPColoursInitialized = false;
+wxColour* ImageFileView::m_pColours[256];
+
+ImageFileView::ImageFileView(void)
+ : wxView(), m_canvas(NULL), m_frame(NULL), m_bMinSpecified(false), m_bMaxSpecified(false)
+{
+ if (! m_bPColoursInitialized) {
+ for (int i = 0; i < 256; i++)
+ m_pColours[i] = new wxColour (i, i, i);
+
+ m_bPColoursInitialized = true;
+ }
+}
+
+ImageFileView::~ImageFileView(void)
+{
+}
+
+void
+ImageFileView::OnProperties (wxCommandEvent& event)
+{
+ double min, max, mean, mode, median, stddev;
+ const ImageFile& rIF = dynamic_cast<ImageFileDocument*>(GetDocument())->getImageFile();
+ const string& rFilename = rIF.getFilename();
+ rIF.statistics (min, max, mean, mode, median, stddev);
+ ostringstream os;
+ os << "file: " << rFilename << "\nmin: "<<min<<"\nmax: "<<max<<"\nmean: "<<mean<<"\nmode: "<<mode<<"\nstddev: "<<stddev;
+ wxMessageBox(os.str().c_str(), "Image Properties", wxOK | wxICON_INFORMATION, m_frame);
+}
+
+
+ImageFileCanvas*
+ImageFileView::CreateCanvas (wxView *view, wxFrame *parent)
+{
+ ImageFileCanvas* pCanvas;
+ int width, height;
+ parent->GetClientSize(&width, &height);
+
+ pCanvas = new ImageFileCanvas (dynamic_cast<ImageFileView*>(view), parent, wxPoint(0, 0), wxSize(width, height), 0);
+
+ pCanvas->SetScrollbars(20, 20, 50, 50);
+ pCanvas->SetBackgroundColour(*wxWHITE);
+ pCanvas->Clear();
+
+ return pCanvas;
+}
+
+wxFrame*
+ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view)
+{
+ wxDocChildFrame *subframe = new wxDocChildFrame(doc, view, dynamic_cast<CTSimApp*>(wxTheApp)->getMainFrame(), -1, "ImageFile Frame", wxPoint(10, 10), wxSize(300, 300), wxDEFAULT_FRAME_STYLE);
+
+ wxMenu *file_menu = new wxMenu;
+
+ file_menu->Append(wxID_NEW, "&New...");
+ file_menu->Append(wxID_OPEN, "&Open...");
+ file_menu->Append(wxID_CLOSE, "&Close");
+ file_menu->Append(wxID_SAVE, "&Save");
+ file_menu->Append(wxID_SAVEAS, "Save &As...");
+
+ file_menu->AppendSeparator();
+ file_menu->Append(IFMENU_FILE_PROPERTIES, "P&roperties");
+
+ file_menu->AppendSeparator();
+ file_menu->Append(wxID_PRINT, "&Print...");
+ file_menu->Append(wxID_PRINT_SETUP, "Print &Setup...");
+ file_menu->Append(wxID_PREVIEW, "Print Pre&view");
+
+ wxMenu *help_menu = new wxMenu;
+ help_menu->Append(MAINMENU_HELP_ABOUT, "&About");
+
+ wxMenuBar *menu_bar = new wxMenuBar;
+
+ menu_bar->Append(file_menu, "&File");
+ menu_bar->Append(help_menu, "&Help");
+
+ subframe->SetMenuBar(menu_bar);
+
+ subframe->Centre(wxBOTH);
+
+ return subframe;
+}
+
+
+bool
+ImageFileView::OnCreate(wxDocument *doc, long WXUNUSED(flags) )
+{
+ m_frame = CreateChildFrame(doc, this);
+
+ m_bMinSpecified = false;
+ m_bMaxSpecified = false;
+
+ int width, height;
+ m_frame->GetClientSize(&width, &height);
+ m_frame->SetTitle("ImageFileView");
+ m_canvas = CreateCanvas(this, m_frame);
+
+#ifdef __X__
+ int x, y; // X requires a forced resize
+ m_frame->GetSize(&x, &y);
+ m_frame->SetSize(-1, -1, x, y);
+#endif
+
+ m_frame->Show(true);
+ Activate(true);
+
+ return true;
+}
+
+void
+ImageFileView::OnDraw (wxDC* dc)
+{
+ const ImageFile& rIF = dynamic_cast<ImageFileDocument*>(GetDocument())->getImageFile();
+ dc->Blit(static_cast<wxCoord>(0), static_cast<wxCoord>(0), static_cast<wxCoord>(rIF.nx()), static_cast<wxCoord>(rIF.ny()), &m_memoryDC, static_cast<wxCoord>(0), static_cast<wxCoord>(0));
+}
+
+
+void
+ImageFileView::OnUpdate(wxView *WXUNUSED(sender), wxObject *WXUNUSED(hint) )
+{
+ const ImageFile& rIF = dynamic_cast<ImageFileDocument*>(GetDocument())->getImageFile();
+ ImageFileArrayConst v = rIF.getArray();
+ int nx = rIF.nx();
+ int ny = rIF.ny();
+ if (v != NULL && nx != 0 && ny != 0) {
+ if (! m_bMinSpecified || ! m_bMaxSpecified) {
+ double min, max;
+ rIF.getMinMax (min, max);
+ if (! m_bMinSpecified)
+ m_minPixel = min;
+ if (! m_bMaxSpecified)
+ m_maxPixel = max;
+ }
+ double scaleWidth = m_maxPixel - m_minPixel;
+
+ m_pBitmap = new wxBitmap (rIF.nx(), rIF.ny());
+ m_memoryDC.SelectObject (*m_pBitmap);
+ m_memoryDC.BeginDrawing();
+
+ wxPen pen (*m_pColours[128], 1, wxDOT);
+ m_memoryDC.SetPen (pen);
+ for (int ix = 0; ix < nx; ix++) {
+ for (int iy = 0; iy < ny; iy++) {
+ double scaleValue = ((v[ix][iy] - m_minPixel) / scaleWidth) * 255;
+ int intensity = static_cast<int>(scaleValue + 0.5);
+ intensity = clamp (intensity, 0, 255);
+ pen.SetColour (*m_pColours[intensity]);
+ m_memoryDC.SetPen(pen);
+ m_memoryDC.DrawPoint(ix, ny - iy);
+ }
+ }
+ }
+ m_memoryDC.EndDrawing();
+
+ if (m_canvas) {
+ // m_canvas->SetScrollbars (50, 50, nx / 50, ny / 50, 0, 0, false);
+ m_canvas->Refresh();
+ }
+
+#ifdef __WXMSW__
+ if (m_canvas)
+ m_canvas->Refresh();
+#else
+ if (m_canvas) {
+ wxClientDC dc(m_canvas);
+ dc.Clear();
+ OnDraw (&dc);
+ }
+#endif
+}
+
+bool
+ImageFileView::OnClose (bool deleteWindow)
+{
+ if (!GetDocument()->Close())
+ return false;
+
+ m_canvas->Clear();
+ m_canvas->m_pView = NULL;
+ m_canvas = NULL;
+ wxString s(wxTheApp->GetAppName());
+ if (m_frame)
+ m_frame->SetTitle(s);
+ SetFrame(NULL);
+
+ Activate(false);
+
+ if (deleteWindow) {
+ delete m_frame;
+ return true;
+ }
+ return true;
+}
+
+
+
+// ProjectionFileCanvas
+
+ProjectionFileCanvas::ProjectionFileCanvas (ProjectionFileView* v, wxFrame *frame, const wxPoint& pos, const wxSize& size, const long style)
+ : wxScrolledWindow(frame, -1, pos, size, style)
+{
+ m_pView = v;
+}
+
+void
+ProjectionFileCanvas::OnDraw(wxDC& dc)
+{
+ if (m_pView)
+ m_pView->OnDraw(& dc);
+}
+
+// ProjectionFileView
+
+IMPLEMENT_DYNAMIC_CLASS(ProjectionFileView, wxView)
+
+BEGIN_EVENT_TABLE(ProjectionFileView, wxView)
+ EVT_MENU(PJMENU_FILE_PROPERTIES, ProjectionFileView::OnProperties)
+ EVT_MENU(PJMENU_FILE_RECONSTRUCT, ProjectionFileView::OnReconstruct)
+END_EVENT_TABLE()
+
+bool ProjectionFileView::m_bPColoursInitialized = false;
+wxColour* ProjectionFileView::m_pColours[256];
+
+ProjectionFileView::ProjectionFileView(void)
+ : wxView(), m_canvas(NULL), m_frame(NULL)
+{
+ if (! m_bPColoursInitialized) {
+ for (int i = 0; i < 256; i++)
+ m_pColours[i] = new wxColour (i, i, i);
+
+ m_bPColoursInitialized = true;
+ }
+}
+
+ProjectionFileView::~ProjectionFileView(void)
+{
+}
+
+void
+ProjectionFileView::OnProperties (wxCommandEvent& event)
+{
+ const Projections& rProj = dynamic_cast<ProjectionFileDocument*>(GetDocument())->getProjections();
+ const string& rFilename = rProj.getFilename();
+ ostringstream os;
+ os << "file: " << rFilename << "\nDetectors: " << rProj.nDet() << "\nViews: " << rProj.nView();
+ wxMessageBox(os.str().c_str(), "Projection Properties", wxOK | wxICON_INFORMATION, m_frame);
+}
+
+
+void
+ProjectionFileView::OnReconstruct (wxCommandEvent& event)
+{
+ const Projections& rProj = dynamic_cast<ProjectionFileDocument*>(GetDocument())->getProjections();
+ const string& rFilename = rProj.getFilename();
+ ostringstream os;
+ os << "Reconstruct file " << rFilename;
+ wxMessageBox(os.str().c_str(), "Reconstruction Dialog", wxOK | wxICON_INFORMATION, m_frame);
+}
+
+
+ProjectionFileCanvas*
+ProjectionFileView::CreateCanvas (wxView *view, wxFrame *parent)
+{
+ ProjectionFileCanvas* pCanvas;
+ int width, height;
+ parent->GetClientSize(&width, &height);
+
+ pCanvas = new ProjectionFileCanvas (dynamic_cast<ProjectionFileView*>(view), parent, wxPoint(0, 0), wxSize(width, height), 0);
+
+ pCanvas->SetScrollbars(20, 20, 50, 50);
+ pCanvas->SetBackgroundColour(*wxWHITE);
+ pCanvas->Clear();
+
+ return pCanvas;
+}
+
+wxFrame*
+ProjectionFileView::CreateChildFrame(wxDocument *doc, wxView *view)
+{
+ wxDocChildFrame *subframe = new wxDocChildFrame(doc, view, dynamic_cast<CTSimApp*>(wxTheApp)->getMainFrame(), -1, "ProjectionFile Frame", wxPoint(10, 10), wxSize(300, 300), wxDEFAULT_FRAME_STYLE);
+
+ wxMenu *file_menu = new wxMenu;
+
+ file_menu->Append(wxID_NEW, "&New...");
+ file_menu->Append(wxID_OPEN, "&Open...");
+ file_menu->Append(wxID_CLOSE, "&Close");
+ file_menu->Append(wxID_SAVE, "&Save");
+ file_menu->Append(wxID_SAVEAS, "Save &As...");
+
+ file_menu->AppendSeparator();
+ file_menu->Append(PJMENU_FILE_RECONSTRUCT, "R&econstuct");
+ file_menu->Append(PJMENU_FILE_PROPERTIES, "P&roperties");
+
+ file_menu->AppendSeparator();
+ file_menu->Append(wxID_PRINT, "&Print...");
+ file_menu->Append(wxID_PRINT_SETUP, "Print &Setup...");
+ file_menu->Append(wxID_PREVIEW, "Print Pre&view");
+
+ wxMenu *help_menu = new wxMenu;
+ help_menu->Append(MAINMENU_HELP_ABOUT, "&About");
+
+ wxMenuBar *menu_bar = new wxMenuBar;
+
+ menu_bar->Append(file_menu, "&File");
+ menu_bar->Append(help_menu, "&Help");
+
+ subframe->SetMenuBar(menu_bar);
+
+ subframe->Centre(wxBOTH);
+
+ return subframe;
+}
+
+
+bool
+ProjectionFileView::OnCreate(wxDocument *doc, long WXUNUSED(flags) )
+{
+ m_frame = CreateChildFrame(doc, this);
+
+ int width, height;
+ m_frame->GetClientSize(&width, &height);
+ m_frame->SetTitle("ProjectionFileView");
+ m_canvas = CreateCanvas(this, m_frame);
+
+#ifdef __X__
+ int x, y; // X requires a forced resize
+ m_frame->GetSize(&x, &y);
+ m_frame->SetSize(-1, -1, x, y);
+#endif
+
+ m_frame->Show(true);
+ Activate(true);
+
+ return true;
+}
+
+void
+ProjectionFileView::OnDraw (wxDC* dc)
+{
+ const Projections& rProj = dynamic_cast<ProjectionFileDocument*>(GetDocument())->getProjections();
+ dc->Blit(static_cast<wxCoord>(0), static_cast<wxCoord>(0), static_cast<wxCoord>(rProj.nDet()), static_cast<wxCoord>(rProj.nView()), &m_memoryDC, static_cast<wxCoord>(0), static_cast<wxCoord>(0));
+}
+
+
+void
+ProjectionFileView::OnUpdate(wxView *WXUNUSED(sender), wxObject *WXUNUSED(hint) )
+{
+ const Projections& rProj = dynamic_cast<ProjectionFileDocument*>(GetDocument())->getProjections();
+ const DetectorArray& detarray = rProj.getDetectorArray (0);
+ const DetectorValue* detval = detarray.detValues();
+ double min = detval[0];
+ double max = detval[0];
+ for (int iy = 0; iy < rProj.nView(); iy++) {
+ const DetectorArray& detarray = rProj.getDetectorArray (iy);
+ detval = detarray.detValues();
+ for (int ix = 0; ix < rProj.nDet(); ix++) {
+ if (min > detval[ix])
+ min = detval[ix];
+ else if (max < detval[ix])
+ max = detval[ix];
+ }
+ }
+
+ double scaleWidth = max - min;
+
+ m_pBitmap = new wxBitmap (rProj.nDet(), rProj.nView());
+ m_memoryDC.SelectObject (*m_pBitmap);
+ m_memoryDC.BeginDrawing();
+
+ wxPen pen (*m_pColours[128], 1, wxDOT);
+ m_memoryDC.SetPen (pen);
+ for (int iy = 0; iy < rProj.nView(); iy++) {
+ const DetectorArray& detarray = rProj.getDetectorArray (iy);
+ detval = detarray.detValues();
+ for (int ix = 0; ix < rProj.nDet(); ix++) {
+ double scaleValue = ((detval[ix] - min) / scaleWidth) * 255;
+ int intensity = static_cast<int>(scaleValue + 0.5);
+ intensity = clamp (intensity, 0, 255);
+ pen.SetColour (*m_pColours[intensity]);
+ m_memoryDC.SetPen(pen);
+ m_memoryDC.DrawPoint(ix, iy);
+ }
+ }
+ m_memoryDC.EndDrawing();
+
+ if (m_canvas) {
+ // m_canvas->SetScrollbars (50, 50, nx / 50, ny / 50, 0, 0, false);
+ m_canvas->Refresh();
+ }
+
+#ifdef __WXMSW__
+ if (m_canvas)
+ m_canvas->Refresh();
+#else
+ if (m_canvas) {
+ wxClientDC dc(m_canvas);
+ dc.Clear();
+ OnDraw (&dc);
+ }
+#endif
+}
+
+bool
+ProjectionFileView::OnClose (bool deleteWindow)
+{
+ if (!GetDocument()->Close())
+ return false;
+
+ m_canvas->Clear();
+ m_canvas->m_pView = NULL;
+ m_canvas = NULL;
+ wxString s(wxTheApp->GetAppName());
+ if (m_frame)
+ m_frame->SetTitle(s);
+ SetFrame(NULL);
+
+ Activate(false);
+
+ if (deleteWindow) {
+ delete m_frame;
+ return true;
+ }
+ return true;
+}
+
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: view.h
+** Purpose: Header file for View & Canvas routines of CTSim program
+** Programmer: Kevin Rosenberg
+** Date Started: July 2000
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: views.h,v 1.1 2000/07/13 07:01:59 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
+******************************************************************************/
+
+#ifdef __GNUG__
+// #pragma interface
+#endif
+
+#ifndef __VIEWSH__
+#define __VIEWSH__
+
+#include "wx/wx.h"
+#include "imagefile.h"
+
+class ImageFileCanvas;
+class ImageFileView : public wxView
+{
+ DECLARE_DYNAMIC_CLASS(ImageFileView)
+
+ wxMemoryDC m_memoryDC;
+ wxBitmap* m_pBitmap;
+
+private:
+ ImageFileCanvas *CreateCanvas(wxView *view, wxFrame *parent);
+ wxFrame *CreateChildFrame(wxDocument *doc, wxView *view);
+
+ ImageFileCanvas *m_canvas;
+ wxFrame *m_frame;
+ bool m_bMinSpecified;
+ bool m_bMaxSpecified;
+ double m_minPixel;
+ double m_maxPixel;
+
+ static bool m_bPColoursInitialized;
+ static wxColour* m_pColours[256];
+
+public:
+ ImageFileView(void);
+ ~ImageFileView(void);
+
+ bool OnCreate(wxDocument *doc, long flags);
+ void OnDraw(wxDC* dc);
+ void OnUpdate(wxView *sender, wxObject *hint = NULL);
+ bool OnClose (bool deleteWindow = true);
+ void OnProperties (wxCommandEvent& event);
+
+ DECLARE_EVENT_TABLE()
+};
+
+class ImageFileCanvas: public wxScrolledWindow
+{
+public:
+ ImageFileView* m_pView;
+
+ ImageFileCanvas (ImageFileView* v, wxFrame *frame, const wxPoint& pos, const wxSize& size, const long style);
+ virtual void OnDraw(wxDC& dc);
+ void OnMouseEvent(wxMouseEvent& event);
+
+ DECLARE_EVENT_TABLE()
+};
+
+
+class MyTextWindow: public wxTextCtrl
+{
+public:
+ wxView *m_pView;
+
+ MyTextWindow(wxView *v, wxFrame *frame, const wxPoint& pos, const wxSize& size, const long style);
+};
+
+
+class TextEditView: public wxView
+{
+ DECLARE_DYNAMIC_CLASS(TextEditView)
+private:
+ wxFrame *CreateChildFrame(wxDocument *doc, wxView *view);
+
+public:
+ wxFrame *frame;
+ MyTextWindow *textsw;
+
+ TextEditView(): wxView() { frame = (wxFrame *) NULL; textsw = (MyTextWindow *) NULL; }
+ ~TextEditView(void) {}
+
+ bool OnCreate(wxDocument *doc, long flags);
+ void OnDraw(wxDC* dc);
+ void OnUpdate(wxView *sender, wxObject *hint = (wxObject *) NULL);
+ bool OnClose(bool deleteWindow = TRUE);
+};
+
+class ProjectionFileCanvas;
+class ProjectionFileView : public wxView
+{
+ DECLARE_DYNAMIC_CLASS(ProjectionFileView)
+
+ wxMemoryDC m_memoryDC;
+ wxBitmap* m_pBitmap;
+
+private:
+ ProjectionFileCanvas *CreateCanvas(wxView *view, wxFrame *parent);
+ wxFrame *CreateChildFrame(wxDocument *doc, wxView *view);
+
+ ProjectionFileCanvas *m_canvas;
+ wxFrame *m_frame;
+
+ static bool m_bPColoursInitialized;
+ static wxColour* m_pColours[256];
+
+public:
+ ProjectionFileView(void);
+ ~ProjectionFileView(void);
+
+ bool OnCreate(wxDocument *doc, long flags);
+ void OnDraw(wxDC* dc);
+ void OnUpdate(wxView *sender, wxObject *hint = NULL);
+ bool OnClose (bool deleteWindow = true);
+ void OnProperties (wxCommandEvent& event);
+ void OnReconstruct (wxCommandEvent& event);
+
+ DECLARE_EVENT_TABLE()
+};
+
+class ProjectionFileCanvas: public wxScrolledWindow
+{
+public:
+ ProjectionFileView* m_pView;
+
+ ProjectionFileCanvas (ProjectionFileView* v, wxFrame *frame, const wxPoint& pos, const wxSize& size, const long style);
+ virtual void OnDraw(wxDC& dc);
+};
+
+
+#endif
--- /dev/null
+bin_PROGRAMS = pjrec phm2pj pj2if @lamprograms@ ifinfo phm2if if-1 if-2 if2img
+bin_SCRIPTS = sample-ctsim.sh
+EXTRA_PROGRAMS = pjrec-lam phm2pj-lam phm2if-lam
+INCLUDES=@my_includes@
+EXTRA_DIST=Makefile.nt mpiworld.cpp
+SOURCE_DEPEND=../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a ../libctgraphics/libctgraphics.a
+
+pjrec_SOURCES = pjrec.cpp
+pjrec_LDADD=@ctlibs@
+pjrec_DEPENDENCIES=$(SOURCE_DEPEND)
+phm2pj_SOURCES=phm2pj.cpp
+phm2pj_LDADD=@ctlibs@
+phm2pj_DEPENDENCIES=$(SOURCE_DEPEND)
+phm2if_SOURCES = phm2if.cpp
+phm2if_LDADD=@ctlibs@
+phm2if_DEPENDENCIES=$(SOURCE_DEPEND)
+if2img_SOURCES = if2img.cpp
+if2img_LDADD=@ctlibs@
+if2img_DEPENDENCIES=$(SOURCE_DEPEND)
+pj2if_SOURCES = pj2if.cpp
+pj2if_LDADD=@ctlibs@
+pj2if_DEPENDENCIES=$(SOURCE_DEPEND)
+if_1_SOURCES=if-1.cpp
+if_1_LDADD=@ctlibs@
+if_1_DEPENDENCIES=$(SOURCE_DEPEND)
+if_2_SOURCES=if-2.cpp
+if_2_LDADD=@ctlibs@
+if_2_DEPENDENCIES=$(SOURCE_DEPEND)
+ifinfo_SOURCES = ifinfo.cpp
+ifinfo_LDADD=@ctlibs@
+ifinfo_DEPENDENCIES=$(SOURCE_DEPEND)
+
+pjrec_lam_SOURCES=pjrec.cpp
+pjrec_lam_LDADD=@ctlamlibs@
+phm2if_lam_SOURCES=phm2if.cpp
+phm2if_lam_LDADD=@ctlamlibs@
+phm2pj_lam_SOURCES=phm2pj.cpp
+phm2pj_lam_LDADD=@ctlamlibs@
+
+realclean:
+ rm -f *.pgm *.if *~ *.pj
+
+if USE_LAM
+CC_LAM = $(lamdir)/bin/balky
+LAM_EXTRA_SRC = mpiworld.cpp
+
+pjrec-lam: pjrec.cpp mpiworld.cpp ../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a
+ $(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI pjrec.cpp -o pjrec-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@
+
+phm2pj-lam: phm2pj.cpp mpiworld.cpp ../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a
+ $(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI phm2pj.cpp -o phm2pj-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@
+
+phm2if-lam: phm2if.cpp mpiworld.cpp ../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a
+ $(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI phm2if.cpp -o phm2if-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@
+
+endif
+
+shared: pjrec.cpp phm2pj.cpp
+ $(CC) @DEFS@ $(CFLAGS) $(INCLUDES) -DNO_MAIN -shared phm2pj.cpp pjrec.cpp @ctlibs@ $(LDFLAGS) -o ir.so
+
+
+
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: if-1.cpp
+** Purpose: Manipulate a single image file
+** Programmer: Kevin Rosenberg
+** Date Started: Aug 1984
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: if-1.cpp,v 1.1 2000/07/13 07:01:35 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
+******************************************************************************/
+
+/* FILE
+ * if-1.c Filter a single IF file
+ */
+
+#include "ct.h"
+
+enum {O_LOG, O_EXP, O_SQRT, O_SQR, O_INVERT, O_VERBOSE, O_HELP, O_VERSION};
+
+static struct option my_options[] =
+{
+ {"invert", 0, 0, O_INVERT},
+ {"verbose", 0, 0, O_VERBOSE},
+ {"log", 0, 0, O_LOG},
+ {"exp", 0, 0, O_EXP},
+ {"sqr", 0, 0, O_SQR},
+ {"sqrt", 0, 0, O_SQRT},
+ {"help", 0, 0, O_HELP},
+ {"version", 0, 0, O_VERSION},
+ {0, 0, 0, 0}
+};
+
+void
+if1_usage (const char *program)
+{
+ cout << "usage: " << fileBasename(program) << " infile outfile [OPTIONS]" << endl;
+ cout << "Generate a IF file from a IF file" << endl;
+ cout << endl;
+ cout << " --invert Invert image" << endl;
+ cout << " --log Natural logrithm of image" << endl;
+ cout << " --exp Natural exponential of image" << endl;
+ cout << " --sqr Square of image" << endl;
+ cout << " --sqrt Square root of image" << endl;
+ cout << " --verbose Verbose modem" << endl;
+ cout << " --version Print version" << endl;
+ cout << " --help Print this help message" << endl;
+}
+
+int
+if1_main (int argc, char *const argv[])
+{
+ ImageFile *im_in;
+ ImageFile *im_out;
+ char *in_file;
+ char *out_file;
+ int opt_verbose = 0;
+ int opt_invert = 0;
+ int opt_log = 0;
+ int opt_exp = 0;
+ int opt_sqr = 0;
+ int opt_sqrt = 0;
+
+ while (1)
+ {
+ int c = getopt_long (argc, argv, "", my_options, NULL);
+
+ if (c == -1)
+ break;
+
+ switch (c)
+ {
+ case O_INVERT:
+ opt_invert = 1;
+ break;
+ case O_LOG:
+ opt_log = 1;
+ break;
+ case O_SQR:
+ opt_sqr = 1;
+ break;
+ case O_SQRT:
+ opt_sqrt = 1;
+ break;
+ case O_EXP:
+ opt_exp = 1;
+ break;
+ case O_VERBOSE:
+ opt_verbose = 1;
+ break;
+ case O_VERSION:
+#ifdef VERSION
+ cout << "Version " << VERSION << endl;
+#else
+ cout << "Unknown version number" << endl;
+#endif
+ return (0);
+ case O_HELP:
+ case '?':
+ if1_usage(argv[0]);
+ return (0);
+ default:
+ if1_usage(argv[0]);
+ return (1);
+ }
+ }
+
+ if (optind + 2 != argc)
+ {
+ if1_usage(argv[0]);
+ return (1);
+ }
+
+ in_file = argv[optind];
+ out_file = argv[optind + 1];
+
+
+ string histString;
+
+ if (opt_invert || opt_log || opt_exp || opt_sqr || opt_sqrt) {
+ int ix, iy;
+
+ im_in = new ImageFile ();
+ im_in->fileRead (in_file);
+ int nx = im_in->nx();
+ int ny = im_in->ny();
+ im_out = new ImageFile (nx, ny);
+
+ ImageFileArray vIn = im_in->getArray();
+ ImageFileArray vOut = im_out->getArray();
+
+ if (opt_invert) {
+ for (ix = 0; ix < nx; ix++)
+ for (iy = 0; iy < ny; iy++)
+ vOut[ix][iy] = - vIn[ix][iy];
+ histString = "Invert transformation";
+ }
+ if (opt_log) {
+ for (ix = 0; ix < nx; ix++)
+ for (iy = 0; iy < ny; iy++)
+ vOut[ix][iy] = log (vIn[ix][iy]);
+ histString = "Log transformation";
+ }
+ if (opt_exp) {
+ for (ix = 0; ix < nx; ix++)
+ for (iy = 0; iy < ny; iy++)
+ vOut[ix][iy] = exp (vIn[ix][iy]);
+ histString = "Exp transformation";
+ }
+ if (opt_sqr) {
+ for (ix = 0; ix < nx; ix++)
+ for (iy = 0; iy < ny; iy++)
+ vOut[ix][iy] = vIn[ix][iy] * vIn[ix][iy];
+ histString = "Sqr transformation";
+ }
+ if (opt_sqrt) {
+ for (ix = 0; ix < nx; ix++)
+ for (iy = 0; iy < ny; iy++)
+ vOut[ix][iy] = sqrt (vIn[ix][iy]);
+ histString = "Sqrt transformation";
+ }
+
+ im_out->labelsCopy (*im_in);
+ im_out->labelAdd (Array2dFileLabel::L_HISTORY, histString.c_str());
+ im_out->fileWrite (out_file);
+ }
+
+ return (0);
+}
+
+#ifndef NO_MAIN
+int
+main (int argc, char *const argv[])
+{
+ int retval = 1;
+
+ try {
+ retval = if1_main(argc, argv);
+ } catch (exception e) {
+ cerr << "Exception: " << e.what() << endl;
+ } catch (...) {
+ cerr << "Unknown exception" << endl;
+ }
+
+ return (retval);
+}
+#endif
+
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: if-2.cpp
+** Purpose: Manipulate two image files
+** Programmer: Kevin Rosenberg
+** Date Started: May 2000
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: if-2.cpp,v 1.1 2000/07/13 07:01:35 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_ADD, O_SUB, O_MUL, O_COMP, O_ROW_PLOT, O_COLUMN_PLOT, O_VERBOSE, O_HELP, O_VERSION};
+
+static struct option my_options[] =
+{
+ {"add", 0, 0, O_ADD},
+ {"sub", 0, 0, O_SUB},
+ {"mul", 0, 0, O_MUL},
+ {"comp", 0, 0, O_COMP},
+ {"column-plot", 1, 0, O_COLUMN_PLOT},
+ {"row-plot", 1, 0, O_ROW_PLOT},
+ {"verbose", 0, 0, O_VERBOSE},
+ {"help", 0, 0, O_HELP},
+ {"version", 0, 0, O_VERSION},
+ {0, 0, 0, 0}
+};
+
+void
+if2_usage (const char *program)
+{
+ cout << "usage: " << fileBasename(program) << " infile1 infile2 outfile [OPTIONS]" << endl;
+ cout << "Perform functions on two input image files" << endl;
+ cout << endl;
+ cout << " infile1 Name of first input IF file" << endl;
+ cout << " infile2 Name of second input IF file" << endl;
+ cout << " outfile Name of output IF file" << endl;
+ cout << " --add Add images" << endl;
+ cout << " --sub Subtract image 2 from image 1" << endl;
+ cout << " --mul Multiply images" << endl;
+ cout << " --comp Compare images" << endl;
+ cout << " --column-plot n Plot column\n";
+ cout << " --row-plot n Plot row\n";
+ cout << " --verbose Verbose modem" << endl;
+ cout << " --version Print version" << endl;
+ cout << " --help Print this help message" << endl;
+}
+
+int
+if2_main (int argc, char *const argv[])
+{
+ ImageFile* pim_in1;
+ ImageFile* pim_in2;
+ ImageFile* pim_out = NULL;
+ string in_file1;
+ string in_file2;
+ string out_file;
+ int opt_verbose = 0;
+ int opt_add = 0;
+ int opt_sub = 0;
+ int opt_mul = 0;
+ int opt_comp = 0;
+ int opt_outputFile = 0;
+ int opt_rowPlot = -1;
+ int opt_columnPlot = -1;
+ Timer timerProgram;
+
+ while (1) {
+ char* endptr;
+ int c = getopt_long (argc, argv, "", my_options, NULL);
+
+ if (c == -1)
+ break;
+
+ switch (c) {
+ case O_ADD:
+ opt_add = 1;
+ opt_outputFile = 1;
+ break;
+ case O_SUB :
+ opt_sub = 1;
+ opt_outputFile = 1;
+ break;
+ case O_MUL:
+ opt_mul = 1;
+ opt_outputFile = 1;
+ break;
+ case O_ROW_PLOT:
+ opt_rowPlot = strtol(optarg, &endptr, 10);
+ if (endptr != optarg + strlen(optarg)) {
+ if2_usage(argv[0]);
+ }
+ case O_COLUMN_PLOT:
+ opt_columnPlot = strtol(optarg, &endptr, 10);
+ if (endptr != optarg + strlen(optarg)) {
+ if2_usage(argv[0]);
+ }
+ case O_COMP:
+ opt_comp = 1;
+ break;
+ case O_VERBOSE:
+ opt_verbose = 1;
+ break;
+ case O_VERSION:
+#ifdef VERSION
+ cout << "Version " << VERSION << endl;
+#else
+ cout << "Unknown version number" << endl;
+#endif
+ return (0);
+ case O_HELP:
+ case '?':
+ if2_usage(argv[0]);
+ return (0);
+ default:
+ if2_usage(argv[0]);
+ return (1);
+ }
+ }
+
+ if (opt_outputFile && (optind + 3 != argc)) {
+ if2_usage(argv[0]);
+ return (1);
+ }
+ else if (! opt_outputFile && optind + 2 != argc) {
+ if2_usage(argv[0]);
+ return (1);
+ }
+
+ in_file1 = argv[optind];
+ in_file2 = argv[optind + 1];
+ if (opt_outputFile)
+ out_file = argv[optind + 2];
+
+ pim_in1 = new ImageFile ();
+ pim_in2 = new ImageFile ();
+ ImageFile& im_in1 = *pim_in1;
+ ImageFile& im_in2 = *pim_in2;
+
+ if (! im_in1.fileRead(in_file1) || ! im_in2.fileRead(in_file2)) {
+ sys_error (ERR_WARNING, "Error reading an image");
+ return (1);
+ }
+
+ if (im_in1.nx() != im_in2.nx() || im_in1.ny() != im_in2.ny()) {
+ sys_error (ERR_SEVERE, "Error: Size of image 1 (%d,%d) and image 2 (%d,%d) do not match",
+ im_in1.nx(), im_in1.ny(), im_in2.nx(), im_in2.ny());
+ return(1);
+ }
+ if (im_in1.nx() < 0 || im_in1.ny() < 0) {
+ sys_error (ERR_SEVERE, "Error: Size of image < 0");
+ return(1);
+ }
+
+ ImageFileArray v1 = im_in1.getArray();
+ ImageFileArray v2 = im_in2.getArray();
+ ImageFileArray vout = NULL;
+
+ if (opt_outputFile) {
+ pim_out = new ImageFile (im_in1.nx(), im_in1.ny());
+ vout = pim_out->getArray();
+ }
+
+ string strOperation;
+ int nx = im_in1.nx();
+ int ny = im_in1.ny();
+ int nx2 = im_in2.nx();
+ int ny2 = im_in2.ny();
+
+ if (opt_add) {
+ strOperation = "Add Images";
+ for (int ix = 0; ix < nx; ix++) {
+ ImageFileColumn in1 = v1[ix];
+ ImageFileColumn in2 = v2[ix];
+ ImageFileColumn out = vout[ix];
+ for (int iy = 0; iy < ny; iy++)
+ *out++ = *in1++ + *in2++;
+ }
+ } else if (opt_sub) {
+ strOperation = "Subtract Images";
+ for (int ix = 0; ix < nx; ix++) {
+ ImageFileColumn in1 = v1[ix];
+ ImageFileColumn in2 = v2[ix];
+ ImageFileColumn out = vout[ix];
+ for (int iy = 0; iy < ny; iy++)
+ *out++ = *in1++ - *in2++;
+ }
+ } else if (opt_mul) {
+ strOperation = "Multiply Images";
+ for (int ix = 0; ix < nx; ix++) {
+ ImageFileColumn in1 = v1[ix];
+ ImageFileColumn in2 = v2[ix];
+ ImageFileColumn out = vout[ix];
+ for (int iy = 0; iy < ny; iy++)
+ *out++ = *in1++ * *in2++;
+ }
+ }
+ if (opt_comp) {
+ double d, r, e;
+ im_in1.comparativeStatistics (im_in2, d, r, e);
+ cout << "d=" << d << ", r=" << r << ", e=" << e << endl;
+ }
+ if (opt_columnPlot > 0) {
+ if (opt_columnPlot >= nx || opt_columnPlot >= nx2) {
+ sys_error (ERR_SEVERE, "column-plot > nx");
+ return (1);
+ }
+ double plot_xaxis [nx];
+ for (int i = 0; i < nx; i++)
+ plot_xaxis[i] = i;
+#if HAVE_SGP
+#if 0
+#else
+ ezset ("clear.");
+ ezset ("xticks major 5.");
+ ezset ("xlabel Column");
+ ezset ("ylabel Pixel");
+ ezset ("curves 2");
+ ezset ("box.");
+ ezset ("grid.");
+ ezplot (v1[opt_columnPlot], plot_xaxis, im_in1.ny());
+ ezplot (v2[opt_columnPlot], plot_xaxis, im_in2.ny());
+#endif
+ char str[256];
+ cout << "Press enter to continue" << flush;
+ fgets(str, sizeof(str), stdin);
+ sgp2_close (sgp2_get_active_win());
+#endif
+ }
+
+ if (opt_rowPlot > 0) {
+ if (opt_rowPlot >= ny || opt_rowPlot >= ny2) {
+ sys_error (ERR_SEVERE, "row_plot > ny");
+ return (1);
+ }
+ double plot_xaxis [ny];
+ double v1Row[nx], v2Row[nx2];
+
+ for (int i = 0; i < ny; i++)
+ plot_xaxis[i] = i;
+ for (int i = 0; i < nx; i++)
+ v1Row[i] = v1[opt_rowPlot][i];
+ for (int i = 0; i < nx2; i++)
+ v2Row[i] = v2[opt_rowPlot][i];
+
+#if HAVE_SGP
+#if 0
+#else
+ ezset ("clear.");
+ ezset ("xticks major 5.");
+ ezset ("xlabel Column");
+ ezset ("ylabel Pixel");
+ ezset ("curves 2");
+ ezset ("box.");
+ ezset ("grid.");
+ ezplot (v1Row, plot_xaxis, im_in1.nx());
+ ezplot (v2Row, plot_xaxis, im_in2.nx());
+#endif
+ char str[256];
+ cout << "Press enter to continue" << flush;
+ fgets(str, sizeof(str), stdin);
+ sgp2_close (sgp2_get_active_win());
+#endif
+ }
+
+ if (opt_outputFile) {
+ pim_out->labelsCopy (im_in1, "if-2 file 1: ");
+ pim_out->labelsCopy (im_in2, "if-2 file 2: ");
+ pim_out->labelAdd (Array2dFileLabel::L_HISTORY, strOperation.c_str(), timerProgram.timerEnd());
+ pim_out->fileWrite(out_file);
+ }
+
+ return (0);
+}
+
+#ifndef NO_MAIN
+int
+main (int argc, char *const argv[])
+{
+ int retval = 1;
+
+ try {
+ retval = if2_main(argc, argv);
+ } catch (exception e) {
+ cerr << "Exception: " << e.what() << endl;
+ } catch (...) {
+ cerr << "Unknown exception" << endl;
+ }
+
+ return (retval);
+}
+#endif
+
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: if2img.cpp
+** Purpose: Convert an image file to a viewable image
+** Programmer: Kevin Rosenberg
+** Date Started: April 2000
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: if2img.cpp,v 1.1 2000/07/13 07:01:35 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"
+
+
+enum { O_SCALE, O_MIN, O_MAX, O_AUTO, O_CENTER, O_STATS, O_FORMAT, O_LABELS,
+ O_HELP, O_VERBOSE, O_VERSION, O_DEBUG };
+
+static struct option my_options[] =
+{
+ {"scale", 1, 0, O_SCALE},
+ {"min", 1, 0, O_MIN},
+ {"max", 1, 0, O_MAX},
+ {"auto", 1, 0, O_AUTO},
+ {"center", 1, 0, O_CENTER},
+ {"format", 1, 0, O_FORMAT},
+ {"labels", 0, 0, O_LABELS},
+ {"stats", 0, 0, O_STATS},
+ {"verbose", 0, 0, O_VERBOSE},
+ {"debug", 0, 0, O_DEBUG},
+ {"help", 0, 0, O_HELP},
+ {"version", 0, 0, O_VERSION},
+ {0, 0, 0, 0}
+};
+
+enum { O_AUTO_FULL, O_AUTO_STD0_1, O_AUTO_STD0_5, O_AUTO_STD1, O_AUTO_STD2, O_AUTO_STD3 };
+static const char O_AUTO_FULL_STR[]="full";
+static const char O_AUTO_STD0_1_STR[]="std0.1";
+static const char O_AUTO_STD0_5_STR[]="std0.5";
+static const char O_AUTO_STD1_STR[]="std1";
+static const char O_AUTO_STD2_STR[]="std2";
+static const char O_AUTO_STD3_STR[]="std3";
+
+enum { O_CENTER_MEAN, O_CENTER_MODE };
+static const char O_CENTER_MEAN_STR[]="mean";
+static const char O_CENTER_MODE_STR[]="mode";
+
+enum { O_FORMAT_GIF, O_FORMAT_PNG, O_FORMAT_PNG16, O_FORMAT_PGM, O_FORMAT_PGMASC, O_FORMAT_DISP };
+static const char O_FORMAT_GIF_STR[]="gif";
+static const char O_FORMAT_PNG_STR[]="png" ;
+static const char O_FORMAT_PNG16_STR[]="png16" ;
+static const char O_FORMAT_PGM_STR[]="pgm";
+static const char O_FORMAT_PGMASC_STR[]="pgmasc";
+static const char O_FORMAT_DISP_STR[]="disp";
+
+void
+if2img_usage (const char *program)
+{
+ cout << "usage: " << fileBasename(program) << " ifname outfile [OPTIONS]" << endl;
+ cout << "Convert IF file to an image file" << endl;
+ cout << endl;
+ cout << " ifname Name of input file" << endl;
+ cout << " outfile Name of output file" << endl;
+ cout << " --format Output format" << endl;
+ cout << " pgm PGM (portable graymap) format (default)" << endl;
+ cout << " pgmasc PGM (portable graymap) ASCII format" << endl;
+#ifdef HAVE_PNG
+ cout << " png PNG (8-bit) format" << endl;
+ cout << " png16 PNG (16-bit) format" << endl;
+#endif
+#if HAVE_G2
+ cout << " gif GIF format" << endl;
+#endif
+ cout << " disp Display on screen" << endl;
+ cout << " --center Center of window" << endl;
+ cout << " mode Mode is center of window (default)" << endl;
+ cout << " mean Mean is center of window" << endl;
+ cout << " --auto Set auto window" << endl;
+ cout << " full Use full window (default)" << endl;
+ cout << " std0.1 Use 0.1 standard deviation about center" << endl;
+ cout << " std0.5 Use 0.5 standard deviation about center" << endl;
+ cout << " std1 Use one standard deviation about center" << endl;
+ cout << " std2 Use two standard deviations about center" << endl;
+ cout << " std3 Use three standard deviations about center" << endl;
+ cout << " --scale Scaling factor for output size" << endl;
+ cout << " --min Set minimum intensity" << endl;
+ cout << " --max Set maximum intensity" << endl;
+ cout << " --stats Print image statistics" << endl;
+ cout << " --labels Print image labels" << endl;
+ cout << " --debug Set debug mode" << endl;
+ cout << " --verbose Set verbose mode" << endl;
+ cout << " --version Print version" << endl;
+ cout << " --help Print this help message" << endl;
+}
+
+
+int
+if2img_main (int argc, char *const argv[])
+{
+ ImageFile* pim = NULL;
+ double densmin = HUGE_VAL, densmax = -HUGE_VAL;
+ char *in_file, *out_file;
+ int opt_verbose = 0;
+ int opt_scale = 1;
+ int opt_auto = O_AUTO_FULL;
+ int opt_set_max = 0;
+ int opt_set_min = 0;
+ int opt_stats = 0;
+ int opt_debug = 0;
+ int opt_center = O_CENTER_MODE;
+ int opt_format = O_FORMAT_PGM;
+ int opt_labels = 0;
+
+ while (1)
+ {
+ int c = getopt_long (argc, argv, "", my_options, NULL);
+ char *endptr = NULL;
+ char *endstr;
+
+ if (c == -1)
+ break;
+
+ switch (c)
+ {
+ case O_MIN:
+ opt_set_min = 1;
+ densmin = strtod(optarg, &endptr);
+ endstr = optarg + strlen(optarg);
+ if (endptr != endstr)
+ {
+ sys_error (ERR_SEVERE, "Error setting --min to %s", optarg);
+ if2img_usage(argv[0]);
+ return (1);
+ }
+ break;
+ case O_MAX:
+ opt_set_max = 1;
+ densmax = strtod(optarg, &endptr);
+ endstr = optarg + strlen(optarg);
+ if (endptr != endstr)
+ {
+ sys_error (ERR_SEVERE, "Error setting --max to %s", optarg);
+ if2img_usage(argv[0]);
+ return (1);
+ }
+ break;
+ case O_SCALE:
+ opt_scale = strtol(optarg, &endptr, 10);
+ endstr = optarg + strlen(optarg);
+ if (endptr != endstr)
+ {
+ sys_error (ERR_SEVERE, "Error setting --scale to %s", optarg);
+ if2img_usage(argv[0]);
+ return (1);
+ }
+ break;
+ case O_AUTO:
+ if (strcmp(optarg, O_AUTO_FULL_STR) == 0)
+ opt_auto = O_AUTO_FULL;
+ else if (strcmp(optarg, O_AUTO_STD1_STR) == 0)
+ opt_auto = O_AUTO_STD1;
+ else if (strcmp(optarg, O_AUTO_STD0_5_STR) == 0)
+ opt_auto = O_AUTO_STD0_5;
+ else if (strcmp(optarg, O_AUTO_STD0_1_STR) == 0)
+ opt_auto = O_AUTO_STD0_1;
+ else if (strcmp(optarg, O_AUTO_STD2_STR) == 0)
+ opt_auto = O_AUTO_STD2;
+ else if (strcmp(optarg, O_AUTO_STD3_STR) == 0)
+ opt_auto = O_AUTO_STD3;
+ else
+ {
+ sys_error (ERR_SEVERE, "Invalid auto mode %s", optarg);
+ if2img_usage(argv[0]);
+ return (1);
+ }
+ break;
+ case O_CENTER:
+ if (strcmp(optarg, O_CENTER_MEAN_STR) == 0)
+ opt_center = O_CENTER_MEAN;
+ else if (strcmp(optarg, O_CENTER_MODE_STR) == 0)
+ opt_center = O_CENTER_MODE;
+ else
+ {
+ sys_error (ERR_SEVERE, "Invalid center mode %s", optarg);
+ if2img_usage(argv[0]);
+ return (1);
+ }
+ break;
+ case O_FORMAT:
+ if (strcmp(optarg, O_FORMAT_PGM_STR) == 0)
+ opt_format = O_FORMAT_PGM;
+ else if (strcmp(optarg, O_FORMAT_PGMASC_STR) == 0)
+ opt_format = O_FORMAT_PGMASC;
+#if HAVE_PNG
+ else if (strcmp(optarg, O_FORMAT_PNG_STR) == 0)
+ opt_format = O_FORMAT_PNG;
+ else if (strcmp(optarg, O_FORMAT_PNG16_STR) == 0)
+ opt_format = O_FORMAT_PNG16;
+#endif
+#if HAVE_GIF
+ else if (strcmp(optarg, O_FORMAT_GIF_STR) == 0)
+ opt_format = O_FORMAT_GIF;
+#endif
+ else if (strcmp(optarg, O_FORMAT_DISP_STR) == 0)
+ opt_format = O_FORMAT_DISP;
+ else {
+ sys_error (ERR_SEVERE, "Invalid format mode %s", optarg);
+ if2img_usage(argv[0]);
+ return (1);
+ }
+ break;
+ case O_LABELS:
+ opt_labels = 1;
+ break;
+ case O_VERBOSE:
+ opt_verbose = 1;
+ break;
+ case O_DEBUG:
+ opt_debug = 1;
+ break;
+ case O_STATS:
+ opt_stats = 1;
+ break;
+ case O_VERSION:
+#ifdef VERSION
+ cout << "Version " << VERSION;
+#else
+ cout << "Unknown version number" << endl;
+#endif
+ return (0);
+ case O_HELP:
+ case '?':
+ if2img_usage(argv[0]);
+ return (0);
+ default:
+ if2img_usage(argv[0]);
+ return (1);
+ }
+ }
+
+ if ((opt_format == O_FORMAT_DISP && optind + 1 != argc)
+ || (opt_format != O_FORMAT_DISP && optind + 2 != argc))
+ {
+ if2img_usage(argv[0]);
+ return (1);
+ }
+
+ in_file = argv[optind];
+
+ if (opt_format != O_FORMAT_DISP)
+ out_file = argv[optind+1];
+ else out_file = NULL;
+
+ pim = new ImageFile ();
+ ImageFile& im = *pim;
+ if (! im.fileRead(in_file)) {
+ sys_error (ERR_FATAL, "File %s does not exist", in_file);
+ return (1);
+ }
+
+ if (opt_labels)
+ im.printLabels(cout);
+
+ if (opt_stats || (! (opt_set_max && opt_set_min))) {
+ double min, max, mean, mode, median, stddev;
+ double window = 0;
+ im.statistics(min, max, mean, mode, median, stddev);
+
+ if (opt_auto == O_AUTO_FULL) {
+ if (! opt_set_max)
+ densmax = max;
+ if (! opt_set_min)
+ densmin = min;
+ }
+ if (opt_stats || opt_auto != O_AUTO_FULL) {
+
+ if (opt_auto == O_AUTO_FULL)
+ ;
+ else if (opt_auto == O_AUTO_STD1)
+ window = stddev;
+ else if (opt_auto == O_AUTO_STD0_1)
+ window = stddev * 0.1;
+ else if (opt_auto == O_AUTO_STD0_5)
+ window = stddev * 0.5;
+ else if (opt_auto == O_AUTO_STD2)
+ window = stddev * 2;
+ else if (opt_auto == O_AUTO_STD3)
+ window = stddev * 3;
+ else {
+ sys_error (ERR_SEVERE, "Internal Error: Invalid auto mode %d", opt_auto);
+ return (1);
+ }
+ }
+ if (opt_stats) {
+ cout <<"nx: " << im.nx() << endl;
+ cout <<"ny: " << im.ny() << endl;
+ cout <<"min: " << min << endl;
+ cout <<"max: " << max << endl;
+ cout <<"mean: " << mean << endl;
+ cout <<"mode: " << mode << endl;
+ cout <<"stddev: " << stddev << endl;
+ }
+ if (opt_auto != O_AUTO_FULL) {
+ double center;
+
+ if (opt_center == O_CENTER_MODE)
+ center = mode;
+ else if (opt_center == O_CENTER_MEAN)
+ center = mean;
+ else {
+ sys_error (ERR_SEVERE, "Internal Error: Invalid center mode %d", opt_center);
+ return (1);
+ }
+ if (! opt_set_max)
+ densmax = center + window;
+ if (! opt_set_min)
+ densmin = center - window;
+ }
+ }
+
+ if (opt_stats) {
+ cout << "min display: " << densmin << endl;
+ cout << "max display: " << densmax << endl;
+ }
+
+ if (opt_format == O_FORMAT_PGM)
+ im.writeImagePGM (out_file, opt_scale, opt_scale, densmin, densmax);
+ else if (opt_format == O_FORMAT_PGMASC)
+ im.writeImagePGMASCII (out_file, opt_scale, opt_scale, densmin, densmax);
+#if HAVE_PNG
+ else if (opt_format == O_FORMAT_PNG)
+ im.writeImagePNG (out_file, 8, opt_scale, opt_scale, densmin, densmax);
+ else if (opt_format == O_FORMAT_PNG16)
+ im.writeImagePNG (out_file, 16, opt_scale, opt_scale, densmin, densmax);
+#endif
+#if HAVE_GIF
+ else if (opt_format == O_FORMAT_GIF)
+ im.writeImageGIF (out_file, opt_scale, opt_scale, densmin, densmax);
+#endif
+ else if (opt_format == O_FORMAT_DISP) {
+#if HAVE_SGP
+ im.displayScaling (opt_scale, densmin, densmax);
+ cout << "Press enter to continue\n";
+ cio_kb_getc();
+ sgp2_close(sgp2_get_active_win());
+#endif
+ }
+ else
+ {
+ sys_error (ERR_SEVERE, "Internal Error: Invalid format mode %d", opt_format);
+ return (1);
+ }
+ return (0);
+}
+
+
+#ifndef NO_MAIN
+int
+main (int argc, char *const argv[])
+{
+ int retval = 1;
+
+ try {
+ retval = if2img_main(argc, argv);
+ } catch (exception e) {
+ cerr << "Exception: " << e.what() << endl;
+ } catch (...) {
+ cerr << "Unknown exception" << endl;
+ }
+
+ return (retval);
+}
+#endif
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: ifinfo.cpp
+** Purpose: Display information about an image file
+** Programmer: Kevin Rosenberg
+** Date Started: April 2000
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: ifinfo.cpp,v 1.1 2000/07/13 07:01:35 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
+******************************************************************************/
+
+/* FILE
+ * ifinfo.c Display info on sdf files
+ */
+
+#include "ct.h"
+
+enum { O_LABELS, O_STATS, O_NO_STATS, O_NO_LABELS, O_VERBOSE, O_HELP, O_VERSION, O_DEBUG };
+
+static struct option my_options[] =
+{
+ {"labels", 0, 0, O_LABELS},
+ {"no-labels", 0, 0, O_NO_LABELS},
+ {"stats", 0, 0, O_STATS},
+ {"no-stats", 0, 0, O_NO_STATS},
+ {"debug", 0, 0, O_DEBUG},
+ {"verbose", 0, 0, O_VERBOSE},
+ {"help", 0, 0, O_HELP},
+ {"version", 0, 0, O_VERSION},
+ {0, 0, 0, 0}
+};
+
+void
+ifinfo_usage (const char *program)
+{
+ cout << "usage: " << fileBasename(program) << " image-filename [OPTIONS]" << endl;
+ cout << "Imagefile information" << endl;
+ cout << endl;
+ cout << " infile Name of input IF file" << endl;
+ cout << " --display Display image" << endl;
+ cout << " --labels Print image labels (default)" << endl;
+ cout << " --no-labels Do not print image labels" << endl;
+ cout << " --stats Print image statistics (default)" << endl;
+ cout << " --no-stats Do not print image statistics" << endl;
+ cout << " --debug Debug mode" << endl;
+ cout << " --verbose Verbose mode" << endl;
+ cout << " --version Print version" << endl;
+ cout << " --help Print this help message" << endl;
+}
+
+int
+ifinfo_main (int argc, char *const argv[])
+{
+ ImageFile *im = NULL;
+ string in_file;
+ int opt_verbose = 0;
+ int opt_stats = 1;
+ int opt_labels = 1;
+ int opt_debug = 0;
+
+ while (1)
+ {
+ int c = getopt_long (argc, argv, "", my_options, NULL);
+
+ if (c == -1)
+ break;
+
+ switch (c)
+ {
+ case O_LABELS:
+ opt_labels = 1;
+ break;
+ case O_STATS:
+ opt_stats = 1;
+ break;
+ case O_NO_LABELS:
+ opt_labels = 0;
+ break;
+ case O_NO_STATS:
+ opt_stats = 0;
+ break;
+ case O_VERBOSE:
+ opt_verbose = 1;
+ break;
+ case O_DEBUG:
+ opt_debug = 0;
+ break;
+ case O_VERSION:
+#ifdef VERSION
+ cout << "Version " << VERSION << endl;
+#else
+ cout << "Unknown version number" << endl;
+#endif
+ return (0);
+ case O_HELP:
+ case '?':
+ ifinfo_usage(argv[0]);
+ return (0);
+ default:
+ ifinfo_usage(argv[0]);
+ return (1);
+ }
+ }
+
+ if (optind + 1 != argc) {
+ ifinfo_usage (argv[0]);
+ return (1);
+ }
+
+ in_file = argv[optind];
+
+ im = new ImageFile ();
+ if (! im->fileRead (in_file)) {
+ sys_error (ERR_WARNING, "Unable to read file %s", in_file.c_str());
+ return (1);
+ }
+
+ if (opt_labels)
+ im->printLabels (cout);
+
+ if (opt_stats) {
+ cout << "Size: (" << im->nx() << "," << im->ny() << ")" << endl;
+ im->printStatistics (cout);
+ }
+
+ return (0);
+}
+
+#ifndef NO_MAIN
+int
+main (int argc, char *const argv[])
+{
+ int retval = 1;
+
+ try {
+ retval = ifinfo_main(argc, argv);
+ } catch (exception e) {
+ cerr << "Exception: " << e.what() << endl;
+ } catch (...) {
+ cerr << "Unknown exception" << endl;
+ }
+
+ return (retval);
+}
+#endif
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: mpiworld.cpp
+** Purpose: MPI Support class
+** Programmer: Kevin Rosenberg
+** Date Started: June 2000
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: mpiworld.cpp,v 1.1 2000/07/13 07:01:35 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 <mpi++.h>
+#include <mpiworld.h>
+
+
+MPIWorld::MPIWorld (int& argc, char**& argv)
+{
+ MPI::Init (argc, argv);
+ m_comm = MPI::COMM_WORLD.Dup();
+ m_nProcessors = m_comm.Get_size();
+ m_myRank = m_comm.Get_rank();
+ m_vLocalWorkUnits.reserve (m_nProcessors);
+ m_vStartWorkUnit.reserve (m_nProcessors);
+ m_vEndWorkUnit.reserve (m_nProcessors);
+}
+
+
+void
+MPIWorld::setTotalWorkUnits(int totalWorkUnits)
+{
+ if (m_nProcessors < 1)
+ return;
+
+ int baseLocalWorkUnits = totalWorkUnits / m_nProcessors;
+ int remainderWorkUnits = totalWorkUnits % m_nProcessors;
+
+ int currWorkUnits = 0;
+ for (int iProc = 0; iProc < m_nProcessors; iProc++) {
+ m_vLocalWorkUnits[iProc] = baseLocalWorkUnits;
+ if (iProc < remainderWorkUnits)
+ m_vLocalWorkUnits[iProc]++;
+
+ m_vStartWorkUnit[iProc] = currWorkUnits;
+ m_vEndWorkUnit[iProc] = m_vStartWorkUnit[iProc] + m_vLocalWorkUnits[iProc] - 1;
+
+ currWorkUnits += m_vLocalWorkUnits[iProc];
+ }
+
+}
+
+void
+MPIWorld::BcastString (string& str)
+{
+ int len;
+
+ if (m_myRank == 0)
+ len = str.length();
+ m_comm.Bcast (&len, 1, MPI::INT, 0);
+ char buf [ len + 1];
+
+ if (m_myRank == 0)
+ strcpy (buf, str.c_str());
+
+ m_comm.Bcast (buf, len + 1, MPI::CHAR, 0);
+
+ if (m_myRank > 0)
+ str = buf;
+}
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: phm2if.cpp
+** Purpose: Convert an phantom object to an image file
+** Programmer: Kevin Rosenberg
+** Date Started: 1984
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: phm2if.cpp,v 1.1 2000/07/13 07:01:35 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_NSAMPLE, O_FILTER, O_TRACE, O_VERBOSE, O_HELP,
+ O_PHMFILE, O_FILTER_DOMAIN, O_FILTER_BW, O_FILTER_PARAM, O_DEBUG, O_VERSION };
+
+static struct option my_options[] =
+{
+ {"phantom", 1, 0, O_PHANTOM},
+ {"phmfile", 1, 0, O_PHMFILE},
+ {"desc", 1, 0, O_DESC},
+ {"nsample", 1, 0, O_NSAMPLE},
+ {"filter", 1, 0, O_FILTER},
+ {"filter-domain", 1, 0, O_FILTER_DOMAIN},
+ {"filter-bw", 1, 0, O_FILTER_BW},
+ {"filter-param", 1, 0, O_FILTER_PARAM},
+ {"trace", 1, 0, O_TRACE},
+ {"verbose", 0, 0, O_VERBOSE},
+ {"debug", 0, 0, O_DEBUG},
+ {"help", 0, 0, O_HELP},
+ {"version", 0, 0, O_VERSION},
+ {0, 0, 0, 0}
+};
+
+void
+phm2if_usage (const char *program)
+{
+ cout << "phm2if_usage: " << fileBasename(program) << " outfile nx ny [--phantom phantom-name] [--phmfile filename] [--filter filter-name] [OPTIONS]" << endl;
+ cout << "Generate phantom image from a predefined --phantom or a --phmfile" << endl;
+ cout << endl;
+ cout << " outfile Name of output file for image" << endl;
+ cout << " nx Number of pixels X-axis" << endl;
+ cout << " ny Number of pixels Y-axis" << endl;
+ cout << " --phantom Phantom to use for projection" << endl;
+ cout << " herman Herman head phantom" << endl;
+ cout << " bherman Bordered 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 Generate Phantom from phantom file" << endl;
+ cout << " --filter Generate Phantom from a filter function" << endl;
+ cout << " abs_bandlimit Abs * Bandlimiting" << endl;
+ cout << " abs_sinc Abs * Sinc" << endl;
+ cout << " abs_cos Abs * Cosine" << endl;
+ cout << " abs_hamming Abs * Hamming" << endl;
+ cout << " shepp Shepp-Logan" << endl;
+ cout << " bandlimit Bandlimiting" << endl;
+ cout << " sinc Sinc" << endl;
+ cout << " cos Cosine" << endl;
+ cout << " triangle Triangle" << endl;
+ cout << " hamming Hamming" << endl;
+ cout << " --filter-param Alpha level for Hamming filter" << endl;
+ cout << " --filter-domain Set domain of filter" << endl;
+ cout << " spatial Spatial domain (default)" << endl;
+ cout << " freq Frequency domain" << endl;
+ cout << " --filter-bw Filter bandwidth (default = 1)" << endl;
+ cout << " --desc Description of raysum" << endl;
+ cout << " --nsample Number of samples per axis per pixel (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" << endl;
+ cout << " rays Trace rays" << endl;
+ cout << " plot Trace plot" << endl;
+ cout << " clipping Trace clipping" << endl;
+ cout << " --debug Debug mode" << endl;
+ cout << " --verbose Verbose mode" << endl;
+ cout << " --version Print version" << endl;
+ cout << " --help Print this help message" << endl;
+}
+
+#ifdef HAVE_MPI
+void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* imGlobal, ImageFile* imLocal, const int opt_debug);
+#endif
+
+int
+phm2if_main (int argc, char* argv[])
+{
+ ImageFile* imGlobal = NULL;
+ Phantom phm;
+ int opt_nx = 0, opt_ny = 0;
+ int opt_nsample = 1;
+ string optPhmName = Phantom::PHM_HERMAN_STR;
+ string optFilterName;
+ string optDomainName = SignalFilter::DOMAIN_SPATIAL_STR;
+ char *opt_outfile = NULL;
+ int opt_debug = 0;
+ string opt_desc;
+ string opt_phmFileName;
+ char *endstr, *endptr;
+ double opt_filter_param = 1;
+ double opt_filter_bw = 1.;
+ int opt_trace = TRACE_NONE;
+ bool opt_verbose = false;
+#ifdef HAVE_MPI
+ ImageFile* imLocal = NULL;
+ MPIWorld mpiWorld (argc, argv);
+#endif
+
+ Timer timerProgram;
+
+#ifdef HAVE_MPI
+ if (mpiWorld.getRank() == 0) {
+#endif
+ while (1) {
+ int c = getopt_long(argc, argv, "", my_options, NULL);
+ char *endptr = NULL;
+ char *endstr;
+
+ if (c == -1)
+ break;
+
+ switch (c) {
+ case O_PHANTOM:
+ optPhmName = optarg;
+ break;
+ case O_PHMFILE:
+ opt_phmFileName = optarg;
+ break;
+ case O_VERBOSE:
+ opt_verbose = true;
+ break;
+ case O_DEBUG:
+ opt_debug = 1;
+ break;
+ case O_TRACE:
+ if ((opt_trace = convertTraceNameToID(optarg)) == TRACE_INVALID) {
+ phm2if_usage(argv[0]);
+ return (1);
+ }
+ break;
+ case O_FILTER:
+ optFilterName = optarg;
+ break;
+ case O_FILTER_DOMAIN:
+ optDomainName = optarg;
+ break;
+ case O_DESC:
+ opt_desc = optarg;
+ break;
+ case O_FILTER_BW:
+ opt_filter_bw = strtod(optarg, &endptr);
+ endstr = optarg + strlen(optarg);
+ if (endptr != endstr) {
+ sys_error(ERR_SEVERE,"Error setting --filter-bw to %s\n", optarg);
+ phm2if_usage(argv[0]);
+ return (1);
+ }
+ break;
+ case O_FILTER_PARAM:
+ opt_filter_param = strtod(optarg, &endptr);
+ endstr = optarg + strlen(optarg);
+ if (endptr != endstr) {
+ sys_error(ERR_SEVERE,"Error setting --filter-param to %s\n", optarg);
+ phm2if_usage(argv[0]);
+ return (1);
+ }
+ break;
+ case O_NSAMPLE:
+ opt_nsample = strtol(optarg, &endptr, 10);
+ endstr = optarg + strlen(optarg);
+ if (endptr != endstr) {
+ sys_error(ERR_SEVERE,"Error setting --nsample to %s\n", optarg);
+ phm2if_usage(argv[0]);
+ return (1);
+ }
+ break;
+ case O_VERSION:
+#ifdef VERSION
+ cout << "Version " << VERSION << endl;
+#else
+ cerr << "Unknown version number" << endl;
+#endif
+ case O_HELP:
+ case '?':
+ phm2if_usage(argv[0]);
+ return (0);
+ default:
+ phm2if_usage(argv[0]);
+ return (1);
+ }
+ }
+
+ if (optPhmName == "" && optFilterName == "" && opt_phmFileName == "") {
+ cerr << "No phantom defined" << endl << endl;
+ phm2if_usage(argv[0]);
+ return (1);
+ }
+
+ if (optind + 3 != argc) {
+ phm2if_usage(argv[0]);
+ return (1);
+ }
+ opt_outfile = argv[optind];
+ opt_nx = strtol(argv[optind+1], &endptr, 10);
+ endstr = argv[optind+1] + strlen(argv[optind+1]);
+ if (endptr != endstr) {
+ sys_error(ERR_SEVERE,"Error setting nx to %s\n", argv[optind+1]);
+ phm2if_usage(argv[0]);
+ return (1);
+ }
+ opt_ny = strtol(argv[optind+2], &endptr, 10);
+ endstr = argv[optind+2] + strlen(argv[optind+2]);
+ if (endptr != endstr) {
+ sys_error(ERR_SEVERE,"Error setting ny to %s\n", argv[optind+2]);
+ phm2if_usage(argv[0]);
+ return (1);
+ }
+
+ ostringstream oss;
+ oss << "phm2if: nx=" << opt_nx << ", ny=" << opt_ny << ", nsample=" << opt_nsample << ", ";
+ if (opt_phmFileName != "")
+ oss << "phantom=" << opt_phmFileName;
+ else if (optPhmName != "")
+ oss << "phantom=" << optPhmName;
+ else if (optFilterName != "") {
+ oss << "filter=" << optFilterName << " - " << optDomainName;
+ }
+ if (opt_desc != "")
+ oss << ": " << opt_desc;
+ opt_desc = oss.str();
+
+ if (optPhmName != "") {
+ phm.createFromPhantom (optPhmName.c_str());
+ if (phm.fail()) {
+ cout << phm.failMessage() << endl << endl;
+ phm2if_usage(argv[0]);
+ return (1);
+ }
+ }
+
+ if (opt_phmFileName != "") {
+ phm.createFromFile(opt_phmFileName.c_str());
+#ifdef HAVE_MPI
+ if (mpiWorld.getRank() == 0)
+ cerr << "Can't use phantom from file in MPI mode" << endl;
+ return (1);
+#endif
+ }
+
+ if (opt_verbose)
+ cout << "Rasterize Phantom to Image" << endl << endl;
+#ifdef HAVE_MPI
+ }
+#endif
+
+#ifdef HAVE_MPI
+ TimerCollectiveMPI timerBcast (mpiWorld.getComm());
+ mpiWorld.BcastString (optPhmName);
+ 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);
+ mpiWorld.getComm().Bcast (&opt_nx, 1, MPI::INT, 0);
+ mpiWorld.getComm().Bcast (&opt_ny, 1, MPI::INT, 0);
+ mpiWorld.getComm().Bcast (&opt_nsample, 1, MPI::INT, 0);
+ mpiWorld.getComm().Bcast (&opt_filter_param, 1, MPI::DOUBLE, 0);
+ mpiWorld.getComm().Bcast (&opt_filter_bw, 1, MPI::DOUBLE, 0);
+
+ mpiWorld.BcastString (optFilterName);
+ mpiWorld.BcastString (optDomainName);
+
+ if (opt_verbose)
+ timerBcast.timerEndAndReport ("Time to broadcast variables");
+
+ mpiWorld.setTotalWorkUnits (opt_nx);
+
+ if (mpiWorld.getRank() > 0 && optPhmName != "")
+ phm.createFromPhantom (optPhmName.c_str());
+
+ if (mpiWorld.getRank() == 0) {
+ imGlobal = new ImageFile (opt_nx, opt_ny);
+ }
+ imLocal = new ImageFile (opt_nx, opt_ny);
+#else
+ imGlobal = new ImageFile (opt_nx, opt_ny);
+#endif
+
+ ImageFileArray v;
+#ifdef HAVE_MPI
+ if (mpiWorld.getRank() == 0)
+ v = imGlobal->getArray ();
+
+ if (phm.getComposition() == P_UNIT_PULSE) {
+ if (mpiWorld.getRank() == 0) {
+ v[opt_nx/2][opt_ny/2] = 1.;
+ }
+ } else if (optFilterName != "") {
+ if (mpiWorld.getRank() == 0) {
+ imGlobal->filterResponse (optDomainName.c_str(), opt_filter_bw, optFilterName.c_str(), opt_filter_param);
+ }
+ } else {
+ TimerCollectiveMPI timerRasterize (mpiWorld.getComm());
+ phm.convertToImagefile (*imLocal, opt_nsample, opt_trace, mpiWorld.getMyStartWorkUnit(), mpiWorld.getMyLocalWorkUnits());
+ if (opt_verbose)
+ timerRasterize.timerEndAndReport ("Time to rasterize phantom");
+
+ TimerCollectiveMPI timerGather (mpiWorld.getComm());
+ mpi_gather_image (mpiWorld, imGlobal, imLocal, opt_debug);
+ if (opt_verbose)
+ timerGather.timerEndAndReport ("Time to gather image");
+ }
+#else
+ v = imGlobal->getArray ();
+ if (phm.getComposition() == P_UNIT_PULSE) {
+ v[opt_nx/2][opt_ny/2] = 1.;
+ } else if (optFilterName != "") {
+ imGlobal->filterResponse (optDomainName.c_str(), opt_filter_bw, optFilterName.c_str(), opt_filter_param);
+ } else {
+#if HAVE_SGP
+ if (opt_trace >= TRACE_PHM)
+ phm.show();
+#endif
+ phm.convertToImagefile (*imGlobal, opt_nsample, opt_trace);
+ }
+#endif
+
+#ifdef HAVE_MPI
+ if (mpiWorld.getRank() == 0)
+#endif
+ {
+ double calctime = timerProgram.timerEnd ();
+ imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, opt_desc.c_str(), calctime);
+ imGlobal->fileWrite (opt_outfile);
+ if (opt_verbose)
+ cout << "Time to rasterized phantom: " << calctime << " seconds" << endl;
+
+ if (opt_trace >= TRACE_PHM) {
+ double dmin, dmax;
+ int nscale;
+
+ printf ("Enter display size scale (nominal = 1): ");
+ scanf ("%d", &nscale);
+ printf ("Enter minimum and maximum densities (min, max): ");
+ scanf ("%lf %lf", &dmin, &dmax);
+ imGlobal->displayScaling (nscale, dmin, dmax);
+ }
+ }
+
+ return (0);
+}
+
+
+
+#ifdef HAVE_MPI
+void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* imGlobal, ImageFile* imLocal, const int opt_debug)
+{
+ ImageFileArray vLocal = imLocal->getArray();
+ ImageFileArray vGlobal = NULL;
+ int nyLocal = imLocal->ny();
+
+ if (mpiWorld.getRank() == 0)
+ vGlobal = imGlobal->getArray();
+
+ for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++)
+ mpiWorld.getComm().Send(vLocal[iw], nyLocal, imLocal->getMPIDataType(), 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;
+ mpiWorld.getComm().Recv(vGlobal[iw], nyLocal, imLocal->getMPIDataType(), iProc, 0, status);
+ }
+ }
+ }
+
+}
+#endif
+
+#ifndef NO_MAIN
+int
+main (int argc, char* argv[])
+{
+ int retval = 1;
+
+ try {
+ retval = phm2if_main(argc, argv);
+ } catch (exception e) {
+ cerr << "Exception: " << e.what() << endl;
+ } catch (...) {
+ cerr << "Unknown exception" << endl;
+ }
+
+ return (retval);
+}
+#endif
--- /dev/null
+/*****************************************************************************
+** 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/07/13 07:01:35 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 << " bherman Bordered 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;
+ string optGeometryName = Scanner::GEOMETRY_PARALLEL_STR;
+ char *opt_outfile = NULL;
+ string opt_desc;
+ string optPhmFileName;
+ int opt_ndet;
+ int opt_nview;
+ int opt_nray = 1;
+ int opt_trace = 0;
+ string optPhmName = Phantom::PHM_HERMAN_STR;
+ 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:
+ optPhmName = optarg;
+ break;
+ case O_PHMFILE:
+ optPhmFileName = optarg;
+ break;
+ case O_VERBOSE:
+ opt_verbose = 1;
+ break;
+ case O_DEBUG:
+ opt_debug = 1;
+ break;
+ break;
+ case O_TRACE:
+ if ((opt_trace = convertTraceNameToID(optarg)) == TRACE_INVALID) {
+ 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
+ return (0);
+ case O_HELP:
+ case '?':
+ phm2pj_usage(argv[0]);
+ return (0);
+ default:
+ phm2pj_usage(argv[0]);
+ return (1);
+ }
+ }
+
+ if (optPhmName == "" && optPhmFileName == "") {
+ cerr << "No phantom defined" << endl << 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 << "phm2pj: NDet=" << opt_ndet << ", Nview=" << opt_nview << ", NRay=" << opt_nray << ", RotAngle=" << opt_rotangle << ", ";
+ if (optPhmFileName.length()) {
+ desc << "PhantomFile=" << optPhmFileName;
+ } else if (optPhmName != "") {
+ desc << "Phantom=" << optPhmName;
+ }
+ if (opt_desc.length()) {
+ desc << ": " << opt_desc;
+ }
+ opt_desc = desc.str();
+
+ if (optPhmName != "") {
+ phm.createFromPhantom (optPhmName.c_str());
+ if (phm.fail()) {
+ cout << phm.failMessage() << endl << endl;
+ phm2pj_usage(argv[0]);
+ return (1);
+ }
+ }
+
+ if (optPhmFileName != "") {
+#ifdef HAVE_MPI
+ cerr << "Can not read phantom from file in MPI mode" << endl;
+ return (1);
+#endif
+ phm.createFromFile (optPhmFileName.c_str());
+ }
+
+#ifdef HAVE_MPI
+ }
+#endif
+
+#ifdef HAVE_MPI
+ TimerCollectiveMPI timerBcast(mpiWorld.getComm());
+ mpiWorld.BcastString (optPhmName);
+ 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_verbose, 1, MPI::INT, 0);
+ mpiWorld.getComm().Bcast (&opt_debug, 1, MPI::INT, 0);
+ mpiWorld.getComm().Bcast (&opt_trace, 1, MPI::INT, 0);
+ timerBcast.timerEndAndReport ("Time to broadcast variables");
+
+ if (mpiWorld.getRank() > 0 && optPhmName != "")
+ phm.createFromPhantom (optPhmName.c_str());
+#endif
+
+ opt_rotangle *= PI;
+ Scanner scanner (phm, optGeometryName.c_str(), opt_ndet, opt_nview, opt_nray, opt_rotangle);
+ if (scanner.fail()) {
+ cout << "Scanner Creation Error: " << scanner.failMessage() << endl;
+ return (1);
+ }
+#ifdef HAVE_MPI
+ mpiWorld.setTotalWorkUnits (opt_nview);
+
+ Projections pjGlobal;
+ if (mpiWorld.getRank() == 0)
+ pjGlobal.initFromScanner (scanner);
+
+ if (opt_verbose)
+ pjGlobal.printScanInfo();
+
+ Projections pjLocal (scanner);
+ pjLocal.setNView (mpiWorld.getMyLocalWorkUnits());
+
+ if (opt_debug)
+ cout << "pjLocal->nview = " << pjLocal.nView() << " (process " << mpiWorld.getRank() << ")" << endl;;
+
+ TimerCollectiveMPI timerProject (mpiWorld.getComm());
+ scanner.collectProjections (pjLocal, phm, mpiWorld.getMyStartWorkUnit(), opt_trace);
+ if (opt_verbose)
+ timerProject.timerEndAndReport ("Time to collect projections");
+
+ TimerCollectiveMPI 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[])
+{
+ int retval = 1;
+
+ try {
+ retval = phm2pj_main(argc, argv);
+ } catch (exception e) {
+ cerr << "Exception: " << e.what() << endl;
+ } catch (...) {
+ cerr << "Unknown exception" << endl;
+ }
+
+ return (retval);
+}
+#endif
+
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: pj2if.cpp
+** Purpose: Convert an projection data file to an image file
+** Programmer: Kevin Rosenberg
+** Date Started: April 2000
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: pj2if.cpp,v 1.1 2000/07/13 07:01:35 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
+******************************************************************************/
+
+/* FILE
+ * pj2if.c Convert Raysum to image
+ *
+ * DATE
+ * Apr 1999
+ */
+
+#include "ct.h"
+#include "timer.h"
+
+
+enum { O_VERBOSE, O_HELP, O_VERSION };
+
+static struct option my_options[] =
+{
+ {"verbose", 0, 0, O_VERBOSE},
+ {"help", 0, 0, O_HELP},
+ {"version", 0, 0, O_VERSION},
+ {0, 0, 0, 0}
+};
+
+
+void
+pj2if_usage (const char *program)
+{
+ cout << "usage: " << fileBasename(program) << " in-proj-file out-if-file [OPTIONS]" << endl;
+ cout << "Converts a projection file to a IF file" << endl;
+ cout << endl;
+ cout << " --verbose Verbose mode" << endl;
+ cout << " --version Print version" << endl;
+ cout << " --help Print this help message" << endl;
+}
+
+
+
+int
+pj2if_main (const int argc, char *const argv[])
+{
+ char *pj_name, *im_name;
+ bool optVerbose = false;
+ extern int optind;
+ Timer timerProgram;
+
+ while (1)
+ {
+ int c = getopt_long (argc, argv, "", my_options, NULL);
+ if (c == -1)
+ break;
+
+ switch (c)
+ {
+ case O_VERBOSE:
+ optVerbose = true;
+ break;
+ case O_VERSION:
+#ifdef VERSION
+ cout << "Version " << VERSION << endl;
+#else
+ cout << "Unknown version number" << endl;
+#endif
+ return (0);
+ case O_HELP:
+ case '?':
+ pj2if_usage(argv[0]);
+ return (0);
+ default:
+ pj2if_usage(argv[0]);
+ return (1);
+ }
+ }
+
+ if (argc - optind != 2) {
+ pj2if_usage(argv[0]);
+ return (1);
+ }
+
+ pj_name = argv[optind];
+ im_name = argv[optind + 1];
+
+ Projections pj;
+ if (! pj.read (pj_name)) {
+ sys_error (ERR_SEVERE, "Can not open projection file %s", pj_name);
+ return (1);
+ }
+
+ if (optVerbose)
+ pj.printScanInfo();
+
+ ImageFile im (pj.nDet(), pj.nView());
+
+ ImageFileArray v = im.getArray();
+
+ for (int iy = 0; iy < pj.nView(); iy++) {
+ const DetectorArray& detarray = pj.getDetectorArray (iy);
+ const DetectorValue* detval = detarray.detValues();
+ for (int ix = 0; ix < pj.nDet(); ix++) {
+ v[ix][iy] = detval[ix];
+ }
+ }
+
+ im.labelAdd (pj.getLabel());
+ im.labelAdd (Array2dFileLabel::L_HISTORY, "Conversion from .pj to .if", timerProgram.timerEnd());
+ im.fileWrite (im_name);
+
+ return(0);
+}
+
+
+#ifndef NO_MAIN
+int
+main (const int argc, char *const argv[])
+{
+ int retval = 1;
+
+ try {
+ retval = pj2if_main(argc, argv);
+ } catch (exception e) {
+ cerr << "Exception: " << e.what() << endl;
+ } catch (...) {
+ cerr << "Unknown exception" << endl;
+ }
+
+ return (retval);
+}
+#endif
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: pjrec.cpp
+** Purpose: Reconstruct an image from projections
+** Programmer: Kevin Rosenberg
+** Date Started: Aug 1984
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: pjrec.cpp,v 1.1 2000/07/13 07:01:35 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_INTERP, O_FILTER, O_FILTER_METHOD, O_ZEROPAD, O_FILTER_PARAM, O_BACKPROJ, O_PREINTERPOLATION_FACTOR, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION};
+
+static struct option my_options[] =
+{
+ {"interp", 1, 0, O_INTERP},
+ {"preinterpolation-factor", 1, 0, O_PREINTERPOLATION_FACTOR},
+ {"filter", 1, 0, O_FILTER},
+ {"filter-method", 1, 0, O_FILTER_METHOD},
+ {"zeropad", 1, 0, O_ZEROPAD},
+ {"filter-param", 1, 0, O_FILTER_PARAM},
+ {"backproj", 1, 0, O_BACKPROJ},
+ {"trace", 1, 0, O_TRACE},
+ {"debug", 0, 0, O_DEBUG},
+ {"verbose", 0, 0, O_VERBOSE},
+ {"help", 0, 0, O_HELP},
+ {"version", 0, 0, O_VERSION},
+ {0, 0, 0, 0}
+};
+
+
+void
+pjrec_usage (const char *program)
+{
+ cout << "usage: " << fileBasename(program) << " raysum-file image-file nx-image ny-image [OPTIONS]" << endl;
+ cout << "Image reconstruction from raysum projections" << endl;
+ cout << endl;
+ cout << " raysum-file Input raysum file" << endl;
+ cout << " image-file Output image file in SDF2D format" << endl;
+ cout << " nx-image Number of columns in output image" << endl;
+ cout << " ny-image Number of rows in output image" << endl;
+ cout << " --interp Interpolation method during backprojection" << endl;
+ cout << " nearest Nearest neighbor interpolation" << endl;
+ cout << " linear Linear interpolation" << endl;
+#if HAVE_BSPLINE_INTERP
+ cout << " bspline B-spline interpolation" << endl;
+#endif
+ cout << " --preinterpolate Preinterpolation factor (default = 1)\n";
+ cout << " Used only with frequency-based filtering\n";
+ cout << " --filter Filter name" << endl;
+ cout << " abs_bandlimit Abs * Bandlimiting (default)" << endl;
+ cout << " abs_sinc Abs * Sinc" << endl;
+ cout << " abs_cos Abs * Cosine" << endl;
+ cout << " abs_hamming Abs * Hamming" << endl;
+ cout << " shepp Shepp-Logan" << endl;
+ cout << " bandlimit Bandlimiting" << endl;
+ cout << " sinc Sinc" << endl;
+ cout << " cos Cosine" << endl;
+ cout << " triangle Triangle" << endl;
+ cout << " hamming Hamming" << endl;
+ cout << " --filter-method Filter method before backprojections\n";;
+ cout << " convolution Spatial filtering (default)\n";
+ cout << " fourier Frequency filtering with discete fourier\n";
+ cout << " fourier_table Frequency filtering with table lookup fourier\n";
+ cout << " fft Fast Fourier Transform\n";
+#if HAVE_FFTW
+ cout << " fftw Fast Fourier Transform West library\n";
+ cout << " rfftw Fast Fourier Transform West (real-mode) library\n";
+#endif
+ cout << " --zeropad n Set zeropad level (default = 0)\n";
+ cout << " set n to number of powers to two to pad\n";
+ cout << " --backproj Backprojection Method" << endl;
+ cout << " trig Trigometric functions at every point" << endl;
+ cout << " table Trigometric functions with precalculated table" << endl;
+ cout << " diff Difference method" << endl;
+ cout << " diff2 Optimized difference method (default)" << endl;
+ cout << " idiff2 Optimized difference method with integer math" << endl;
+ cout << " idiff3 Highly-optimized difference method with integer math" << endl;
+ cout << " --filter-param Alpha level for Hamming filter" << endl;
+ cout << " --trace Set tracing to level" << endl;
+ cout << " none No tracing (default)" << endl;
+ cout << " text Text level tracing" << endl;
+ cout << " phm Trace phantom" << endl;
+ cout << " rays Trace allrays" << endl;
+ cout << " plot Trace plotting" << endl;
+ cout << " clipping Trace clipping" << endl;
+ cout << " --verbose Turn on verbose mode" << endl;
+ cout << " --debug Turn on debug mode" << endl;
+ cout << " --version Print version" << endl;
+ cout << " --help Print this help message" << endl;
+}
+
+
+#ifdef HAVE_MPI
+static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int debug);
+static void ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* imGlobal);
+#endif
+
+
+int
+pjrec_main (int argc, char * argv[])
+{
+ Projections projGlobal;
+ ImageFile* imGlobal = NULL;
+ char* filenameProj = NULL;
+ char* filenameImage = NULL;
+ string remark;
+ char *endptr;
+ int optVerbose = 0;
+ int optDebug = 0;
+ int optZeroPad = 0;
+ int optTrace = TRACE_NONE;
+ double optFilterParam = -1;
+ string optFilterName = SignalFilter::FILTER_ABS_BANDLIMIT_STR;
+ string optFilterMethodName = SignalFilter::FILTER_METHOD_CONVOLUTION_STR;
+ string optInterpName = Backprojector::INTERP_LINEAR_STR;
+ string optBackprojName = Backprojector::BPROJ_IDIFF2_STR;
+ // string optFilterMethodName = SignalFilter::FILTER_METHOD_CONVOLUTION_STR;
+ int optPreinterpolationFactor = 1;
+ int nx, ny;
+#ifdef HAVE_MPI
+ ImageFile* imLocal;
+ int mpi_nview, mpi_ndet;
+ double mpi_detinc, mpi_rotinc, mpi_phmlen;
+ MPIWorld mpiWorld (argc, argv);
+#endif
+
+ Timer timerProgram;
+
+#ifdef HAVE_MPI
+ if (mpiWorld.getRank() == 0) {
+#endif
+ while (1) {
+ int c = getopt_long(argc, argv, "", my_options, NULL);
+ char *endptr = NULL;
+
+ if (c == -1)
+ break;
+
+ switch (c)
+ {
+ case O_INTERP:
+ optInterpName = optarg;
+ break;
+ case O_PREINTERPOLATION_FACTOR:
+ optPreinterpolationFactor = strtol(optarg, &endptr, 10);
+ if (endptr != optarg + strlen(optarg)) {
+ pjrec_usage(argv[0]);
+ return(1);
+ }
+ break;
+ case O_FILTER:
+ optFilterName = optarg;
+ break;
+ case O_FILTER_METHOD:
+ optFilterMethodName = optarg;
+ break;
+ case O_BACKPROJ:
+ optBackprojName = optarg;
+ break;
+ case O_FILTER_PARAM:
+ optFilterParam = strtod(optarg, &endptr);
+ if (endptr != optarg + strlen(optarg)) {
+ pjrec_usage(argv[0]);
+ return(1);
+ }
+ break;
+ case O_ZEROPAD:
+ optZeroPad = strtol(optarg, &endptr, 10);
+ if (endptr != optarg + strlen(optarg)) {
+ pjrec_usage(argv[0]);
+ return(1);
+ }
+ break;
+ case O_VERBOSE:
+ optVerbose = 1;
+ break;
+ case O_DEBUG:
+ optDebug = 1;
+ break;
+ case O_TRACE:
+ if ((optTrace = convertTraceNameToID(optarg)) == TRACE_INVALID) {
+ pjrec_usage(argv[0]);
+ return (1);
+ }
+ break;
+ case O_VERSION:
+#ifdef VERSION
+ cout << "Version " << VERSION << endl;
+#else
+ cout << "Unknown version number" << endl;
+#endif
+ return (0);
+ case O_HELP:
+ case '?':
+ pjrec_usage(argv[0]);
+ return (0);
+ default:
+ pjrec_usage(argv[0]);
+ return (1);
+ }
+ }
+
+ if (optind + 4 != argc) {
+ pjrec_usage(argv[0]);
+ return (1);
+ }
+
+ filenameProj = argv[optind];
+
+ filenameImage = argv[optind + 1];
+
+ nx = strtol(argv[optind + 2], &endptr, 10);
+ ny = strtol(argv[optind + 3], &endptr, 10);
+
+ ostringstream filterDesc;
+ if (optFilterParam >= 0)
+ filterDesc << optFilterName << ": alpha=" << optFilterParam;
+ else
+ filterDesc << optFilterName;
+
+ ostringstream label;
+ label << "pjrec: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << optInterpName << ", preinterpolation=" << optPreinterpolationFactor << ", " << optBackprojName;
+ remark = label.str();
+
+ if (optVerbose)
+ cout << "Remark: " << remark << endl;
+#ifdef HAVE_MPI
+ }
+#endif
+
+#ifdef HAVE_MPI
+ if (mpiWorld.getRank() == 0) {
+ projGlobal.read (filenameProj);
+ if (optVerbose)
+ projGlobal.printScanInfo();
+
+ mpi_ndet = projGlobal.nDet();
+ mpi_nview = projGlobal.nView();
+ mpi_detinc = projGlobal.detInc();
+ mpi_phmlen = projGlobal.phmLen();
+ mpi_rotinc = projGlobal.rotInc();
+ }
+
+ TimerCollectiveMPI timerBcast (mpiWorld.getComm());
+ mpiWorld.BcastString (optBackprojName);
+ mpiWorld.BcastString (optFilterName);
+ mpiWorld.BcastString (optFilterMethodName);
+ mpiWorld.BcastString (optInterpName);
+ mpiWorld.getComm().Bcast (&optVerbose, 1, MPI::INT, 0);
+ mpiWorld.getComm().Bcast (&optDebug, 1, MPI::INT, 0);
+ mpiWorld.getComm().Bcast (&optTrace, 1, MPI::INT, 0);
+ mpiWorld.getComm().Bcast (&optFilterParam, 1, MPI::DOUBLE, 0);
+ mpiWorld.getComm().Bcast (&optZeroPad, 1, MPI::INT, 0);
+ mpiWorld.getComm().Bcast (&optPreinterpolationFactor, 1, MPI::INT, 0);
+ mpiWorld.getComm().Bcast (&mpi_ndet, 1, MPI::INT, 0);
+ mpiWorld.getComm().Bcast (&mpi_nview, 1, MPI::INT, 0);
+ mpiWorld.getComm().Bcast (&mpi_detinc, 1, MPI::DOUBLE, 0);
+ mpiWorld.getComm().Bcast (&mpi_phmlen, 1, MPI::DOUBLE, 0);
+ mpiWorld.getComm().Bcast (&mpi_rotinc, 1, MPI::DOUBLE, 0);
+ mpiWorld.getComm().Bcast (&nx, 1, MPI::INT, 0);
+ mpiWorld.getComm().Bcast (&ny, 1, MPI::INT, 0);
+ if (optVerbose)
+ timerBcast.timerEndAndReport ("Time to broadcast variables");
+
+ mpiWorld.setTotalWorkUnits (mpi_nview);
+
+ Projections projLocal (mpiWorld.getMyLocalWorkUnits(), mpi_ndet);
+ projLocal.setDetInc (mpi_detinc);
+ projLocal.setPhmLen (mpi_phmlen);
+ projLocal.setRotInc (mpi_rotinc);
+
+ TimerCollectiveMPI timerScatter (mpiWorld.getComm());
+ ScatterProjectionsMPI (mpiWorld, projGlobal, projLocal, optDebug);
+ if (optVerbose)
+ timerScatter.timerEndAndReport ("Time to scatter projections");
+
+ if (mpiWorld.getRank() == 0) {
+ imGlobal = new ImageFile (nx, ny);
+ }
+
+ imLocal = new ImageFile (nx, ny);
+#else
+ projGlobal.read (filenameProj);
+ if (optVerbose)
+ projGlobal.printScanInfo();
+
+ imGlobal = new ImageFile (nx, ny);
+#endif
+
+#ifdef HAVE_MPI
+ TimerCollectiveMPI timerReconstruct (mpiWorld.getComm());
+ projLocal.reconstruct (*imLocal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optPreinterpolationFactor, optBackprojName.c_str(), optTrace);
+ if (optVerbose)
+ timerReconstruct.timerEndAndReport ("Time to reconstruct");
+
+ TimerCollectiveMPI timerReduce (mpiWorld.getComm());
+ ReduceImageMPI (mpiWorld, imLocal, imGlobal);
+ if (optVerbose)
+ timerReduce.timerEndAndReport ("Time to reduce image");
+#else
+ projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optPreinterpolationFactor, optBackprojName.c_str(), optTrace);
+#endif
+
+#ifdef HAVE_MPI
+ if (mpiWorld.getRank() == 0)
+#endif
+ {
+ double calcTime = timerProgram.timerEnd();
+ imGlobal->labelAdd (projGlobal.getLabel());
+ imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, remark.c_str(), calcTime);
+ imGlobal->fileWrite (filenameImage);
+ if (optVerbose)
+ cout << "Run time: " << calcTime << " seconds" << endl;
+ }
+#ifdef HAVE_MPI
+ MPI::Finalize();
+#endif
+
+ return (0);
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////
+// MPI Support Routines
+//
+//////////////////////////////////////////////////////////////////////////////////////
+
+#ifdef HAVE_MPI
+static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int optDebug)
+{
+ if (mpiWorld.getRank() == 0) {
+ for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) {
+ for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) {
+ DetectorArray& detarray = projGlobal.getDetectorArray( iw );
+ int nDet = detarray.nDet();
+ DetectorValue* detval = detarray.detValues();
+
+ double viewAngle = detarray.viewAngle();
+ mpiWorld.getComm().Send(&nDet, 1, MPI::INT, iProc, 0);
+ mpiWorld.getComm().Send(&viewAngle, 1, MPI::DOUBLE, iProc, 0);
+ mpiWorld.getComm().Send(detval, nDet, MPI::FLOAT, iProc, 0);
+ }
+ }
+ }
+
+ for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++) {
+ MPI::Status status;
+ int nDet;
+ double viewAngle;
+ DetectorValue* detval = projLocal.getDetectorArray(iw).detValues();
+
+ mpiWorld.getComm().Recv(&nDet, 1, MPI::INT, 0, 0, status);
+ mpiWorld.getComm().Recv(&viewAngle, 1, MPI::DOUBLE, 0, 0, status);
+ mpiWorld.getComm().Recv(detval, nDet, MPI::FLOAT, 0, 0, status);
+ projLocal.getDetectorArray(iw).setViewAngle( viewAngle );
+ }
+}
+
+static void
+ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* imGlobal)
+{
+ ImageFileArray vLocal = imLocal->getArray();
+
+ for (int ix = 0; ix < imLocal->nx(); ix++) {
+ void *recvbuf = NULL;
+ if (mpiWorld.getRank() == 0) {
+ ImageFileArray vGlobal = imGlobal->getArray();
+ recvbuf = vGlobal[ix];
+ }
+ mpiWorld.getComm().Reduce (vLocal[ix], recvbuf, imLocal->ny(), imLocal->getMPIDataType(), MPI::SUM, 0);
+ }
+}
+
+#endif
+
+
+#ifndef NO_MAIN
+int
+main (int argc, char* argv[])
+{
+ int retval = 1;
+
+ try {
+ retval = pjrec_main(argc, argv);
+ } catch (exception e) {
+ cerr << "Exception: " << e.what() << endl;
+ } catch (...) {
+ cerr << "Unknown exception" << endl;
+ }
+
+ return (retval);
+}
+#endif
+
projects / ctsim.git / commitdiff