From: Kevin M. Rosenberg Date: Tue, 13 Jun 2000 16:20:31 +0000 (+0000) Subject: r94: finished c++ conversions X-Git-Tag: debian-4.5.3-3~923 X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=commitdiff_plain;h=031437896d0dc6cac70c16e5604b10f5aa4d0767;ds=sidebyside r94: finished c++ conversions --- diff --git a/configure.in b/configure.in index 9e32c1f..8693d57 100644 --- a/configure.in +++ b/configure.in @@ -4,7 +4,7 @@ dnl Must reset CDPATH so that bash's cd does not print to stdout dnl CDPATH= AC_INIT(src/ctrec.cpp) -AM_INIT_AUTOMAKE(ctsim,0.6.0-b2) +AM_INIT_AUTOMAKE(ctsim,0.6.0-b3) AM_CONFIG_HEADER(config.h) dnl Checks for programs. @@ -292,21 +292,21 @@ else fi if test "$png" = "true" ; then - ctlibs_graphics="$ctlibs_graphics -lpng" + ctlibs_tools="$ctlibs_tools -lpng" fi if test "$zlib" = "true" ; then - ctlibs_graphics="$ctlibs_graphics -lz" + ctlibs_tools="$ctlibs_tools -lz" fi if test "$g2" = "true" ; then - ctlibs_graphics="$ctlibs_graphics -lg2" + ctlibs_tools="$ctlibs_tools -lg2" fi if test "$wx_gtk" = "true" ; then - ctlibs_graphics="$ctlibs_graphics -lwx_gtk" + ctlibs_tools="$ctlibs_tools -lwx_gtk" fi dnl Setting projet libraries and includes LDFLAGS="$LDFLAGS -L../libkmath -L../libk -L../libcio -L../libir" -ctlibs="$ctlibs_base -lir $ctlibs_graphics -lkmath -lk -lcio" +ctlibs="$ctlibs_base -lir $ctlibs_graphics -lkmath -lk -lcio $ctlibs_tools" AC_SUBST(ctlibs) if test -n "$lamdir" ; then diff --git a/include/backprojectors.h b/include/backprojectors.h index 923cee1..ea0e85c 100644 --- a/include/backprojectors.h +++ b/include/backprojectors.h @@ -9,8 +9,11 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: backprojectors.h,v 1.1 2000/06/10 23:00:17 kevin Exp $ +** $Id: backprojectors.h,v 1.2 2000/06/13 16:20:31 kevin Exp $ ** $Log: backprojectors.h,v $ +** Revision 1.2 2000/06/13 16:20:31 kevin +** finished c++ conversions +** ** Revision 1.1 2000/06/10 23:00:17 kevin ** *** empty log message *** ** @@ -49,12 +52,12 @@ class Backproject ImageFileArray v; kuint32 nx; kuint32 ny; - double det_inc; - double rot_inc; - int ncent; // index refering to L=0 projection - int ndet; - double xmin, xmax, ymin, ymax; // Retangular coords of phantom - double xinc, yinc; // size of cells + double detInc; + double rotInc; + int iDetCenter; // index refering to L=0 projection + int nDet; + double xMin, xMax, yMin, yMax; // Retangular coords of phantom + double xInc, yInc; // size of cells private: Backproject (const Backproject& rhs); diff --git a/include/cio.h b/include/cio.h index 1d14b89..4e83202 100644 --- a/include/cio.h +++ b/include/cio.h @@ -2,26 +2,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: cio.h,v 1.8 2000/06/05 01:32:45 kevin Exp $ -** $Log: cio.h,v $ -** Revision 1.8 2000/06/05 01:32:45 kevin -** Added C++ compatibility -** -** Revision 1.7 2000/05/08 20:45:10 kevin -** *** empty log message *** -** -** Revision 1.6 2000/05/08 20:00:47 kevin -** ANSI C changes -** -** Revision 1.5 2000/05/07 12:46:19 kevin -** made c++ compatible -** -** Revision 1.4 2000/04/28 18:18:59 kevin -** removed unused files -** -** Revision 1.3 2000/04/28 14:14:16 kevin -** *** empty log message *** -** +** $Id: cio.h,v 1.9 2000/06/13 16:20:31 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 @@ -39,10 +20,6 @@ #ifndef __CIO_H #define __CIO_H -#ifdef __cplusplus -extern "C" { -#endif /* __cplusplus */ - #define C_BLACK 0 /* color codes */ #define C_BLUE 1 @@ -70,23 +47,19 @@ extern "C" { #define SC_BLANK ' ' -/* audio.c */ +/* audio.cpp */ void cio_beep(void); void cio_tone(double freq, double length); -/* crtput.c */ +/* crtput.cpp */ void cio_put_c(int c); void cio_put_cc(int c, int count); void cio_put_str(const char *str); -/* kbget.c */ +/* kbget.cpp */ unsigned int cio_kb_getc(void); void cio_kb_ungetc(unsigned int c); char *cio_kb_gets(char *str, int maxlen); unsigned int cio_kb_waitc(const char *astr, int beep); -#ifdef __cplusplus -} -#endif /* _cplusplus */ - #endif diff --git a/include/ezplot.h b/include/ezplot.h index 4494f01..744a33d 100644 --- a/include/ezplot.h +++ b/include/ezplot.h @@ -2,17 +2,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ezplot.h,v 1.4 2000/05/24 22:48:17 kevin Exp $ -** $Log: ezplot.h,v $ -** Revision 1.4 2000/05/24 22:48:17 kevin -** First functional version of SDF library for X-window -** -** Revision 1.3 2000/05/07 12:46:19 kevin -** made c++ compatible -** -** Revision 1.2 2000/04/28 14:14:16 kevin -** *** empty log message *** -** +** $Id: ezplot.h,v 1.5 2000/06/13 16:20:31 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 @@ -27,6 +17,7 @@ ** along with this program; if not, write to the Free Software ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ******************************************************************************/ + /*----------------------------------------------------------------------*/ /* EZPLOT */ /* */ @@ -35,10 +26,6 @@ #ifndef __H_EZPLOT #define __H_EZPLOT -#ifdef __cplusplus -extern "C" { -#endif /* __cplusplus */ - #include #include @@ -328,27 +315,24 @@ extern struct ezplot_var ez; extern bool ezplot_firstcall; /* set to false on first call to EZSET or EZPLOT */ -/* axis.c */ -int axis_scale (double min, double max, int nint, double *minp, double *maxp, int *nintp); - -/* ezplot.c */ +/* ezplot.cpp */ SGP_ID ezplot(double x[], double y[], int num); void ezinit(void); void ezfree(void); void ezclear(void); -/* ezplot1d.c */ -void ezplot_1d(double *y, int n); - -/* ezset.c */ +/* ezset.cpp */ int ezset(char *command); -/* makefmt.c */ +/* ezplot1d.cpp */ +void ezplot_1d(double *y, int n); + +/* makefmt.cpp */ void make_numfmt(char *fmtstr, int *fldwid, int *nfrac, double min, double max, int nint); -#ifdef __cplusplus -} -#endif /* __cplusplus */ +/* axis.cpp */ +int axis_scale (double min, double max, int nint, double *minp, double *maxp, int *nintp); + #endif diff --git a/include/ir.h b/include/ir.h index 9c8f92e..ffa83fb 100644 --- a/include/ir.h +++ b/include/ir.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ir.h,v 1.25 2000/06/10 22:33:11 kevin Exp $ +** $Id: ir.h,v 1.26 2000/06/13 16:20:31 kevin Exp $ ** ** ** This program is free software; you can redistribute it and/or modify @@ -30,16 +30,6 @@ #define IR_H -struct histo_st { - int *b; /* Histogram array (# of elements in each bin) */ - int nbin; /* Number of histogram bins */ - double xmin, xmax, xinc; /* Limits of histogram boundaries */ -}; - -typedef struct histo_st HISTOGRAM; - -/*---------------------------------------------------------------------------*/ - static const int POINTS_PER_CIRCLE=90; #define MAXREMARK 99 @@ -51,7 +41,10 @@ typedef enum { SEGMENT } PElmType; -struct pelm_st { + +class PhmElement +{ + public: PElmType type; /* pelm type (box, ellipse, etc) */ double atten; /* X-ray attenuation coefficient */ double cx,cy; /* center of pelm */ @@ -63,9 +56,9 @@ struct pelm_st { double radius; /* " */ GRFMTX_2D p_to_o; /* map from phantom to standard pelm coords */ GRFMTX_2D o_to_p; /* map from std pelm coords to phantom coords */ - struct pelm_st *next; /* pointer to next pelm in phantom */ + class PhmElement *next; /* pointer to next pelm in phantom */ }; -typedef struct pelm_st PELM; +typedef class PhmElement PELM; typedef enum { P_PELMS, /* Phantom made of Pelms */ @@ -73,14 +66,18 @@ typedef enum { P_FILTER /* defined only by this type */ } PhmType; -struct phm_st { /* Phantom structure */ + +/* Phantom class */ +class Phantom +{ +public: PELM *pelm_list; /* pelm linked-list */ PhmType type; int n_pelm; /* number of pelms in phantom */ double xmin, xmax, ymin, ymax; /* extent of pelms in pelm coordinates */ double radius; /* " " */ }; -typedef struct phm_st PHANTOM; +typedef class Phantom PHANTOM; /*----------------------------------------------------------------------*/ @@ -93,12 +90,13 @@ typedef struct phm_st PHANTOM; typedef float DETECT_TYPE; -struct detarray_st { +struct DetectorArray +{ DETECT_TYPE *detval; /* Pointer to array of values recorded by detector */ int ndet; /* Number of detectors in array */ double view_angle; /* View angle in radians */ }; -typedef struct detarray_st DETARRAY; +typedef struct DetectorArray DETARRAY; typedef enum { @@ -107,7 +105,10 @@ typedef enum { DETECTOR_EQUILINEAR } ScannerGeometry; -struct detector_st { + +class Detector +{ + public: ScannerGeometry geometry; /* Geometry of detectory */ int ndet; /* Number of detectors in array */ int nview; /* Number of rotated views */ @@ -125,14 +126,17 @@ struct detector_st { double angle; /* Starting angle */ } init; }; -typedef struct detector_st DETECTOR; +typedef class Detector DETECTOR; + -struct raysum_st { +class Projections +{ + public: int fd; /* Disk file descriptor */ int file_mode; /* Current file mode (read or write) */ int header_size; /* Size of disk file header */ int geometry; /* Geometry of scanner */ - struct detarray_st **view; /* Pointer to array of detarray_st pointers */ + struct DetectorArray **view; /* Pointer to array of detarray_st pointers */ char remark[MAXREMARK+1]; /* description of raysum data */ double calctime; /* time required to calculate raysums */ @@ -146,7 +150,8 @@ struct raysum_st { double det_inc; /* increment between detectors */ double phmlen; /* Length of PHANTOM edge (phm is square) */ }; -typedef struct raysum_st RAYSUM; +typedef class Projections RAYSUM; + /*----------------------------------------------------------------------*/ /* USER SYMBOLS */ @@ -279,21 +284,21 @@ const static int RAYSUM_TRACE_ROW_CURR_VIEW=17; const static int RAYSUM_TRACE_ROW_ATTEN=18; + /************************************************************************* * FUNCTION DECLARATIONS ************************************************************************/ -/* convolve.c */ +/* convolve.cpp */ double convolve (const double f1[], const double f2[], const double dx, const int n, const int np, const FunctionSymmetry func_type); -double convolve_both (const double f1[], const double f2[], const double dx, const int n, const int np); -/* dialogs.c */ +/* dialogs.cpp */ int phm_add_pelm_kb (PHANTOM *phm); PHANTOM *phm_select (void); int interpolation_select (void); int filter_select (double *filter_param); -/* filter.c */ +/* filter.cpp */ double *filter_generate (const FilterType filt_type, double bw, double xmin, double xmax, int n, double param, const DomainType domain, int numint); double filter_spatial_response_calc (int filt_type, double x, double bw, double param, int n); double filter_spatial_response_analytic (int filt_type, double x, double bw, double param); @@ -301,7 +306,7 @@ double filter_frequency_response (int filt_type, double u, double bw, double par double sinc (double x, double mult); double integral_abscos(double u, double w); -/* options.c */ +/* options.cpp */ int opt_set_trace(const char *optarg); const char *name_of_phantom(const int phmid); int opt_set_phantom(const char *optarg); @@ -314,7 +319,7 @@ const char *name_of_filter_domain(const DomainType domain); BackprojType opt_set_backproj(const char *optarg); const char *name_of_backproj(const BackprojType backproj); -/* phm.c */ +/* phm.cpp */ PHANTOM *phm_create(const int phmid); PHANTOM *phm_create_from_file(const char *fname); PHANTOM *phm_init(void); @@ -335,12 +340,12 @@ void phm_show(const PHANTOM *phm); void phm_draw(const PHANTOM *phm); #endif -/* phmstd.c */ +/* phmstd.cpp */ void phm_std_herman (PHANTOM *phm); void phm_std_rowland (PHANTOM *phm); void phm_std_rowland_bordered (PHANTOM *phm); -/* raycollect.c */ +/* raycollect.cpp */ int raysum_collect(RAYSUM *rs, const DETECTOR *det, const PHANTOM *phm, const int start_view, const int trace, const int unit_pulse); void rayview(const PHANTOM *phm, DETARRAY *darray, const DETECTOR *det, const double xd1, const double yd1, const double xd2, const double yd2, const double xs1, const double ys1, const double xs2, const double ys2, const int unit_pulse); double phm_ray_attenuation (const PHANTOM *phm, const double x1, const double y1, const double x2, const double y2); @@ -348,11 +353,11 @@ double pelm_ray_attenuation (PELM *pelm, const double x1, const double y1, const int pelm_clip_line (const PELM *pelm, double& x1, double& y1, double& x2, double& y2); void raysum_trace_show_param (const char *label, const char *fmt, int row, int color, ...); -/* scanner.c */ +/* scanner.cpp */ DETECTOR *detector_create(const PHANTOM *phm, const ScannerGeometry geometry, int ndet, int nview, int nsample, const double rot_anglen); void detector_free(DETECTOR *det); -/* rayio.c */ +/* rayio.cpp */ RAYSUM *raysum_create(const char *fname, const int nview, const int ndet); RAYSUM *raysum_create_from_det(const char *fname, const DETECTOR *det); RAYSUM *raysum_open(const char *filename); diff --git a/include/kmath.h b/include/kmath.h index ed8de01..5599327 100644 --- a/include/kmath.h +++ b/include/kmath.h @@ -1,36 +1,15 @@ /***************************************************************************** -** This is part of the CTSim program -** Copyright (C) 1983-2000 Kevin Rosenberg -** -** $Id: kmath.h,v 1.11 2000/06/09 01:35:33 kevin Exp $ -** $Log: kmath.h,v $ -** Revision 1.11 2000/06/09 01:35:33 kevin -** Convert MPI structure to C++ class -** -** Revision 1.10 2000/06/07 00:59:38 kevin -** added imagefiles -** -** Revision 1.9 2000/05/08 20:00:48 kevin -** ANSI C changes -** -** Revision 1.8 2000/05/07 12:46:19 kevin -** made c++ compatible -** -** Revision 1.7 2000/05/04 18:16:34 kevin -** renamed filter definitions +** FILE IDENTIFICATION ** -** Revision 1.6 2000/05/02 20:00:25 kevin -** *** empty log message *** +** Name: kmath.h +** Purpose: Header file containing definitions for numerical app +** Programmer: Kevin Rosenberg +** Date Started: Nov 84 ** -** Revision 1.5 2000/05/02 15:31:39 kevin -** code cleaning -** -** Revision 1.4 2000/04/30 19:17:35 kevin -** Set up include files for conditional INTERACTIVE_GRAPHICS -** -** Revision 1.3 2000/04/28 14:14:16 kevin -** *** empty log message *** +** This is part of the CTSim program +** Copyright (C) 1983-2000 Kevin Rosenberg ** +** $Id: kmath.h,v 1.12 2000/06/13 16:20:31 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 @@ -45,23 +24,13 @@ ** along with this program; if not, write to the Free Software ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ******************************************************************************/ -/****************************************************************************** - * - * PURPOSE - * Header file containing definitions for numerical app's - * Date Started: Nov 84 - * - *****************************************************************************/ #ifndef _H_kmath #define _H_kmath #include #include - -#ifdef __cplusplus -extern "C" { -#endif /* __cplusplus */ +#include #define PI 3.14159265358979323846 #define HALFPI 1.57079632679489661923 /* PI divided by 2 */ @@ -78,33 +47,32 @@ extern "C" { #define DEG_TO_RAD(x) (x*(PI/180.)) #define RAD_TO_DEG(x) (x*(180./PI)) - #define ASSUMEDZERO 1E-10 typedef double GRFMTX_2D[3][3]; typedef double GRFMTX_3D[4][4]; -/* clip.c */ + +template +inline T nearest (double x) +{ return (x > 0 ? static_cast(x+0.5) : static_cast(x-0.5)); } + +template +inline T clamp (T value, T upperBounds, T lowerBounds) +{ return (value >= upperBounds ? upperBounds : (value <= lowerBounds ? lowerBounds : value )); } + + +/* clip.cpp */ int clip_rect(double *x1, double *y1, double *x2, double *y2, const double rect[4]); int clip_segment(double *x1, double *y1, double *x2, double *y2, const double u, const double v); int clip_sector(double *x1, double *y1, double *x2, double *y2, const double u, const double v); int clip_circle(double *x1, double *y1, double *x2, double *y2, const double cx, const double cy, const double radius, double t1, double t2); int clip_triangle(double *x1, double *y1, double *x2, double *y2, const double u, const double v, const int clip_xaxis); -/* lnearest.c */ -long int lnearest(double x); - -/* minmax.c */ -double fmax(const double a, const double b); -void minmax_dvector(const double array[], const int pts, double *xmin, double *xmax); - -/* norm_ang.c */ +/* norm_ang.cpp */ double norm_ang(double theta); -/* simpson.c */ -double simpson(const double xmin, const double xmax, const double *y, const int np); - -/* xform.c */ +/* xform.cpp */ void indent_mtx2(GRFMTX_2D m); void xlat_mtx2(GRFMTX_2D m, const double x, const double y); void scale_mtx2(GRFMTX_2D m, const double sx, const double sy); @@ -115,8 +83,11 @@ void rotate2d(double x[], double y[], int pts, double angle); void xlat2d(double x[], double y[], int pts, double xoffset, double yoffset); void scale2d(double x[], double y[], int pts, double xfact, double yfact); -#ifdef __cplusplus -} -#endif /* __cplusplus */ +/* simpson.cpp */ +double simpson(const double xmin, const double xmax, const double *y, const int np); + +/* minmax.cpp */ +void minmax_dvector(const double array[], const int pts, double *xmin, double *xmax); + #endif diff --git a/include/kstddef.h b/include/kstddef.h index 00fbddd..dfb3543 100644 --- a/include/kstddef.h +++ b/include/kstddef.h @@ -1,50 +1,15 @@ /***************************************************************************** -** This is part of the CTSim program -** Copyright (C) 1983-2000 Kevin Rosenberg -** -** $Id: kstddef.h,v 1.15 2000/06/09 01:35:33 kevin Exp $ -** $Log: kstddef.h,v $ -** Revision 1.15 2000/06/09 01:35:33 kevin -** Convert MPI structure to C++ class -** -** Revision 1.14 2000/06/07 03:49:54 kevin -** *** empty log message *** -** -** Revision 1.13 2000/06/07 00:59:38 kevin -** added imagefiles -** -** Revision 1.12 2000/06/05 01:32:45 kevin -** Added C++ compatibility -** -** Revision 1.11 2000/06/03 06:29:08 kevin -** *** empty log message *** -** -** Revision 1.10 2000/05/16 04:33:17 kevin -** Updated documentation -** -** Revision 1.9 2000/05/11 01:04:44 kevin -** Added Microsoft Windows compatibility -** -** Revision 1.8 2000/05/08 20:00:48 kevin -** ANSI C changes +** FILE IDENTIFICATION ** -** Revision 1.7 2000/05/07 12:46:19 kevin -** made c++ compatible +** File Name: kstddef.h +** Author: Kevin Rosenberg +** Purpose: Header file containing KRL standard C definitions +** Date Started: Dec. 83 ** -** Revision 1.6 2000/05/03 19:51:41 kevin -** function renaming for phantoms and phantom elements -** -** Revision 1.5 2000/05/02 20:00:25 kevin -** *** empty log message *** -** -** Revision 1.4 2000/04/28 18:00:55 kevin -** remove unused files -** -** Revision 1.3 2000/04/28 17:38:16 kevin -** Removed unused files +** This is part of the CTSim program +** Copyright (C) 1983-2000 Kevin Rosenberg ** -** Revision 1.2 2000/04/28 14:14:16 kevin -** *** empty log message *** +** $Id: kstddef.h,v 1.16 2000/06/13 16:20:31 kevin Exp $ ** ** ** This program is free software; you can redistribute it and/or modify @@ -61,22 +26,6 @@ ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ******************************************************************************/ -/****************************************************************************** - * - * FILE IDENTIFICATION - * - * File Name: STDDEF.H - * Author: Kevin Rosenberg - * Purpose: Header file containing KRL standard C definitions - * Date Started: Dec. 83 - * - * DESCRIPTION - * Header file contains KRL standard C definitions - * - * MODIFICATION LOG - * - *****************************************************************************/ - #ifndef STDDEF_H #define STDDEF_H @@ -88,17 +37,9 @@ #define snprintf _snprintf #endif -#if !defined(bool) && !defined(__cplusplus) -typedef int bool; /* Boolean variable type */ -#endif - #define STR_MAX_LEN 255 #define STR_SIZE STR_MAX_LEN+1 -#if !defined(__cplusplus) -typedef unsigned char string[STR_SIZE]; -#endif - #include #include #include @@ -138,12 +79,6 @@ typedef unsigned char string[STR_SIZE]; #define ABS(x) ((x) < 0 ? -(x) : (x)) #define SQR(x) ((x) * (x)) -#ifndef MAX -#define MAX(a,b) ((a) > (b) ? (a) : (b)) -#endif -#ifndef MIN -#define MIN(a,b) ((a) <= (b) ? (a) : (b)) -#endif #define ISWAP(a,b) {int i; i = a; a = b; b = i;} #define CLIP(n,lb,ub) if (n < lb) n = lb; else if (n > ub) n = ub @@ -243,56 +178,25 @@ typedef unsigned char kuint8; typedef double kfloat64; #endif -#ifdef __cplusplus -extern "C" { -#endif /* __cplusplus */ - - -/* allocnum.c */ -float *alloc_float(int n); -double *alloc_double(int n); -int *alloc_int(int n); - -/* fexist.c */ -int file_exists(const char *fname); - -/* kbasename.c */ -char *kbasename(const char *filename); - -/* iclip.c */ -int iclip(int n, int lb, int ub); - -/* s_head.c */ -char *str_skip_head(const char *str, const char *charlist); +/* filefuncs.cpp */ +bool file_exists(const char* fname); +const char* fileBasename(const char* filename); +FILE *sys_fopen(const char *filename, const char *mode, const char *progname); -/* s_lower.c */ +/* strfuncs.cpp */ +const char* str_skip_head(char* const str, const char* const charlist); +char* str_skip_head(char* str, char* charlist); char *str_lower(char *s); - -/* s_rmtail.c */ char *str_wrm_tail(char *str); -char *str_rm_tail(char *str, const char *charlist); - -/* s_save.c */ -char *str_save(const char *s); - -/* s_upper.c */ +char *str_rm_tail(char *str, const char* const charlist); char *str_upper(char *str); -/* sysalloc.c */ -void *sys_alloc(const int nbytes, const char *name); - -/* syserror.c */ +/* syserror.cpp */ void sys_error(int severity, const char *msg, ...); void sys_verror(int severity, const char *msg, va_list arg); void sys_error_level(int severity); -/* sysfopen.c */ -FILE *sys_fopen(const char *filename, const char *mode, const char *progname); - -/* sysfree.c */ -void sys_free(void *ptr, const char *name); - -/* timedate.c */ +/* timedate.cpp */ DATE *td_get_date(DATE *d); TIME *td_get_time(TIME *t); double td_current_sec(void); @@ -310,7 +214,7 @@ char *td_str_cdate(DATE *d); char *td_month_name(int n); char *td_day_name(int n); -/* netorder.c */ +/* netorder.cpp */ void *strreverse (void *dest, const void *src, size_t n); int read_nint16 (kuint16 *n, int fd); int write_nint16 (kuint16 const *n, int fd); @@ -323,8 +227,4 @@ int write_nfloat64 (double const *d, int fd); void ConvertNetworkOrder (void* buffer, size_t bytes); void ConvertReverseNetworkOrder (void* buffer, size_t bytes); -#ifdef __cplusplus -} -#endif /* __cplusplus */ - #endif diff --git a/include/mpiworld.h b/include/mpiworld.h index 74042b3..2889821 100644 --- a/include/mpiworld.h +++ b/include/mpiworld.h @@ -60,7 +60,7 @@ class MPIWorld int getMyLocalWorkUnits (void) const { return m_vLocalWorkUnits [m_myRank]; } - MPI::Intracomm getComm() + MPI::Intracomm& getComm() { return m_comm; } private: @@ -72,4 +72,3 @@ private: MPI::Intracomm m_comm; }; - diff --git a/include/pol.h b/include/pol.h index 2d89324..49a45c2 100644 --- a/include/pol.h +++ b/include/pol.h @@ -2,14 +2,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: pol.h,v 1.3 2000/05/07 12:46:19 kevin Exp $ -** $Log: pol.h,v $ -** Revision 1.3 2000/05/07 12:46:19 kevin -** made c++ compatible -** -** Revision 1.2 2000/04/28 14:14:16 kevin -** *** empty log message *** -** +** $Id: pol.h,v 1.4 2000/06/13 16:20:31 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 @@ -27,10 +20,6 @@ #ifndef __H_POL #define __H_POL -#ifdef __cplusplus -extern "C" { -#endif /* __cplusplus */ - /* codes for pol_usefile */ #define P_USE_STR 1 /* use string as input source */ @@ -119,8 +108,4 @@ void pol_ungetch(int c); int get_inputline(FILE *fp); void set_inputline(char *line); -#ifdef __cplusplus -} -#endif /* __cplusplus */ - #endif diff --git a/include/sgp.h b/include/sgp.h index 3267139..7189c96 100644 --- a/include/sgp.h +++ b/include/sgp.h @@ -2,26 +2,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: sgp.h,v 1.7 2000/05/24 22:48:17 kevin Exp $ -** $Log: sgp.h,v $ -** Revision 1.7 2000/05/24 22:48:17 kevin -** First functional version of SDF library for X-window -** -** Revision 1.6 2000/05/11 14:07:00 kevin -** Added support for Windows NT -** -** Revision 1.5 2000/05/07 12:46:19 kevin -** made c++ compatible -** -** Revision 1.4 2000/04/30 19:17:35 kevin -** Set up include files for conditional INTERACTIVE_GRAPHICS -** -** Revision 1.3 2000/04/28 18:35:21 kevin -** removed unused files -** -** Revision 1.2 2000/04/28 14:14:16 kevin -** *** empty log message *** -** +** $Id: sgp.h,v 1.8 2000/06/13 16:20:31 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 @@ -56,11 +37,6 @@ #include "g2_X11.h" #endif -#ifdef __cplusplus -extern "C" { -#endif /* __cplusplus */ - - /* device names */ #define CRTDEV 1 @@ -193,11 +169,11 @@ struct point {double x, y, z;}; /*-------------------------------------------------------------------------*/ -/* circle.c */ +/* circle.cpp */ void sgp2_draw_circle(const double r); void sgp2_draw_arc(double start, double stop, const double r); -/* ctm.c */ +/* ctm.cpp */ void ctm_xlat_pre_2(double x, double y); void ctm_xlat_post_2(double x, double y); void ctm_scale_pre_2(double sx, double sy); @@ -215,10 +191,10 @@ void mult_gmtx_2(GRFMTX_2D a, GRFMTX_2D b, GRFMTX_2D c); void invert_gmtx_2(GRFMTX_2D a, GRFMTX_2D b); double determ_gmtx_2(GRFMTX_2D a); -/* drawbox.c */ +/* drawbox.cpp */ void sgp2_draw_rect (double xmin, double ymin, double xmax, double ymax); -/* sgp.c */ +/* sgp.cpp */ SGP_ID sgp2_init (int xsize, int ysize, const char *title); void sgp2_close (SGP_ID gid); void sgp2_set_active_win (SGP_ID); @@ -251,7 +227,7 @@ void ctm_set_2(GRFMTX_2D m); void ctm_pre_mult_2(GRFMTX_2D m); void ctm_post_mult_2(GRFMTX_2D m); -/* sgpdrive.c */ +/* sgpdrive.cpp */ int _sgp2_init_dev(SGP_ID gid); int initdevice(int dev, int mode, int xsize, int ysize); int opendevice(int dev); @@ -275,15 +251,11 @@ void _sgp2_dev_text(SGP_ID gid, char *message); void termgrf2(void); void flushdevice(int dev); -/* sgptext.c */ +/* sgptext.cpp */ void wrtsymbol(int sym, int x, int y, DEVICE *dev); void wrtchar(int ch, int x, int y, CHARSPEC *cspec, DEVICE *dev); void wrttext(char txtstr[], int x, int y, CHARSPEC *cspec, DEVICE *dev); void crtcolor(int mode, int *f, int *b); -#ifdef __cplusplus -} -#endif /* __cplusplus */ - #endif diff --git a/src/Makefile.am b/src/Makefile.am index 8420182..7fda6a4 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -37,13 +37,13 @@ if USE_LAM CC_LAM = $(lamdir)/bin/balky LAM_EXTRA_SRC = mpiworld.cpp -ctrec-lam: ctrec.cpp +ctrec-lam: ctrec.cpp mpiworld.cpp $(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI ctrec.cpp -o ctrec-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@ -phm2rs-lam: phm2rs.cpp +phm2rs-lam: phm2rs.cpp mpiworld.cpp $(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI phm2rs.cpp -o phm2rs-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@ -phm2if-lam: phm2if.cpp +phm2if-lam: phm2if.cpp mpiworld.cpp $(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI phm2if.cpp -o phm2if-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@ endif diff --git a/src/ctrec.cpp b/src/ctrec.cpp index 19278c9..c9c3ead 100644 --- a/src/ctrec.cpp +++ b/src/ctrec.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ctrec.cpp,v 1.7 2000/06/10 22:33:11 kevin Exp $ +** $Id: ctrec.cpp,v 1.8 2000/06/13 16:20:31 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 @@ -26,6 +26,8 @@ ******************************************************************************/ #include "ct.h" +#include "timer.h" + enum {O_INTERP, O_FILTER, O_FILTER_PARAM, O_BACKPROJ, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION}; @@ -43,10 +45,11 @@ static struct option my_options[] = {0, 0, 0, 0} }; + void ctrec_usage (const char *program) { - cout << "usage: " << kbasename(program) << " raysum-file image-file nx-image ny-image [OPTIONS]" << endl; + 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; @@ -92,19 +95,19 @@ ctrec_usage (const char *program) #ifdef HAVE_MPI -static void mpi_scatter_rs (MPIWorld& mpiWorld, RAYSUM *rs_global, RAYSUM *rs_local, const int debug); +static void ScatterProjectionsMPI (MPIWorld& mpiWorld, RAYSUM *rsGlobal, RAYSUM *rsLocal, const int debug); +static void ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* imGlobal); #endif int ctrec_main (int argc, char * argv[]) { - ImageFile *im_global = NULL; - RAYSUM *rs_global = NULL; + ImageFile *imGlobal = NULL; + RAYSUM *rsGlobal = NULL; char *rs_name, *im_filename = NULL; char remark[MAXREMARK]; char filt_name[80]; - double time_start = 0, time_end = 0; char *endptr; int opt_verbose = 0; int opt_debug = 0; @@ -116,19 +119,14 @@ ctrec_main (int argc, char * argv[]) BackprojType opt_backproj = O_BPROJ_DIFF2; int nx, ny; #ifdef HAVE_MPI - ImageFile *im_local; - RAYSUM *rs_local; + ImageFile *imLocal; + RAYSUM *rsLocal; int mpi_nview, mpi_ndet; double mpi_detinc, mpi_rotinc, mpi_phmlen; - double mpi_t1, mpi_t2, mpi_t, mpi_t_g; MPIWorld mpiWorld (argc, argv); #endif -#ifdef HAVE_MPI - time_start = MPI::Wtime(); -#else - time_start = td_current_sec(); -#endif + Timer timerProgram; #ifdef HAVE_MPI if (mpiWorld.getRank() == 0) { @@ -217,26 +215,26 @@ ctrec_main (int argc, char * argv[]) nx, ny, filt_name, name_of_interpolation (opt_interp), name_of_backproj(opt_backproj)); if (opt_verbose) - fprintf (stdout, "%s\n", remark); + cout << "Remark: " << remark << endl; #ifdef HAVE_MPI } #endif #ifdef HAVE_MPI if (mpiWorld.getRank() == 0) { - rs_global = raysum_open (rs_name); - raysum_read (rs_global); + rsGlobal = raysum_open (rs_name); + raysum_read (rsGlobal); if (opt_verbose) - raysum_print_info(rs_global); + raysum_print_info(rsGlobal); - mpi_ndet = rs_global->ndet; - mpi_nview = rs_global->nview; - mpi_detinc = rs_global->det_inc; - mpi_phmlen = rs_global->phmlen; - mpi_rotinc = rs_global->rot_inc; + mpi_ndet = rsGlobal->ndet; + mpi_nview = rsGlobal->nview; + mpi_detinc = rsGlobal->det_inc; + mpi_phmlen = rsGlobal->phmlen; + mpi_rotinc = rsGlobal->rot_inc; } - mpi_t1 = MPI::Wtime(); + TimerCollectiveMPI timerBcast (mpiWorld.getComm()); 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); @@ -252,144 +250,119 @@ ctrec_main (int argc, char * argv[]) 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 (opt_verbose) { - mpi_t2 = MPI::Wtime(); - mpi_t = mpi_t2 - mpi_t1; - mpiWorld.getComm().Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0); - if (mpiWorld.getRank() == 0) - printf("Time to Bcast vars = %f secs, Max time = %f\n", mpi_t, mpi_t_g); - } + if (opt_verbose) + timerBcast.timerEndAndReport ("Time to broadcast variables"); mpiWorld.setTotalWorkUnits (mpi_nview); - rs_local = raysum_create (NULL, mpiWorld.getMyLocalWorkUnits(), mpi_ndet); - rs_local->ndet = mpi_ndet; - rs_local->nview = mpi_nview; - rs_local->det_inc = mpi_detinc; - rs_local->phmlen = mpi_phmlen; - rs_local->rot_inc = mpi_rotinc; + rsLocal = raysum_create (NULL, mpiWorld.getMyLocalWorkUnits(), mpi_ndet); + rsLocal->ndet = mpi_ndet; + rsLocal->nview = mpi_nview; + rsLocal->det_inc = mpi_detinc; + rsLocal->phmlen = mpi_phmlen; + rsLocal->rot_inc = mpi_rotinc; + TimerCollectiveMPI timerScatter (mpiWorld.getComm()); + ScatterProjectionsMPI (mpiWorld, rsGlobal, rsLocal, opt_debug); if (opt_verbose) - mpi_t1 = MPI::Wtime(); - mpi_scatter_rs(mpiWorld, rs_global, rs_local, opt_debug); - if (opt_verbose) { - mpi_t2 = MPI::Wtime(); - mpi_t = mpi_t2 - mpi_t1; - mpiWorld.getComm().Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0); - if (mpiWorld.getRank() == 0) - printf("Time to scatter rs = %f secs, Max time = %f sec\n", mpi_t, mpi_t_g); - } + timerScatter.timerEndAndReport ("Time to scatter projections"); if (mpiWorld.getRank() == 0) { - im_global = new ImageFile (im_filename, nx, ny); - im_global->fileCreate(); + imGlobal = new ImageFile (im_filename, nx, ny); + imGlobal->fileCreate(); } - im_local = new ImageFile (nx, ny); + imLocal = new ImageFile (nx, ny); #else - rs_global = raysum_open (rs_name); - raysum_read (rs_global); + rsGlobal = raysum_open (rs_name); + raysum_read (rsGlobal); if (opt_verbose) - raysum_print_info(rs_global); + raysum_print_info(rsGlobal); - im_global = new ImageFile (im_filename, nx, ny); - im_global->fileCreate(); -#endif - -#ifdef HAVE_MPI - mpi_t1 = MPI::Wtime(); - proj_reconst (*im_local, rs_local, opt_filter, opt_filter_param, - opt_interp, opt_interp_param, opt_backproj, opt_trace); - - mpi_t2 = MPI::Wtime(); - mpi_t = mpi_t2 - mpi_t1; - mpiWorld.getComm().Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0); - if (mpiWorld.getRank() == 0 && opt_verbose) - printf("Time to reconstruct = %f, Max time = %f\n", mpi_t, mpi_t_g); -#else - proj_reconst (*im_global, rs_global, opt_filter, opt_filter_param, - opt_interp, opt_interp_param, opt_backproj, opt_trace); + imGlobal = new ImageFile (im_filename, nx, ny); + imGlobal->fileCreate(); #endif #ifdef HAVE_MPI + TimerCollectiveMPI timerReconstruct (mpiWorld.getComm()); + proj_reconst (*imLocal, rsLocal, opt_filter, opt_filter_param, opt_interp, opt_interp_param, opt_backproj, opt_trace); if (opt_verbose) - mpi_t1 = MPI::Wtime(); - - int nxLocal = im_local->nx(); - int nyLocal = im_local->ny(); - ImageFileArray vLocal = im_local->getArray(); - ImageFileArray vGlobal = NULL; - if (mpiWorld.getRank() == 0) - vGlobal = im_global->getArray(); - - for (int ix = 0; ix < nxLocal; ix++) { - void *recvbuf = NULL; - if (mpiWorld.getRank() == 0) - recvbuf = vGlobal[ix]; + timerReconstruct.timerEndAndReport ("Time to reconstruct"); - mpiWorld.getComm().Reduce(vLocal[ix], recvbuf, nyLocal, im_local->getMPIDataType(), MPI::SUM, 0); - } - - if (opt_verbose) { - mpi_t2 = MPI::Wtime(); - mpi_t = mpi_t2 - mpi_t1; - mpiWorld.getComm().Reduce (&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0); - if (mpiWorld.getRank() == 0) - printf("Time to reduce image = %f secs, max time = %f\n", mpi_t, mpi_t_g); - } - - if (mpiWorld.getRank() == 0) - time_end = MPI::Wtime(); + TimerCollectiveMPI timerReduce (mpiWorld.getComm()); + ReduceImageMPI (mpiWorld, imLocal, imGlobal); + if (opt_verbose) + timerReduce.timerEndAndReport ("Time to reduce image"); #else - time_end = td_current_sec(); + proj_reconst (*imGlobal, rsGlobal, opt_filter, opt_filter_param, opt_interp, opt_interp_param, opt_backproj, opt_trace); #endif #ifdef HAVE_MPI if (mpiWorld.getRank() == 0) #endif { - raysum_close (rs_global); - double calctime = time_end - time_start; - im_global->arrayDataWrite (); - im_global->labelAdd (Array2dFileLabel::L_HISTORY, rs_global->remark, rs_global->calctime); - im_global->labelAdd (Array2dFileLabel::L_HISTORY, remark, calctime); - im_global->fileClose (); + raysum_close (rsGlobal); + double calctime = timerProgram.timerEnd(); + imGlobal->arrayDataWrite (); + imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, rsGlobal->remark, rsGlobal->calctime); + imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, remark, calctime); + imGlobal->fileClose (); if (opt_verbose) - cout << "Time active = " << calctime << " sec" << endl; + cout << "Run time: " << calctime << " seconds" << endl; } - #ifdef HAVE_MPI - MPI::Finalize(); + MPI::Finalize(); #endif - return (0); + return (0); } + +////////////////////////////////////////////////////////////////////////////////////// +// MPI Support Routines +// +////////////////////////////////////////////////////////////////////////////////////// + #ifdef HAVE_MPI -static void mpi_scatter_rs (MPIWorld& mpiWorld, RAYSUM *rs_global, RAYSUM *rs_local, const int opt_debug) +static void ScatterProjectionsMPI (MPIWorld& mpiWorld, RAYSUM *rsGlobal, RAYSUM *rsLocal, const int opt_debug) { if (mpiWorld.getRank() == 0) { for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) { for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) { - mpiWorld.getComm().Send(&rs_global->view[iw]->ndet, 1, MPI::INT, iProc, 0); - mpiWorld.getComm().Send(&rs_global->view[iw]->view_angle, 1, MPI::DOUBLE, iProc, 0); - mpiWorld.getComm().Send(rs_global->view[iw]->detval, rs_global->ndet, MPI::FLOAT, iProc, 0); + mpiWorld.getComm().Send(&rsGlobal->view[iw]->ndet, 1, MPI::INT, iProc, 0); + mpiWorld.getComm().Send(&rsGlobal->view[iw]->view_angle, 1, MPI::DOUBLE, iProc, 0); + mpiWorld.getComm().Send(rsGlobal->view[iw]->detval, rsGlobal->ndet, MPI::FLOAT, iProc, 0); } } } for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++) { MPI::Status status; + mpiWorld.getComm().Recv(&rsLocal->view[iw]->ndet, 1, MPI::INT, 0, 0, status); + mpiWorld.getComm().Recv(&rsLocal->view[iw]->view_angle, 1, MPI::DOUBLE, 0, 0, status); + mpiWorld.getComm().Recv(rsLocal->view[iw]->detval, rsLocal->ndet, MPI::FLOAT, 0, 0, status); + } + rsLocal->nview = mpiWorld.getMyLocalWorkUnits(); +} + +static void +ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* imGlobal) +{ + ImageFileArray vLocal = imLocal->getArray(); - mpiWorld.getComm().Recv(&rs_local->view[iw]->ndet, 1, MPI::INT, 0, 0, status); - mpiWorld.getComm().Recv(&rs_local->view[iw]->view_angle, 1, MPI::DOUBLE, 0, 0, status); - mpiWorld.getComm().Recv(rs_local->view[iw]->detval, rs_local->ndet, MPI::FLOAT, 0, 0, status); + 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); } - rs_local->nview = mpiWorld.getMyLocalWorkUnits(); } #endif + #ifndef NO_MAIN int main (int argc, char* argv[]) diff --git a/src/if-1.cpp b/src/if-1.cpp index e2993a5..3150cf9 100644 --- a/src/if-1.cpp +++ b/src/if-1.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: if-1.cpp,v 1.5 2000/06/09 11:03:08 kevin Exp $ +** $Id: if-1.cpp,v 1.6 2000/06/13 16:20:31 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 @@ -49,7 +49,7 @@ static struct option my_options[] = void if1_usage (const char *program) { - fprintf(stdout, "if1_usage: %s infile outfile [OPTIONS]\n", kbasename(program)); + fprintf(stdout, "if1_usage: %s infile outfile [OPTIONS]\n", fileBasename(program)); fprintf(stdout, "Generate a IF file from a IF file\n"); fprintf(stdout, "\n"); fprintf(stdout, " --invert Invert image\n"); diff --git a/src/if-2.cpp b/src/if-2.cpp index 51fc23d..29aef01 100644 --- a/src/if-2.cpp +++ b/src/if-2.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: if-2.cpp,v 1.3 2000/06/09 11:03:08 kevin Exp $ +** $Id: if-2.cpp,v 1.4 2000/06/13 16:20:31 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 @@ -48,7 +48,7 @@ static struct option my_options[] = void sdf2_usage (const char *program) { - fprintf(stdout, "sdf2_usage: %s infile1 infile2 outfile [OPTIONS]\n", kbasename(program)); + fprintf(stdout, "sdf2_usage: %s infile1 infile2 outfile [OPTIONS]\n", fileBasename(program)); fprintf(stdout, "Generate an SDF2D file from two input SDF2D files\n"); fprintf(stdout, "\n"); fprintf(stdout, " infile1 Name of first input SDF file\n"); diff --git a/src/if2img.cpp b/src/if2img.cpp index a99b801..6971fbb 100644 --- a/src/if2img.cpp +++ b/src/if2img.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: if2img.cpp,v 1.3 2000/06/09 11:03:08 kevin Exp $ +** $Id: if2img.cpp,v 1.4 2000/06/13 16:20:31 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 @@ -83,7 +83,7 @@ enum { O_FORMAT_GIF, O_FORMAT_PNG, O_FORMAT_PNG16, O_FORMAT_PGM, O_FORMAT_PGMASC void if2img_usage (const char *program) { - fprintf(stdout, "usage: %s sdfname outfile [OPTIONS]\n", kbasename(program)); + fprintf(stdout, "usage: %s sdfname outfile [OPTIONS]\n", fileBasename(program)); fprintf(stdout, "Convert IF file to an image file\n"); fprintf(stdout, "\n"); fprintf(stdout, " sdfname Name of input SDF file\n"); @@ -449,15 +449,11 @@ void sdf2d_to_pgm (ImageFile& im, char *outfile, int nxcell, int nycell, double densmin, double densmax) { FILE *fp; - int irow; - unsigned char *rowp; int nx = im.nx(); int ny = im.ny(); ImageFileArray v = im.getArray(); - rowp = (unsigned char *) sys_alloc(nx * nxcell,"sdf2d_to_img"); - if (rowp == NULL) - return; + unsigned char* rowp = new unsigned char [nx * nxcell]; if ((fp = fopen (outfile, "wb")) == NULL) return; @@ -466,33 +462,21 @@ sdf2d_to_pgm (ImageFile& im, char *outfile, int nxcell, int nycell, double densm fprintf(fp, "%d %d\n", nx, ny); fprintf(fp, "255\n"); - for (irow = ny - 1; irow >= 0; irow--) - { - int icol, ir; - - for (icol = 0; icol < nx; icol++) - { - double dens; - int p; - int pos = icol * nxcell; - dens = (v[icol][irow] - densmin) / (densmax - densmin); - if (dens < 0) - dens = 0; - else if (dens > 1) - dens = 1; - for (p = pos; p < pos + nxcell; p++) - { - rowp[p] = (unsigned int) (dens * 255.); - } - } - for (ir = 0; ir < nycell; ir++) - { - int ic; - for (ic = 0; ic < nx * nxcell; ic++) - fprintf(fp, "%c ", rowp[ic]); - } + for (int irow = ny - 1; irow >= 0; irow--) { + for (int icol = 0; icol < nx; icol++) { + int pos = icol * nxcell; + double dens = (v[icol][irow] - densmin) / (densmax - densmin); + dens = clamp (dens, 0., 1.); + for (int p = pos; p < pos + nxcell; p++) { + rowp[p] = static_cast (dens * 255.); + } + } + for (int ir = 0; ir < nycell; ir++) { + for (int ic = 0; ic < nx * nxcell; ic++) + fprintf(fp, "%c ", rowp[ic]); } - sys_free(rowp, "if2img"); + } + delete rowp; fclose(fp); } @@ -501,13 +485,11 @@ void sdf2d_to_pgmasc (ImageFile& im, char *outfile, int nxcell, int nycell, double densmin, double densmax) { FILE *fp; - int irow; - unsigned char *rowp; int nx = im.nx(); int ny = im.ny(); ImageFileArray v = im.getArray(); - rowp = (unsigned char *) sys_alloc(nx * nxcell,"sdf2d_to_img"); + unsigned char* rowp = new unsigned char [nx * nxcell]; if (rowp == NULL) return; @@ -518,34 +500,22 @@ sdf2d_to_pgmasc (ImageFile& im, char *outfile, int nxcell, int nycell, double de fprintf(fp, "%d %d\n", nx, ny); fprintf(fp, "255\n"); - for (irow = ny - 1; irow >= 0; irow--) - { - int icol, ir; - - for (icol = 0; icol < nx; icol++) - { - double dens; - int p; - int pos = icol * nxcell; - dens = (v[icol][irow] - densmin) / (densmax - densmin); - if (dens < 0) - dens = 0; - else if (dens > 1) - dens = 1; - for (p = pos; p < pos + nxcell; p++) - { - rowp[p] = (unsigned int) (dens * 255.); - } - } - for (ir = 0; ir < nycell; ir++) - { - int ic; - for (ic = 0; ic < nx * nxcell; ic++) - fprintf(fp, "%d ", rowp[ic]); - fprintf(fp, "\n"); - } + for (int irow = ny - 1; irow >= 0; irow--) { + for (int icol = 0; icol < nx; icol++) { + int pos = icol * nxcell; + double dens = (v[icol][irow] - densmin) / (densmax - densmin); + dens = clamp (dens, 0., 1.); + for (int p = pos; p < pos + nxcell; p++) { + rowp[p] = static_cast (dens * 255.); + } + } + for (int ir = 0; ir < nycell; ir++) { + for (int ic = 0; ic < nx * nxcell; ic++) + fprintf(fp, "%d ", rowp[ic]); + fprintf(fp, "\n"); } - sys_free(rowp, "if2img"); + } + delete rowp; fclose(fp); } @@ -558,16 +528,12 @@ sdf2d_to_png (ImageFile& im, char *outfile, int bitdepth, int nxcell, int nycell FILE *fp; png_structp png_ptr; png_infop info_ptr; - int irow; - unsigned char *rowp; double max_out_level = (1 << bitdepth) - 1; int nx = im.nx(); int ny = im.ny(); ImageFileArray v = im.getArray(); - rowp = (unsigned char *) sys_alloc(nx * nxcell * (bitdepth / 8),"sdf2d_to_img"); - if (rowp == NULL) - return; + unsigned char* rowp = new unsigned char [nx * nxcell * (bitdepth / 8)]; if ((fp = fopen (outfile, "wb")) == NULL) return; @@ -577,59 +543,46 @@ sdf2d_to_png (ImageFile& im, char *outfile, int bitdepth, int nxcell, int nycell return; info_ptr = png_create_info_struct(png_ptr); - if (! info_ptr) - { - png_destroy_write_struct(&png_ptr, (png_infopp) NULL); - fclose(fp); - return; - } + if (! info_ptr) { + png_destroy_write_struct(&png_ptr, (png_infopp) NULL); + fclose(fp); + return; + } - if (setjmp(png_ptr->jmpbuf)) - { - png_destroy_write_struct(&png_ptr, &info_ptr); - fclose(fp); - return; - } + if (setjmp(png_ptr->jmpbuf)) { + png_destroy_write_struct(&png_ptr, &info_ptr); + fclose(fp); + return; + } png_init_io(png_ptr, fp); png_set_IHDR(png_ptr, info_ptr, nx * nxcell, ny * nycell, bitdepth, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_DEFAULT); png_write_info(png_ptr, info_ptr); - for (irow = ny - 1; irow >= 0; irow--) - { - png_bytep row_pointer = rowp; - int icol, ir; - - for (icol = 0; icol < nx; icol++) - { - double dens; - int p; - unsigned int outval; - int pos = icol * nxcell; - dens = (v[icol][irow] - densmin) / (densmax - densmin); - if (dens < 0) - dens = 0; - else if (dens > 1) - dens = 1; - outval = (unsigned int) (dens * max_out_level); - for (p = pos; p < pos + nxcell; p++) - { - if (bitdepth == 8) - rowp[p] = outval; - else { - int rowpos = p * 2; - rowp[rowpos] = (outval >> 8) & 0xFF; - rowp[rowpos+1] = (outval & 0xFF); - } - } - } - for (ir = 0; ir < nycell; ir++) - { - png_write_rows (png_ptr, &row_pointer, 1); + for (int irow = ny - 1; irow >= 0; irow--) { + png_bytep row_pointer = rowp; + + for (int icol = 0; icol < nx; icol++) { + int pos = icol * nxcell; + double dens = (v[icol][irow] - densmin) / (densmax - densmin); + dens = clamp (dens, 0., 1.); + unsigned int outval = static_cast (dens * max_out_level); + + for (int p = pos; p < pos + nxcell; p++) { + if (bitdepth == 8) + rowp[p] = outval; + else { + int rowpos = p * 2; + rowp[rowpos] = (outval >> 8) & 0xFF; + rowp[rowpos+1] = (outval & 0xFF); } + } } - sys_free(rowp, "if2img"); + for (int ir = 0; ir < nycell; ir++) + png_write_rows (png_ptr, &row_pointer, 1); + } + delete rowp; png_write_end(png_ptr, info_ptr); png_destroy_write_struct(&png_ptr, &info_ptr); @@ -650,54 +603,39 @@ sdf2d_to_gif (ImageFile& im, char *outfile, int nxcell, int nycell, double densm gdImagePtr gif; FILE *out; int gs_indices[N_GRAYSCALE]; - int i; - int irow; - int lastrow; - unsigned char *rowp; int nx = im.nx(); int ny = im.ny(); ImageFileArray v = im.getArray(); - rowp = (unsigned char *) sys_alloc(nx * nxcell,"sdf2d_to_img"); + usnigned char* rowp = new unsigned char [nx * nxcell]; if (rowp == NULL) return; gif = gdImageCreate(nx * nxcell, ny * nycell); - for (i = 0; i < N_GRAYSCALE; i++) + for (int i = 0; i < N_GRAYSCALE; i++) gs_indices[i] = gdImageColorAllocate(gif, i, i, i); - lastrow = ny * nycell - 1; - for (irow = 0; irow < ny; irow++) - { - int icol, ir; - int rpos = irow * nycell; - for (ir = rpos; ir < rpos + nycell; ir++) - { - for (icol = 0; icol < nx; icol++) - { - double dens; - int ic; - int cpos = icol * nxcell; - dens = (v[icol][irow] - densmin) / (densmax - densmin); - if (dens < 0) - dens = 0; - else if (dens > 1) - dens = 1; - for (ic = cpos; ic < cpos + nxcell; ic++) - { - rowp[ic] = (unsigned int) (dens * (double) (N_GRAYSCALE - 1)); - gdImageSetPixel(gif, ic, lastrow - ir, gs_indices[rowp[ic]]); - } - } + int lastrow = ny * nycell - 1; + for (int irow = 0; irow < ny; irow++) { + int rpos = irow * nycell; + for (int ir = rpos; ir < rpos + nycell; ir++) { + for (int icol = 0; icol < nx; icol++) { + int cpos = icol * nxcell; + double dens = (v[icol][irow] - densmin) / (densmax - densmin); + dens = clamp(dens, 0., 1.); + for (int ic = cpos; ic < cpos + nxcell; ic++) { + rowp[ic] = (unsigned int) (dens * (double) (N_GRAYSCALE - 1)); + gdImageSetPixel(gif, ic, lastrow - ir, gs_indices[rowp[ic]]); } + } } - sys_free(rowp, "if2img"); + } + delete rowp; - if ((out = fopen(outfile,"w")) == NULL) - { - fprintf(stderr,"Error opening output file %s for writing\n", outfile); - exit(1); - } + if ((out = fopen(outfile,"w")) == NULL) { + fprintf(stderr,"Error opening output file %s for writing\n", outfile); + exit(1); + } gdImageGif(gif,out); fclose(out); gdImageDestroy(gif); diff --git a/src/ifinfo.cpp b/src/ifinfo.cpp index a78a085..7c771f7 100644 --- a/src/ifinfo.cpp +++ b/src/ifinfo.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ifinfo.cpp,v 1.5 2000/06/09 11:03:08 kevin Exp $ +** $Id: ifinfo.cpp,v 1.6 2000/06/13 16:20:31 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 @@ -49,7 +49,7 @@ static struct option my_options[] = void ifinfo_usage (const char *program) { - fprintf(stdout, "usage: %s infile [OPTIONS]\n", kbasename(program)); + fprintf(stdout, "usage: %s infile [OPTIONS]\n", fileBasename(program)); fprintf(stdout, "Imagefile information\n"); fprintf(stdout, "\n"); fprintf(stdout, " infile Name of input SDF file\n"); diff --git a/src/mpitest.c b/src/mpitest.c deleted file mode 100644 index ce118de..0000000 --- a/src/mpitest.c +++ /dev/null @@ -1,21 +0,0 @@ -#include "ct.h" - -int main(const int argc, char * const argv[]) -{ - int i; - -#ifdef MPI_CT - MPI_Init(&argc, (char ***) &argv); -#endif - - for (i = 0; i < argc; i++) - { - printf("%d: %s\n", i, argv[i]); - } - -#ifdef MPI_CT - MPI_Finalize(); -#endif - return(0); - -} diff --git a/src/phm2if.cpp b/src/phm2if.cpp index 63d5a49..208de9a 100644 --- a/src/phm2if.cpp +++ b/src/phm2if.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: phm2if.cpp,v 1.7 2000/06/09 11:03:08 kevin Exp $ +** $Id: phm2if.cpp,v 1.8 2000/06/13 16:20:31 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 @@ -25,15 +25,9 @@ ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA ******************************************************************************/ -/* FILE - * phm2sdf.c Generate a SDF image from a phantom - * - * HISTORY - * 1984 - Final DOS verion - * 1999 - First UNIX version - */ - #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 }; @@ -57,59 +51,59 @@ static struct option my_options[] = }; void -phm2sdf_usage (const char *program) +phm2if_usage (const char *program) { - fprintf(stdout,"phm2sdf_usage: %s outfile nx ny [--phantom phantom-name] [--phmfile filename] [--filter filter-name] [OPTIONS]\n",kbasename(program)); - fprintf(stdout,"Generate phantom image from a predefined --phantom or a --phmfile\n"); - fprintf(stdout,"\n"); - fprintf(stdout," outfile Name of output file for image\n"); - fprintf(stdout," nx Number of pixels X-axis\n"); - fprintf(stdout," ny Number of pixels Y-axis\n"); - fprintf(stdout," --phantom Phantom to use for projection\n"); - fprintf(stdout," herman Herman head phantom\n"); - fprintf(stdout," rowland Rowland head phantom\n"); - fprintf(stdout," browland Bordered Rowland head phantom\n"); - fprintf(stdout," unitpulse Unit pulse phantom\n"); - fprintf(stdout," --phmfile Generate Phantom from phantom file\n"); - fprintf(stdout," --filter Generate Phantom from a filter function\n"); - fprintf(stdout," abs_bandlimit Abs * Bandlimiting\n"); - fprintf(stdout," abs_sinc Abs * Sinc\n"); - fprintf(stdout," abs_cos Abs * Cosine\n"); - fprintf(stdout," abs_hamming Abs * Hamming\n"); - fprintf(stdout," shepp Shepp-Logan\n"); - fprintf(stdout," bandlimit Bandlimiting\n"); - fprintf(stdout," sinc Sinc\n"); - fprintf(stdout," cos Cosine\n"); - fprintf(stdout," triangle Triangle\n"); - fprintf(stdout," hamming Hamming\n"); - fprintf(stdout," --filter-param Alpha level for Hamming filter\n"); - fprintf(stdout," --filter-domain Set domain of filter\n"); - fprintf(stdout," spatial Spatial domain (default)\n"); - fprintf(stdout," freq Frequency domain\n"); - fprintf(stdout," --filter-bw Filter bandwidth (default = 1)\n"); - fprintf(stdout," --desc Description of raysum\n"); - fprintf(stdout," --nsample Number of samples per axis per pixel (default = 1)\n"); - fprintf(stdout," --trace Trace level to use\n"); - fprintf(stdout," none No tracing (default)\n"); - fprintf(stdout," text Trace text level\n"); - fprintf(stdout," phm Trace phantom\n"); - fprintf(stdout," rays Trace rays\n"); - fprintf(stdout," plot Trace plot\n"); - fprintf(stdout," clipping Trace clipping\n"); - fprintf(stdout," --debug Debug mode\n"); - fprintf(stdout," --verbose Verbose mode\n"); - fprintf(stdout," --version Print version\n"); - fprintf(stdout," --help Print this help message\n"); + 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 << " 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* im_global, ImageFile* im_local, const int opt_debug); +void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* imGlobal, ImageFile* imLocal, const int opt_debug); #endif int -phm2sdf_main (int argc, char* argv[]) +phm2if_main (int argc, char* argv[]) { - ImageFile* im_global = NULL; + ImageFile* imGlobal = NULL; PHANTOM *phm = NULL; int opt_nx = 0, opt_ny = 0; int opt_nsample = 1; @@ -125,23 +119,16 @@ phm2sdf_main (int argc, char* argv[]) double opt_filter_bw = 1.; int opt_trace = TRACE_NONE; int opt_verbose = 0; - double time_start=0, time_end=0; #ifdef HAVE_MPI - double mpi_t1, mpi_t2, mpi_t, mpi_t_g; - ImageFile* im_local = NULL; + ImageFile* imLocal = NULL; MPIWorld mpiWorld (argc, argv); #endif -#ifdef HAVE_MPI - time_start = MPI::Wtime(); -#else - time_start = td_current_sec(); -#endif + Timer timerProgram; #ifdef HAVE_MPI if (mpiWorld.getRank() == 0) { #endif - strcpy(opt_desc, ""); strcpy(opt_phmfilename, ""); while (1) { @@ -155,7 +142,7 @@ phm2sdf_main (int argc, char* argv[]) switch (c) { case O_PHANTOM: if ((opt_phmnum = opt_set_phantom(optarg)) < 0) { - phm2sdf_usage(argv[0]); + phm2if_usage(argv[0]); return (1); } break; @@ -164,7 +151,7 @@ phm2sdf_main (int argc, char* argv[]) phm = phm_create_from_file(opt_phmfilename); #ifdef HAVE_MPI if (mpiWorld.getRank() == 0) - fprintf(stderr, "Can't use phantom from file in MPI mode\n"); + cerr << "Can't use phantom from file in MPI mode" << endl; return (1); #endif break; @@ -176,19 +163,19 @@ phm2sdf_main (int argc, char* argv[]) break; case O_TRACE: if ((opt_trace = opt_set_trace(optarg)) < 0) { - phm2sdf_usage(argv[0]); + phm2if_usage(argv[0]); return (1); } break; case O_FILTER: if ((opt_filter = opt_set_filter(optarg)) < 0) { - phm2sdf_usage (argv[0]); + phm2if_usage (argv[0]); return (1); } break; case O_FILTER_DOMAIN: if ((opt_filter_domain = opt_set_filter_domain(optarg)) < 0) { - phm2sdf_usage (argv[0]); + phm2if_usage (argv[0]); return (1); } break; @@ -200,7 +187,7 @@ phm2sdf_main (int argc, char* argv[]) endstr = optarg + strlen(optarg); if (endptr != endstr) { fprintf(stderr,"Error setting --filter-bw to %s\n", optarg); - phm2sdf_usage(argv[0]); + phm2if_usage(argv[0]); return (1); } break; @@ -209,7 +196,7 @@ phm2sdf_main (int argc, char* argv[]) endstr = optarg + strlen(optarg); if (endptr != endstr) { fprintf(stderr,"Error setting --filter-param to %s\n", optarg); - phm2sdf_usage(argv[0]); + phm2if_usage(argv[0]); return (1); } break; @@ -218,34 +205,34 @@ phm2sdf_main (int argc, char* argv[]) endstr = optarg + strlen(optarg); if (endptr != endstr) { fprintf(stderr,"Error setting --nsample to %s\n", optarg); - phm2sdf_usage(argv[0]); + phm2if_usage(argv[0]); return (1); } break; case O_VERSION: #ifdef VERSION - fprintf(stdout, "Version %s\n", VERSION); + cout << "Version " << VERSION << endl; #else - fprintf(stderr, "Unknown version number"); + cerr << "Unknown version number" << endl; #endif case O_HELP: case '?': - phm2sdf_usage(argv[0]); + phm2if_usage(argv[0]); return (0); default: - phm2sdf_usage(argv[0]); + phm2if_usage(argv[0]); return (1); } } if (phm == NULL && opt_phmnum == -1 && opt_filter == -1) { - fprintf(stderr, "No phantom defined\n"); - phm2sdf_usage(argv[0]); + cerr << "No phantom defined" << endl; + phm2if_usage(argv[0]); return (1); } if (optind + 3 != argc) { - phm2sdf_usage(argv[0]); + phm2if_usage(argv[0]); return (1); } opt_outfile = argv[optind]; @@ -253,14 +240,14 @@ phm2sdf_main (int argc, char* argv[]) endstr = argv[optind+1] + strlen(argv[optind+1]); if (endptr != endstr) { fprintf(stderr,"Error setting nx to %s\n", argv[optind+1]); - phm2sdf_usage(argv[0]); + 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) { fprintf(stderr,"Error setting ny to %s\n", argv[optind+2]); - phm2sdf_usage(argv[0]); + phm2if_usage(argv[0]); return (1); } @@ -286,13 +273,13 @@ phm2sdf_main (int argc, char* argv[]) strncpy(opt_desc, str, sizeof(opt_desc)); if (opt_verbose) - printf("Rasterize Phantom to Image\n\n"); + cout << "Rasterize Phantom to Image" << endl << endl; #ifdef HAVE_MPI } #endif #ifdef HAVE_MPI - mpi_t1 = MPI::Wtime(); + TimerMPI timerBcast (mpiWorld.getComm()); 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); @@ -305,24 +292,19 @@ phm2sdf_main (int argc, char* argv[]) mpiWorld.getComm().Bcast (&opt_filter_param, 1, MPI::DOUBLE, 0); mpiWorld.getComm().Bcast (&opt_filter_bw, 1, MPI::DOUBLE, 0); - if (opt_verbose) { - mpi_t2 = MPI::Wtime(); - mpi_t = mpi_t2 - mpi_t1; - mpiWorld.getComm().Reduce (&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0); - if (mpiWorld.getRank() == 0) - printf("Time to Bcast vars = %f secs, Max time = %f\n", mpi_t, mpi_t_g); - } + if (opt_verbose) + timerBcast.timerEndAndReport ("Time to broadcast variables"); mpiWorld.setTotalWorkUnits (opt_nx); if (mpiWorld.getRank() == 0) { - im_global = new ImageFile (opt_outfile, opt_nx, opt_ny); - im_global->fileCreate(); + imGlobal = new ImageFile (opt_outfile, opt_nx, opt_ny); + imGlobal->fileCreate(); } - im_local = new ImageFile (opt_nx, opt_ny); + imLocal = new ImageFile (opt_nx, opt_ny); #else - im_global = new ImageFile (opt_outfile, opt_nx, opt_ny); - im_global->fileCreate (); + imGlobal = new ImageFile (opt_outfile, opt_nx, opt_ny); + imGlobal->fileCreate (); #endif if (opt_phmnum >= 0) @@ -331,76 +313,60 @@ phm2sdf_main (int argc, char* argv[]) #ifdef HAVE_MPI else { if (mpiWorld.getRank() == 0) - fprintf(stderr, "phmnum < 0\n"); + cerr << "phmnum < 0" << endl; exit(1); } #endif - double calctime = 0; ImageFileArray v; - #ifdef HAVE_MPI if (mpiWorld.getRank() == 0) - v = im_global->getArray (); + v = imGlobal->getArray (); if (phm->type == P_UNIT_PULSE) { if (mpiWorld.getRank() == 0) { v[opt_nx/2][opt_ny/2] = 1.; - calctime = 0; } } else if (opt_filter != -1) { if (mpiWorld.getRank() == 0) { - image_filter_response (*im_global, opt_filter_domain, opt_filter_bw, opt_filter, opt_filter_param, opt_trace); - calctime = 0; + image_filter_response (*imGlobal, opt_filter_domain, opt_filter_bw, opt_filter, opt_filter_param, opt_trace); } } else { - mpiWorld.getComm().Barrier(); + TimerMPI timerRasterize (mpiWorld.getComm()); + phm_to_imagefile (phm, *imLocal, mpiWorld.getMyStartWorkUnit(), mpiWorld.getMyLocalWorkUnits(), opt_nsample, opt_trace); + if (opt_verbose) + timerRasterize.timerEndAndReport ("Time to rasterize phantom"); + TimerMPI timerGather (mpiWorld.getComm()); + mpi_gather_image (mpiWorld, imGlobal, imLocal, opt_debug); if (opt_verbose) - mpi_t1 = MPI::Wtime(); - phm_to_imagefile (phm, *im_local, mpiWorld.getMyStartWorkUnit(), mpiWorld.getMyLocalWorkUnits(), opt_nsample, opt_trace); - if (opt_verbose) { - mpi_t2 = MPI::Wtime(); - mpi_t = mpi_t2 - mpi_t1; - mpiWorld.getComm().Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0); - if (mpiWorld.getRank() == 0) - printf("Time to compile phm = %f secs, Max time = %f\n", mpi_t, mpi_t_g); - } - mpi_gather_image (mpiWorld, im_global, im_local, opt_debug); + timerGather.timerEndAndReport ("Time to gather image"); } #else - v = im_global->getArray (); + v = imGlobal->getArray (); if (phm->type == P_UNIT_PULSE) { v[opt_nx/2][opt_ny/2] = 1.; - calctime = 0; } else if (opt_filter != -1) { - image_filter_response (*im_global, opt_filter_domain, opt_filter_bw, opt_filter, opt_filter_param, opt_trace); - calctime = 0; + image_filter_response (*imGlobal, opt_filter_domain, opt_filter_bw, opt_filter, opt_filter_param, opt_trace); } else { #if HAVE_SGP if (opt_trace >= TRACE_PHM) phm_show(phm); #endif - phm_to_imagefile (phm, *im_global, 0, opt_nx, opt_nsample, opt_trace); + phm_to_imagefile (phm, *imGlobal, 0, opt_nx, opt_nsample, opt_trace); } #endif -#ifdef HAVE_MPI - time_end = MPI::Wtime(); -#else - time_end = td_current_sec(); -#endif - #ifdef HAVE_MPI if (mpiWorld.getRank() == 0) #endif { - im_global->arrayDataWrite (); - calctime = time_end - time_start; - im_global->labelAdd (Array2dFileLabel::L_HISTORY, opt_desc, calctime); - im_global->fileClose (); + imGlobal->arrayDataWrite (); + double calctime = timerProgram.timerEnd (); + imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, opt_desc, calctime); + imGlobal->fileClose (); if (opt_verbose) - fprintf (stdout, "Time to compile image = %.2f sec\n\n", calctime); + cout << "Time to rasterized phantom: " << calctime << " seconds" << endl; if (opt_trace >= TRACE_PHM) { double dmin, dmax; @@ -410,7 +376,7 @@ phm2sdf_main (int argc, char* argv[]) scanf ("%d", &nscale); printf ("Enter minimum and maximum densities (min, max): "); scanf ("%lf %lf", &dmin, &dmax); - image_display_scale (*im_global, nscale, dmin, dmax); + image_display_scale (*imGlobal, nscale, dmin, dmax); } } @@ -422,23 +388,23 @@ phm2sdf_main (int argc, char* argv[]) #ifdef HAVE_MPI -void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* im_global, ImageFile* im_local, const int opt_debug) +void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* imGlobal, ImageFile* imLocal, const int opt_debug) { - ImageFileArray vLocal = im_local->getArray(); + ImageFileArray vLocal = imLocal->getArray(); ImageFileArray vGlobal = NULL; - int nyLocal = im_local->ny(); + int nyLocal = imLocal->ny(); if (mpiWorld.getRank() == 0) - vGlobal = im_global->getArray(); + vGlobal = imGlobal->getArray(); for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++) - mpiWorld.getComm().Send(vLocal[iw], nyLocal, im_local->getMPIDataType(), 0, 0); + 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, im_local->getMPIDataType(), iProc, 0, status); + mpiWorld.getComm().Recv(vGlobal[iw], nyLocal, imLocal->getMPIDataType(), iProc, 0, status); } } } @@ -450,6 +416,6 @@ void mpi_gather_image (MPIWorld& mpiWorld, ImageFile* im_global, ImageFile* im_l int main (int argc, char* argv[]) { - return (phm2sdf_main(argc, argv)); + return (phm2if_main(argc, argv)); } #endif diff --git a/src/phm2rs.cpp b/src/phm2rs.cpp index ac30c59..f87f301 100644 --- a/src/phm2rs.cpp +++ b/src/phm2rs.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: phm2rs.cpp,v 1.4 2000/06/09 01:35:33 kevin Exp $ +** $Id: phm2rs.cpp,v 1.5 2000/06/13 16:20:31 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 @@ -25,8 +25,9 @@ ** 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 }; @@ -50,37 +51,37 @@ static struct option phm2rs_options[] = void phm2rs_usage (const char *program) { - fprintf(stdout,"usage: %s outfile ndet nview [--phantom phantom-name] [--phmfile filename] [OPTIONS]\n", kbasename(program)); - fprintf(stdout,"Calculate raysums (projections) through phantom object, either\n"); - fprintf(stdout,"a predefined --phantom or a --phmfile\n"); - fprintf(stdout,"\n"); - fprintf(stdout," outfile Name of output file for raysums\n"); - fprintf(stdout," ndet Number of detectors\n"); - fprintf(stdout," nview Number of rotated views\n"); - fprintf(stdout," --phantom Phantom to use for projection\n"); - fprintf(stdout," herman Herman head phantom\n"); - fprintf(stdout," rowland Rowland head phantom\n"); - fprintf(stdout," browland Bordered Rowland head phantom\n"); - fprintf(stdout," unitpulse Unit pulse phantom\n"); - fprintf(stdout," --phmfile Get Phantom from phantom file\n"); - fprintf(stdout," --desc Description of raysum\n"); - fprintf(stdout," --nray Number of rays per detector (default = 1)\n"); - fprintf(stdout," --rotangle Degrees to rotate view through, multiple of PI (default = 1)\n"); - fprintf(stdout," --trace Trace level to use\n"); - fprintf(stdout," none No tracing (default)\n"); - fprintf(stdout," text Trace text level\n"); - fprintf(stdout," phm Trace phantom image\n"); - fprintf(stdout," rays Trace rays\n"); - fprintf(stdout," plot Trace plot\n"); - fprintf(stdout," clipping Trace clipping\n"); - fprintf(stdout," --verbose Verbose mode\n"); - fprintf(stdout," --debug Debug mode\n"); - fprintf(stdout," --version Print version\n"); - fprintf(stdout," --help Print this help message\n"); + cout << "usage: " << fileBasename(program) << " outfile ndet nview [--phantom phantom-name] [--phmfile filename] [OPTIONS]" << endl; + cout << "Calculate raysums (projections) through phantom object, either" << endl; + cout << "a predefined --phantom or a --phmfile" << endl; + cout << "" << endl; + cout << " outfile Name of output file for raysums" << endl; + cout << " ndet Number of detectors" << endl; + cout << " nview Number of rotated views" << endl; + cout << " --phantom Phantom to use for projection" << endl; + cout << " herman Herman head phantom" << endl; + cout << " rowland Rowland head phantom" << endl; + cout << " browland Bordered Rowland head phantom" << endl; + cout << " unitpulse Unit pulse phantom" << endl; + cout << " --phmfile Get Phantom from phantom file" << endl; + cout << " --desc Description of raysum" << endl; + cout << " --nray Number of rays per detector (default = 1)" << endl; + cout << " --rotangle Degrees to rotate view through, multiple of PI (default = 1)" << endl; + cout << " --trace Trace level to use" << endl; + cout << " none No tracing (default)" << endl; + cout << " text Trace text level" << endl; + cout << " phm Trace phantom image" << endl; + cout << " rays Trace rays" << endl; + cout << " plot Trace plot" << endl; + cout << " clipping Trace clipping" << endl; + cout << " --verbose Verbose mode" << endl; + cout << " --debug Debug mode" << endl; + cout << " --version Print version" << endl; + cout << " --help Print this help message" << endl; } #ifdef HAVE_MPI -void mpi_gather_rs (MPIWorld& mpiWorld, RAYSUM *rs_global, RAYSUM *rs_local, const int opt_debug); +void GatherProjectionsMPI (MPIWorld& mpiWorld, RAYSUM *rsGlobal, RAYSUM *rsLocal, const int opt_debug); #endif int @@ -88,28 +89,25 @@ phm2rs_main (int argc, char* argv[]) { DETECTOR *det; PHANTOM *phm = NULL; - RAYSUM *rs_global = NULL; - char str[256], *opt_outfile = NULL, opt_desc[MAXREMARK+1]; + RAYSUM *rsGlobal = NULL; + char *opt_outfile = NULL; + char opt_desc[MAXREMARK+1]; char opt_phmfilename[256]; char *endptr, *endstr; int opt_ndet, opt_nview; int opt_nray = 1; - int opt_trace = 0, opt_phmnum = -1; + int opt_trace = 0; + int opt_phmnum = -1; int opt_verbose = 0; int opt_debug = 0; double opt_rotangle = 1; - double time_start = 0, time_end = 0; #ifdef HAVE_MPI - double mpi_t1 = 0, mpi_t2 = 0, mpi_t, mpi_t_g; - RAYSUM *rs_local; + RAYSUM *rsLocal; MPIWorld mpiWorld (argc, argv); - - time_start = MPI::Wtime(); -#else - time_start = td_current_sec(); #endif - + + Timer timerProgram; #ifdef HAVE_MPI if (mpiWorld.getRank() == 0) { @@ -135,7 +133,7 @@ phm2rs_main (int argc, char* argv[]) case O_PHMFILE: #ifdef HAVE_MPI if (mpiWorld.getRank() == 0) - fprintf(stderr, "Can not read phantom from file in MPI mode\n"); + cerr << "Can not read phantom from file in MPI mode" << endl; return (1); #endif strncpy(opt_phmfilename, optarg, sizeof(opt_phmfilename)); @@ -161,7 +159,7 @@ phm2rs_main (int argc, char* argv[]) opt_rotangle = strtod(optarg, &endptr); endstr = optarg + strlen(optarg); if (endptr != endstr) { - fprintf(stderr,"Error setting --rotangle to %s\n", optarg); + cerr << "Error setting --rotangle to " << optarg << endl; phm2rs_usage(argv[0]); return (1); } @@ -170,16 +168,16 @@ phm2rs_main (int argc, char* argv[]) opt_nray = strtol(optarg, &endptr, 10); endstr = optarg + strlen(optarg); if (endptr != endstr) { - fprintf(stderr,"Error setting --nray to %s\n", optarg); + cerr << "Error setting --nray to %s" << optarg << endl; phm2rs_usage(argv[0]); return (1); } break; case O_VERSION: #ifdef VERSION - fprintf(stdout, "Version %s\n", VERSION); + cout << "Version: " << VERSION << endl; #else - fprintf(stderr, "Unknown version number"); + cout << "Unknown version number" << endl; #endif exit(0); case O_HELP: @@ -193,7 +191,7 @@ phm2rs_main (int argc, char* argv[]) } if (phm == NULL) { - fprintf(stderr, "No phantom defined\n"); + cerr << "No phantom defined" << endl; phm2rs_usage(argv[0]); return (1); } @@ -206,18 +204,19 @@ phm2rs_main (int argc, char* argv[]) opt_ndet = strtol(argv[optind+1], &endptr, 10); endstr = argv[optind+1] + strlen(argv[optind+1]); if (endptr != endstr) { - fprintf(stderr,"Error setting --ndet to %s\n", argv[optind+1]); + cerr << "Error setting --ndet to " << argv[optind+1] << endl; phm2rs_usage(argv[0]); return (1); } opt_nview = strtol(argv[optind+2], &endptr, 10); endstr = argv[optind+2] + strlen(argv[optind+2]); if (endptr != endstr) { - fprintf(stderr,"Error setting --nview to %s\n", argv[optind+2]); + cerr << "Error setting --nview to " << argv[optind+2] << endl; phm2rs_usage(argv[0]); return (1); } + char str[256]; snprintf(str, sizeof(str), "Raysum_Collect: NDet=%d, Nview=%d, NRay=%d, RotAngle=%.2f, ", opt_ndet, opt_nview, opt_nray, opt_rotangle); @@ -238,7 +237,6 @@ phm2rs_main (int argc, char* argv[]) #endif #ifdef HAVE_MPI - mpiWorld.getComm().Barrier (); 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); @@ -259,66 +257,43 @@ phm2rs_main (int argc, char* argv[]) mpiWorld.setTotalWorkUnits (opt_nview); if (mpiWorld.getRank() == 0) { - rs_global = raysum_create_from_det (opt_outfile, det); - raysum_alloc_views(rs_global); + rsGlobal = raysum_create_from_det (opt_outfile, det); + raysum_alloc_views(rsGlobal); } - rs_local = raysum_create_from_det (NULL, det); - rs_local->nview = mpiWorld.getMyLocalWorkUnits(); + rsLocal = raysum_create_from_det (NULL, det); + rsLocal->nview = mpiWorld.getMyLocalWorkUnits(); if (opt_debug) - printf("rs_local->nview = %d (process %d)\n", rs_local->nview, mpiWorld.getRank()); + cout << "rsLocal->nview = " << rsLocal->nview << " (process " << mpiWorld.getRank() << ")" << endl;; + TimerMPI timerProject (mpiWorld.getComm()); + raysum_collect (rsLocal, det, phm, mpiWorld.getMyStartWorkUnit(), opt_trace, FALSE); if (opt_verbose) - mpi_t1 = MPI::Wtime(); - raysum_collect (rs_local, det, phm, mpiWorld.getMyStartWorkUnit(), opt_trace, FALSE); - if (opt_verbose) { - mpi_t2 = MPI::Wtime(); - mpi_t = mpi_t2 - mpi_t1; - mpiWorld.getComm().Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0); - if (mpiWorld.getRank() == 0) - printf("Time to collect rs = %f secs, Max = %f secs\n", mpi_t, mpi_t_g); - } + timerProject.timerEndAndReport ("Time to collect projections"); + TimerMPI timerGather (mpiWorld.getComm()); + GatherProjectionsMPI (mpiWorld, rsGlobal, rsLocal, opt_debug); if (opt_verbose) - mpi_t1 = MPI::Wtime(); - - mpi_gather_rs (mpiWorld, rs_global, rs_local, opt_debug); - - if (opt_verbose) { - mpi_t2 = MPI::Wtime(); - mpi_t = mpi_t2 - mpi_t1; - mpiWorld.getComm().Reduce(&mpi_t, &mpi_t_g, 1, MPI::DOUBLE, MPI::MAX, 0); - if (mpiWorld.getRank() == 0) - printf("Time to gather rs = %f secs, Max = %f secs\n", mpi_t, mpi_t_g); - } + timerGather.timerEndAndReport ("Time to gather projections"); #else - rs_global = raysum_create_from_det (opt_outfile, det); - raysum_collect (rs_global, det, phm, 0, opt_trace, FALSE); + rsGlobal = raysum_create_from_det (opt_outfile, det); + raysum_collect (rsGlobal, det, phm, 0, opt_trace, FALSE); #endif -#ifndef HAVE_MPI - time_end = td_current_sec(); -#else - time_end = MPI::Wtime(); +#ifdef HAVE_MPI if (mpiWorld.getRank() == 0) #endif { - rs_global->calctime = time_end - time_start; - strncpy (rs_global->remark, opt_desc, sizeof(rs_global->remark)); + rsGlobal->calctime = timerProgram.timerEnd(); + strncpy (rsGlobal->remark, opt_desc, sizeof(rsGlobal->remark)); + raysum_write (rsGlobal); + raysum_close (rsGlobal); if (opt_verbose) { - fprintf(stdout, "Remark: %s\n", rs_global->remark); - fprintf(stdout, "Time active = %.2f secs\n", rs_global->calctime); + cout << "Remark: " << rsGlobal->remark << endl; + cout << "Run time: " << rsGlobal->calctime << " seconds" << endl; } } -#ifdef HAVE_MPI - if (mpiWorld.getRank() == 0) -#endif - { - raysum_write (rs_global); - raysum_close (rs_global); - } - phm_free (phm); detector_free (det); @@ -327,19 +302,19 @@ phm2rs_main (int argc, char* argv[]) /* FUNCTION - * mpi_gather_rs + * GatherProjectionsMPI * * SYNOPSIS - * Gather's raysums from all processes in rs_local in rs_global in process 0 + * Gather's raysums from all processes in rsLocal in rsGlobal in process 0 */ #ifdef HAVE_MPI -void mpi_gather_rs (MPIWorld& mpiWorld, RAYSUM *rs_global, RAYSUM *rs_local, const int opt_debug) +void GatherProjectionsMPI (MPIWorld& mpiWorld, RAYSUM *rsGlobal, RAYSUM *rsLocal, const int opt_debug) { for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++) { - mpiWorld.getComm().Send(&rs_local->view[iw]->view_angle, 1, MPI::DOUBLE, 0, 0); - mpiWorld.getComm().Send(&rs_local->view[iw]->ndet, 1, MPI::INT, 0, 0); - mpiWorld.getComm().Send(rs_local->view[iw]->detval, rs_local->ndet, MPI::FLOAT, 0, 0); + mpiWorld.getComm().Send(&rsLocal->view[iw]->view_angle, 1, MPI::DOUBLE, 0, 0); + mpiWorld.getComm().Send(&rsLocal->view[iw]->ndet, 1, MPI::INT, 0, 0); + mpiWorld.getComm().Send(rsLocal->view[iw]->detval, rsLocal->ndet, MPI::FLOAT, 0, 0); } if (mpiWorld.getRank() == 0) { @@ -347,9 +322,9 @@ void mpi_gather_rs (MPIWorld& mpiWorld, RAYSUM *rs_global, RAYSUM *rs_local, con for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) { MPI::Status status; - mpiWorld.getComm().Recv(&rs_global->view[iw]->view_angle, 1, MPI::DOUBLE, iProc, 0, status); - mpiWorld.getComm().Recv(&rs_global->view[iw]->ndet, 1, MPI::INT, iProc, 0, status); - mpiWorld.getComm().Recv(rs_global->view[iw]->detval, rs_global->ndet, MPI::FLOAT, iProc, 0, status); + mpiWorld.getComm().Recv(&rsGlobal->view[iw]->view_angle, 1, MPI::DOUBLE, iProc, 0, status); + mpiWorld.getComm().Recv(&rsGlobal->view[iw]->ndet, 1, MPI::INT, iProc, 0, status); + mpiWorld.getComm().Recv(rsGlobal->view[iw]->detval, rsGlobal->ndet, MPI::FLOAT, iProc, 0, status); } } } diff --git a/src/rs2if.cpp b/src/rs2if.cpp index 82ed201..5da8d68 100644 --- a/src/rs2if.cpp +++ b/src/rs2if.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: rs2if.cpp,v 1.4 2000/06/09 11:03:08 kevin Exp $ +** $Id: rs2if.cpp,v 1.5 2000/06/13 16:20:31 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 @@ -49,7 +49,7 @@ static struct option my_options[] = void rs2if_usage (const char *program) { - fprintf(stdout, "usage: %s in-rs-file out-if-file [OPTIONS]\n", kbasename(program)); + fprintf(stdout, "usage: %s in-rs-file out-if-file [OPTIONS]\n", fileBasename(program)); fprintf(stdout, "This program converts a raysum file to a IF file\n"); fprintf(stdout, "\n"); fprintf(stdout, " --verbose Verbose mode\n");