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.
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
** 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 ***
**
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);
** 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
#ifndef __CIO_H
#define __CIO_H
-#ifdef __cplusplus
-extern "C" {
-#endif /* __cplusplus */
-
#define C_BLACK 0 /* color codes */
#define C_BLUE 1
#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
** 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
** along with this program; if not, write to the Free Software
** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
******************************************************************************/
+
/*----------------------------------------------------------------------*/
/* EZPLOT */
/* */
#ifndef __H_EZPLOT
#define __H_EZPLOT
-#ifdef __cplusplus
-extern "C" {
-#endif /* __cplusplus */
-
#include <stdio.h>
#include <stddef.h>
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
** 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
#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
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 */
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 */
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;
/*----------------------------------------------------------------------*/
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 {
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 */
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 */
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 */
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);
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);
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);
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);
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);
/*****************************************************************************
-** 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
** 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 <stdio.h>
#include <math.h>
-
-#ifdef __cplusplus
-extern "C" {
-#endif /* __cplusplus */
+#include <algo.h>
#define PI 3.14159265358979323846
#define HALFPI 1.57079632679489661923 /* PI divided by 2 */
#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<class T>
+inline T nearest (double x)
+{ return (x > 0 ? static_cast<T>(x+0.5) : static_cast<T>(x-0.5)); }
+
+template<class T>
+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);
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
/*****************************************************************************
-** 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
** 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
#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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#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
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);
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);
void ConvertNetworkOrder (void* buffer, size_t bytes);
void ConvertReverseNetworkOrder (void* buffer, size_t bytes);
-#ifdef __cplusplus
-}
-#endif /* __cplusplus */
-
#endif
int getMyLocalWorkUnits (void) const
{ return m_vLocalWorkUnits [m_myRank]; }
- MPI::Intracomm getComm()
+ MPI::Intracomm& getComm()
{ return m_comm; }
private:
MPI::Intracomm m_comm;
};
-
** 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
#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 */
int get_inputline(FILE *fp);
void set_inputline(char *line);
-#ifdef __cplusplus
-}
-#endif /* __cplusplus */
-
#endif
** 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
#include "g2_X11.h"
#endif
-#ifdef __cplusplus
-extern "C" {
-#endif /* __cplusplus */
-
-
/* device names */
#define CRTDEV 1
/*-------------------------------------------------------------------------*/
-/* 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);
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);
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);
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
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
** 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
******************************************************************************/
#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};
{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;
#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;
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) {
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);
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[])
** 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
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");
** 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
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");
** 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
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");
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;
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<unsigned int> (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);
}
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;
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<unsigned int> (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);
}
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;
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<unsigned int> (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);
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);
** 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
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");
+++ /dev/null
-#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);
-
-}
** 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
** 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 };
};
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;
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) {
switch (c) {
case O_PHANTOM:
if ((opt_phmnum = opt_set_phantom(optarg)) < 0) {
- phm2sdf_usage(argv[0]);
+ phm2if_usage(argv[0]);
return (1);
}
break;
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;
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;
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;
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;
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];
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);
}
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);
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)
#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;
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);
}
}
#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);
}
}
}
int
main (int argc, char* argv[])
{
- return (phm2sdf_main(argc, argv));
+ return (phm2if_main(argc, argv));
}
#endif
** 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
** 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 };
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
{
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) {
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));
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);
}
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:
}
if (phm == NULL) {
- fprintf(stderr, "No phantom defined\n");
+ cerr << "No phantom defined" << endl;
phm2rs_usage(argv[0]);
return (1);
}
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);
#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);
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);
/* 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) {
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);
}
}
}
** 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
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");