From: Kevin M. Rosenberg Date: Thu, 13 Jul 2000 07:01:59 +0000 (+0000) Subject: r144: Initial CVS import X-Git-Tag: debian-4.5.3-3~873 X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=commitdiff_plain;h=1fd4f7cc977b9f1499716de10d15656bd50f4816 r144: Initial CVS import --- diff --git a/src/Makefile.am b/src/Makefile.am index 168ab30..9929768 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,63 +1,8 @@ -#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@ diff --git a/src/ctsim.cpp b/src/ctsim.cpp new file mode 100644 index 0000000..0f7bc4e --- /dev/null +++ b/src/ctsim.cpp @@ -0,0 +1,131 @@ +/***************************************************************************** +** 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 \nUsage: ctsim", "About CTSim", wxOK | wxICON_INFORMATION, this); +} + +void +MainFrame::OnExit (wxCommandEvent& WXUNUSED(event) ) +{ + Close(true); +} + diff --git a/src/ctsim.h b/src/ctsim.h new file mode 100644 index 0000000..db4d7c3 --- /dev/null +++ b/src/ctsim.h @@ -0,0 +1,76 @@ +/***************************************************************************** +** 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 diff --git a/src/docs.cpp b/src/docs.cpp new file mode 100644 index 0000000..d5abb54 --- /dev/null +++ b/src/docs.cpp @@ -0,0 +1,117 @@ +/***************************************************************************** +** 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); +} diff --git a/src/docs.h b/src/docs.h new file mode 100644 index 0000000..1cd4c40 --- /dev/null +++ b/src/docs.h @@ -0,0 +1,100 @@ +/***************************************************************************** +** 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 diff --git a/src/if-1.cpp b/src/if-1.cpp deleted file mode 100644 index b1d352b..0000000 --- a/src/if-1.cpp +++ /dev/null @@ -1,203 +0,0 @@ -/***************************************************************************** -** 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 - diff --git a/src/if-2.cpp b/src/if-2.cpp deleted file mode 100644 index 7757f72..0000000 --- a/src/if-2.cpp +++ /dev/null @@ -1,310 +0,0 @@ -/***************************************************************************** -** 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 - diff --git a/src/if2img.cpp b/src/if2img.cpp deleted file mode 100644 index e6ec073..0000000 --- a/src/if2img.cpp +++ /dev/null @@ -1,388 +0,0 @@ -/***************************************************************************** -** 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 diff --git a/src/ifinfo.cpp b/src/ifinfo.cpp deleted file mode 100644 index 6149c9e..0000000 --- a/src/ifinfo.cpp +++ /dev/null @@ -1,161 +0,0 @@ -/***************************************************************************** -** 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 diff --git a/src/phm2if.cpp b/src/phm2if.cpp deleted file mode 100644 index afac816..0000000 --- a/src/phm2if.cpp +++ /dev/null @@ -1,413 +0,0 @@ -/***************************************************************************** -** 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 diff --git a/src/phm2pj.cpp b/src/phm2pj.cpp deleted file mode 100644 index 2734b2d..0000000 --- a/src/phm2pj.cpp +++ /dev/null @@ -1,370 +0,0 @@ -/***************************************************************************** -** 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 - diff --git a/src/pj2if.cpp b/src/pj2if.cpp deleted file mode 100644 index 25b72ca..0000000 --- a/src/pj2if.cpp +++ /dev/null @@ -1,154 +0,0 @@ -/***************************************************************************** -** 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 diff --git a/src/pjrec.cpp b/src/pjrec.cpp deleted file mode 100644 index d6a229c..0000000 --- a/src/pjrec.cpp +++ /dev/null @@ -1,415 +0,0 @@ -/***************************************************************************** -** 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 - diff --git a/src/sample-ctsim.sh.in b/src/sample-ctsim.sh.in deleted file mode 100755 index 1277641..0000000 --- a/src/sample-ctsim.sh.in +++ /dev/null @@ -1,41 +0,0 @@ -#!/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 diff --git a/src/views.cpp b/src/views.cpp new file mode 100644 index 0000000..8a555ce --- /dev/null +++ b/src/views.cpp @@ -0,0 +1,527 @@ +/***************************************************************************** +** 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 +#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(GetDocument())->getImageFile(); + const string& rFilename = rIF.getFilename(); + rIF.statistics (min, max, mean, mode, median, stddev); + ostringstream os; + os << "file: " << rFilename << "\nmin: "<GetClientSize(&width, &height); + + pCanvas = new ImageFileCanvas (dynamic_cast(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(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(GetDocument())->getImageFile(); + dc->Blit(static_cast(0), static_cast(0), static_cast(rIF.nx()), static_cast(rIF.ny()), &m_memoryDC, static_cast(0), static_cast(0)); +} + + +void +ImageFileView::OnUpdate(wxView *WXUNUSED(sender), wxObject *WXUNUSED(hint) ) +{ + const ImageFile& rIF = dynamic_cast(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(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(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(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(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(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(GetDocument())->getProjections(); + dc->Blit(static_cast(0), static_cast(0), static_cast(rProj.nDet()), static_cast(rProj.nView()), &m_memoryDC, static_cast(0), static_cast(0)); +} + + +void +ProjectionFileView::OnUpdate(wxView *WXUNUSED(sender), wxObject *WXUNUSED(hint) ) +{ + const Projections& rProj = dynamic_cast(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(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; +} + diff --git a/src/views.h b/src/views.h new file mode 100644 index 0000000..4e833bb --- /dev/null +++ b/src/views.h @@ -0,0 +1,156 @@ +/***************************************************************************** +** 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 diff --git a/tools/Makefile.am b/tools/Makefile.am new file mode 100644 index 0000000..2b1cf83 --- /dev/null +++ b/tools/Makefile.am @@ -0,0 +1,62 @@ +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 + + + diff --git a/tools/if-1.cpp b/tools/if-1.cpp new file mode 100644 index 0000000..a069219 --- /dev/null +++ b/tools/if-1.cpp @@ -0,0 +1,203 @@ +/***************************************************************************** +** 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 + diff --git a/tools/if-2.cpp b/tools/if-2.cpp new file mode 100644 index 0000000..525330e --- /dev/null +++ b/tools/if-2.cpp @@ -0,0 +1,312 @@ +/***************************************************************************** +** 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 + diff --git a/tools/if2img.cpp b/tools/if2img.cpp new file mode 100644 index 0000000..b00c177 --- /dev/null +++ b/tools/if2img.cpp @@ -0,0 +1,388 @@ +/***************************************************************************** +** 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 diff --git a/tools/ifinfo.cpp b/tools/ifinfo.cpp new file mode 100644 index 0000000..09ddfd6 --- /dev/null +++ b/tools/ifinfo.cpp @@ -0,0 +1,161 @@ +/***************************************************************************** +** 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 diff --git a/tools/mpiworld.cpp b/tools/mpiworld.cpp new file mode 100644 index 0000000..87c9a93 --- /dev/null +++ b/tools/mpiworld.cpp @@ -0,0 +1,84 @@ +/***************************************************************************** +** 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 +#include + + +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; +} diff --git a/tools/phm2if.cpp b/tools/phm2if.cpp new file mode 100644 index 0000000..9b8fd22 --- /dev/null +++ b/tools/phm2if.cpp @@ -0,0 +1,413 @@ +/***************************************************************************** +** 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 diff --git a/tools/phm2pj.cpp b/tools/phm2pj.cpp new file mode 100644 index 0000000..69c6ea5 --- /dev/null +++ b/tools/phm2pj.cpp @@ -0,0 +1,370 @@ +/***************************************************************************** +** 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 + diff --git a/tools/pj2if.cpp b/tools/pj2if.cpp new file mode 100644 index 0000000..1e8ea40 --- /dev/null +++ b/tools/pj2if.cpp @@ -0,0 +1,152 @@ +/***************************************************************************** +** 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 diff --git a/tools/pjrec.cpp b/tools/pjrec.cpp new file mode 100644 index 0000000..dc00484 --- /dev/null +++ b/tools/pjrec.cpp @@ -0,0 +1,415 @@ +/***************************************************************************** +** 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 +