From 44ba9ce559d2d52cbd7bbea6bcd76242840fd3eb Mon Sep 17 00:00:00 2001 From: "Kevin M. Rosenberg" Date: Tue, 20 Jun 2000 17:54:51 +0000 Subject: [PATCH] r115: *** empty log message *** --- ChangeLog | 11 +- INSTALL | 10 +- README | 52 ++--- configure | 2 +- configure.in | 2 +- include/Makefile.am | 2 +- include/byteorder.h | 44 ---- include/ct.h | 17 +- include/ezplot.h | 3 +- include/fnetorderstream.h | 105 ++++++++++ include/imagefile.h | 227 ++++++++++---------- include/projections.h | 9 +- libctgraphics/ezplot.cpp | 14 +- libctgraphics/ezset.cpp | 6 +- libctsim/Makefile.am | 2 +- libctsim/convolve.cpp | 84 -------- libctsim/filter.cpp | 189 +++++++++++++---- libctsim/imagefile.cpp | 6 +- libctsim/projections.cpp | 214 ++++++++----------- libctsim/reconstr.cpp | 27 +-- libctsupport/Makefile.am | 2 +- libctsupport/byteorder.cpp | 343 ------------------------------- libctsupport/fnetorderstream.cpp | 180 ++++++++++++++++ src/sample-ctrec.sh.in | 5 + 24 files changed, 724 insertions(+), 832 deletions(-) delete mode 100644 include/byteorder.h create mode 100644 include/fnetorderstream.h delete mode 100644 libctsim/convolve.cpp delete mode 100644 libctsupport/byteorder.cpp create mode 100644 libctsupport/fnetorderstream.cpp diff --git a/ChangeLog b/ChangeLog index fe5bd6b..dfaae59 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,14 @@ +1.9.4 - 6/20/2000 + Converted projection files to C++ library with frnetorderstream + Converted image files to C++ library with frnetorderstream + Converted filter and convolution to object-oriented code + Changed default image file to float from double -- Changable in imagefile.h + Optimized image file writing on little-endian architectures + Updated README and INSTALL + 1.9.3 - 6/19/2000 - Reorganized source file + Reorganized source files + MPI bug fixed to phm2pj.cpp 1.9.2 - 6/18/2000 Reorganized include files diff --git a/INSTALL b/INSTALL index 5431237..086d8b3 100644 --- a/INSTALL +++ b/INSTALL @@ -22,6 +22,7 @@ apache (http://www.apache.org) --with-cgibin-dir=..., --with-cgibin-url=..., --with-webdata-dir=..., --with-webdata-url=..., and --with-html-dir=... must be set. + CTSim Specific Configuration Help --------------------------------- @@ -45,11 +46,11 @@ Platforms Tested Recent development is with GNU/Linux. I have tested compilation on FreeBSD v4.0, BSD/OS v3.0, and Solaris v8 (gcc 32-bit), and IA64 (gcc) -, and Microsoft Windows 2000 (Visual C++ 6.0 and cygwin). +, and Microsoft Windows 2000 (Visual C++ 6.0, cygwin, and mingw32). -Microsoft Windows Installation ------------------------------- -Run make.bat from the root directory +Visual C++ Installation +----------------------- +Run make.bat from the root directory (Note, make.bat is out of date) CYGWIN v1.1.0 ------------- @@ -59,6 +60,7 @@ should be deleted. getopt_long appears broken, configure.in checks for cygwin to use bundled version of getopt_long. + Basic Installation ================== diff --git a/README b/README index 2028fad..1bf5d38 100644 --- a/README +++ b/README @@ -1,13 +1,14 @@ COPYRIGHT ---------- +========= This program is written by Kevin M. Rosenberg, M.D. It is covered by the GNU General Public License (GPL) which allows copying and modifying this code with restrictions. See the file COPYING for complete details. + HISTORY ------- +======= CTSim development began in 1983 while I was in medical school. It was written using Lattice C and MS-DOS. I used assembly language to write @@ -16,8 +17,11 @@ directly to an IBM EGA video adapter. In 1999, I ported CTsim to GNU/Linux. In April 2000, the source code for CTSim was published on the Internet. +In June 2000, entire code for revised and converted to C++. + + STATUS ------- +====== The official home for CTsim is http://www.ctsim.org. From this site, you can download the CTsim source code and use CTSim online using a @@ -27,15 +31,13 @@ I would be very pleased to have other developers join me in the development of CTSim. Please see the TODO list for the most obvious ideas. -In this release, the main issue that needs addressing is the inclusion -of old graphic library code. This code is not being used by -CTSim. Compilation warnings will occur in these files. Please ignore -these warnings. In the next major release of CTSim, there will be -support for interactive graphics. I just need to write the low-level -drivers to support X-window rather than IBM EGA hardware. +Interactive graphics have not yet been implemented. I am still researching +cross-platform tools to create a graphical interface with interactive +graphics of the simulation and reconstruction processes. + OVERVIEW --------- +======== CTSim simulates the collection of x-rays by a CT scanner. These x-rays of objects are called projections. @@ -43,13 +45,12 @@ of objects are called projections. Phantom objects are defined. Several built-in phantoms are included, as well as an extension to load files of phantom definitions. -CTsim uses .sdf (standard data files) to represent 2-dimensional image -data. It uses a 32-bit floating-point pixel format. SDF is loosely -based on a image file format used by the Jet Propulsion Laboratory -(JPL). It supports multiple text labels per image. +CTsim uses cross-platform compatible file formats for projection data and +image data. + THE PROGRAMS ------------- +============ phm2if - generates an image file of a phantom object @@ -66,29 +67,30 @@ if2img - Converts an image file to a variety of 8-bit and 16-bit image formats ifinfo - Show statistics and history labels of SDF files + TYPICAL USAGE -------------- +============= When evaluating CT simulation, in general, these steps are followed: Create a phantom image and viewable image file - phm2sdf ... - sdf2img ... + phm2if ... + if2img ... Simulate CT data collection and create a viewable image of raw projections phm2pj ... - pj2sdf ... - sdf2img ... + pj2if ... + if2img ... Perform CT reconstruction and create viewable image file ctrec ... - sdf2img ... + if2img ... -These functions are convert output image to a PNG file performed by -the CGI program ctsim.cgi which can be installed as described in -INSTALL file. +There is a sample shell script installed called 'sample-ctrec.sh' that performs +the above commands.n -There is a sample shell script installed called 'sample-ctrec'. +These functions can be invoked via a web interface with a CGI program +as described in the INSTALL file. CLOSING diff --git a/configure b/configure index 5b11285..f99dc3f 100755 --- a/configure +++ b/configure @@ -711,7 +711,7 @@ fi PACKAGE=ctsim -VERSION=1.9.3 +VERSION=1.9.4 if test "`cd $srcdir && pwd`" != "`pwd`" && test -f $srcdir/config.status; then { echo "configure: error: source directory already configured; run "make distclean" there first" 1>&2; exit 1; } diff --git a/configure.in b/configure.in index b63f09a..b48e0f3 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,1.9.3) +AM_INIT_AUTOMAKE(ctsim,1.9.4) AM_CONFIG_HEADER(config.h) dnl Checks for programs. diff --git a/include/Makefile.am b/include/Makefile.am index 0b7826d..087655f 100644 --- a/include/Makefile.am +++ b/include/Makefile.am @@ -1,4 +1,4 @@ -noinst_HEADERS=ct.h ezplot.h pol.h sgp.h array2d.h imagefile.h backprojectors.h mpiworld.h byteorder.h phantom.h timer.h sstream scanner.h projections.h ctsupport.h +noinst_HEADERS=ct.h ezplot.h pol.h sgp.h array2d.h imagefile.h backprojectors.h mpiworld.h fnetorderstream.h phantom.h timer.h sstream scanner.h projections.h ctsupport.h diff --git a/include/byteorder.h b/include/byteorder.h deleted file mode 100644 index 94fe959..0000000 --- a/include/byteorder.h +++ /dev/null @@ -1,44 +0,0 @@ -#ifndef NETORDER_H -#define NETORDER_H - -#include - -/* netorder.cpp */ - -bool read_nint16 (kuint16 *n, int fd); -bool write_nint16 (kuint16 const *n, int fd); -bool read_nint32 (kuint32 *n, int fd); -bool write_nint32 (kuint32 const *n, int fd); -bool read_nfloat32 (float *f, int fd); -bool write_nfloat32 (float const *f, int fd); -bool read_nfloat64 (double *d, int fd); -bool write_nfloat64 (double const *d, int fd); -void ConvertNetworkOrder (void* buffer, size_t bytes); -void ConvertReverseNetworkOrder (void* buffer, size_t bytes); - -class inetorderstream : public istream { - public: - inetorderstream& readInt16 (kuint16& n); - inetorderstream& readInt32 (kuint32& n); - inetorderstream& readFloat32 (kfloat32& n); - inetorderstream& readFloat64 (kfloat64& n); -}; - -class onetorderstream : public ostream { - public: - onetorderstream& writeInt16 (kuint16 n); - onetorderstream& writeInt32 (kuint32 n); - onetorderstream& writeFloat32 (kfloat32 n); - onetorderstream& writeFloat64 (kfloat64 n); -}; - -void read_rnint16 (kuint16& n, istream istr); -void write_rnint16 (kuint16 n, ostream ostr); -void read_rnint32 (kuint32& n, istream istr); -void write_rnint32 (kuint32 n, ostream ostr); -void read_rnfloat32 (kfloat32& n, istream istr); -void write_rnfloat32 (kfloat32 n, ostream ostr); -void read_rnfloat64 (kfloat64& n, istream istr); -void write_rnfloat64 (kfloat64 n, ostream ostr); - -#endif diff --git a/include/ct.h b/include/ct.h index 016a4f8..4dccb37 100644 --- a/include/ct.h +++ b/include/ct.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ct.h,v 1.19 2000/06/19 19:04:05 kevin Exp $ +** $Id: ct.h,v 1.20 2000/06/20 17:54:51 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 @@ -115,7 +115,7 @@ extern "C" { #endif #include "ctsupport.h" -#include "byteorder.h" +#include "fnetorderstream.h" #ifdef HAVE_SGP #include "ezplot.h" @@ -134,6 +134,7 @@ extern "C" { using namespace std; #include "array2d.h" +#include "fnetorderstream.h" #include "imagefile.h" #include "phantom.h" #include "projections.h" @@ -251,29 +252,19 @@ typedef enum { #include "backprojectors.h" +#include "filter.h" /************************************************************************* * FUNCTION DECLARATIONS ************************************************************************/ -// convolve.cpp -double convolve (const double f1[], const double f2[], const double dx, const int n, const int np, const FunctionSymmetry func_type); - // dialogs.cpp bool phm_add_pelem_kb (Phantom& phm); const Phantom& phm_select (Phantom& phm); int interpolation_select (void); int filter_select (double *filter_param); -// 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 filter_frequency_response (int filt_type, double u, double bw, double param); -double sinc (double x, double mult); -double integral_abscos(double u, double w); - // options.cpp int opt_set_trace(const char *optarg); const char *name_of_phantom(const int phmid); diff --git a/include/ezplot.h b/include/ezplot.h index 40e05f4..11ac4e5 100644 --- a/include/ezplot.h +++ b/include/ezplot.h @@ -2,7 +2,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ezplot.h,v 1.7 2000/06/19 19:04:05 kevin Exp $ +** $Id: ezplot.h,v 1.8 2000/06/20 17:54:51 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 @@ -314,6 +314,7 @@ extern bool ezplot_firstcall; /* set to false on first call to EZSET or EZPLOT /* ezplot.cpp */ +SGP_ID ezplot(float x[], double y[], int num); SGP_ID ezplot(double x[], double y[], int num); void ezinit(void); void ezfree(void); diff --git a/include/fnetorderstream.h b/include/fnetorderstream.h new file mode 100644 index 0000000..d341763 --- /dev/null +++ b/include/fnetorderstream.h @@ -0,0 +1,105 @@ +#ifndef NETORDER_H +#define NETORDER_H + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include +using namespace std; + + +inline bool NativeBigEndian (void) +{ +#ifdef WORDS_BIGENDIAN + return true; +#else + return false; +#endif +} + +inline void +SwapBytes2 (void* buffer) +{ + unsigned char* p = static_cast(buffer); + unsigned char temp = p[0]; + p[0] = p[1]; + p[1] = temp; +} + +// 0<->3 1<->2 = 0123 -> 3210 +inline void +SwapBytes4 (void* buffer) +{ + unsigned char* p = static_cast(buffer); + unsigned char temp = p[0]; + p[0] = p[3]; + p[3] = temp; + temp = p[1]; + p[1] = p[2]; + p[2] = temp; +} + +// 0<->7 1<->6 2<->5 3<->4 = 01234567 -> 76543210 +inline void +SwapBytes8 (void* buffer) +{ + unsigned char* p = static_cast(buffer); + unsigned char temp = p[0]; + p[0] = p[7]; + p[7] = temp; + temp = p[1]; + p[1] = p[6]; + p[6] = temp; + temp = p[2]; + p[2] = p[5]; + p[5] = temp; + temp = p[3]; + p[3] = p[4]; + p[4] = temp; +} + + +void ConvertNetworkOrder (void* buffer, size_t bytes); +void ConvertReverseNetworkOrder (void* buffer, size_t bytes); + + +class fnetorderstream : public fstream { + public: + fnetorderstream (const char* filename, int mode) + : fstream (filename, mode) {} + + ~fnetorderstream (void) + {} + + virtual fnetorderstream& writeInt16 (kuint16 n); + virtual fnetorderstream& writeInt32 (kuint32 n); + virtual fnetorderstream& writeFloat32 (kfloat32 n); + virtual fnetorderstream& writeFloat64 (kfloat64 n); + + virtual fnetorderstream& readInt16 (kuint16& n); + virtual fnetorderstream& readInt32 (kuint32& n); + virtual fnetorderstream& readFloat32 (kfloat32& n); + virtual fnetorderstream& readFloat64 (kfloat64& n); +}; + + +class frnetorderstream : public fnetorderstream { + public: + frnetorderstream (const char* filename, int mode) + : fnetorderstream (filename, mode) {} + + virtual frnetorderstream& writeInt16 (kuint16 n); + virtual frnetorderstream& writeInt32 (kuint32 n); + virtual frnetorderstream& writeFloat32 (kfloat32 n); + virtual frnetorderstream& writeFloat64 (kfloat64 n); + + virtual frnetorderstream& readInt16 (kuint16& n); + virtual frnetorderstream& readInt32 (kuint32& n); + virtual frnetorderstream& readFloat32 (kfloat32& n); + virtual frnetorderstream& readFloat64 (kfloat64& n); +}; + +#endif diff --git a/include/imagefile.h b/include/imagefile.h index 5d36219..607f716 100644 --- a/include/imagefile.h +++ b/include/imagefile.h @@ -4,7 +4,10 @@ #include #include #include +#include +#include #include "ctsupport.h" +#include "fnetorderstream.h" #include "array2d.h" using namespace std; @@ -70,7 +73,7 @@ private: kuint16 headersize; string filename; int file_id; - iostream *io; + frnetorderstream *m_pFS; bool bHeaderWritten; bool bDataWritten; @@ -193,7 +196,6 @@ Array2dFile::init (void) { mPixelSize = sizeof(T); signature = ('I' * 256 + 'F'); - file_id = -1; mNX = 0; mNY = 0; headersize = 0; @@ -204,7 +206,7 @@ Array2dFile::init (void) mMinX = mMaxX = mMinY = mMaxY = 0; mOffsetPV = 0; mScalePV = 1; - io = NULL; + m_pFS = NULL; #if 0 const type_info& t_id = typeid(T); @@ -246,10 +248,10 @@ template void Array2dFile::fileClose (void) { - if (file_id >= 0) { + if (m_pFS) { headerWrite (); - close (file_id); - file_id = -1; + m_pFS->close (); + m_pFS = NULL; } } @@ -257,11 +259,11 @@ template bool Array2dFile::fileCreate (void) { - // io = new iostream(filename, ios::out | ios::in | ios::trunc | io::binary); - if ((file_id = open (filename.c_str(), O_RDWR | O_CREAT | O_TRUNC | O_BINARY, S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH)) < 0) { - sys_error (ERR_WARNING, "Error opening file %s for writing [fileCreate]", filename.c_str()); - return (false); - } + m_pFS = new frnetorderstream (filename.c_str(), ios::out | ios::in | ios::trunc | ios::binary); + if (! m_pFS) { + sys_error (ERR_WARNING, "Error opening file %s for writing [fileCreate]", filename.c_str()); + return (false); + } headerWrite(); return (true); } @@ -270,13 +272,14 @@ template bool Array2dFile::fileRead (void) { - // io = new iostream(filename, ios::out | ios::in | io::binary); - if ((file_id = open (filename.c_str(), O_RDONLY | O_BINARY)) < 0) { + m_pFS = new frnetorderstream (filename.c_str(), ios::out | ios::in | ios::binary | ios::nocreate); + if (m_pFS->fail()) { sys_error (ERR_WARNING, "Unable to open file %s [fileRead]", filename.c_str()); return (false); } - headerRead(); + if (! headerRead()) + return false; array.initSetSize(mNX, mNY); @@ -340,35 +343,35 @@ template bool Array2dFile::headerRead (void) { - if (file_id < 0) { + if (! m_pFS) { sys_error (ERR_WARNING, "Tried to read header with file closed [headerRead]"); return (false); } - lseek (file_id, 0, SEEK_SET); + m_pFS->seekg (0); kuint16 file_signature; kuint16 file_mPixelSize; kuint16 file_mPixelType; - read_nint16 (&headersize, file_id); - read_nint16 (&file_signature, file_id); - read_nint16 (&num_labels, file_id); - read_nint16 (&file_mPixelType, file_id); - read_nint16 (&file_mPixelSize, file_id); - read_nint32 (&mNX, file_id); - read_nint32 (&mNY, file_id); - read_nint16 (&axis_increment_known, file_id); - read_nfloat64 (&mIncX, file_id); - read_nfloat64 (&mIncY, file_id); - read_nint16 (&axis_extent_known, file_id); - read_nfloat64 (&mMinX, file_id); - read_nfloat64 (&mMaxX, file_id); - read_nfloat64 (&mMinY, file_id); - read_nfloat64 (&mMaxY, file_id); - read_nfloat64 (&mOffsetPV, file_id); - read_nfloat64 (&mScalePV, file_id); - - int read_headersize = lseek (file_id, 0, SEEK_CUR); + m_pFS->readInt16 (headersize); + m_pFS->readInt16 (file_signature); + m_pFS->readInt16 (num_labels); + m_pFS->readInt16 (file_mPixelType); + m_pFS->readInt16 (file_mPixelSize); + m_pFS->readInt32 (mNX); + m_pFS->readInt32 (mNY); + m_pFS->readInt16 (axis_increment_known); + m_pFS->readFloat64 (mIncX); + m_pFS->readFloat64 (mIncY); + m_pFS->readInt16 (axis_extent_known); + m_pFS->readFloat64 (mMinX); + m_pFS->readFloat64 (mMaxX); + m_pFS->readFloat64 (mMinY); + m_pFS->readFloat64 (mMaxY); + m_pFS->readFloat64 (mOffsetPV); + m_pFS->readFloat64 (mScalePV); + + int read_headersize = m_pFS->tellg(); if (read_headersize != headersize) { sys_error (ERR_WARNING, "Read headersize %d != file headersize %d", read_headersize, headersize); return (false); @@ -393,33 +396,32 @@ template bool Array2dFile::headerWrite (void) { - if (file_id < 0) { + if (! m_pFS) { sys_error (ERR_WARNING, "Tried to write header with file closed"); return (false); } - lseek (file_id, 0, SEEK_SET); - write_nint16 (&headersize, file_id); - write_nint16 (&signature, file_id); - write_nint16 (&num_labels, file_id); - write_nint16 (&mPixelType, file_id); - write_nint16 (&mPixelSize, file_id); - write_nint32 (&mNX, file_id); - write_nint32 (&mNY, file_id); - write_nint16 (&axis_increment_known, file_id); - write_nfloat64 (&mIncX, file_id); - write_nfloat64 (&mIncY, file_id); - write_nint16 (&axis_extent_known, file_id); - write_nfloat64 (&mMinX, file_id); - write_nfloat64 (&mMaxX, file_id); - write_nfloat64 (&mMinY, file_id); - write_nfloat64 (&mMaxY, file_id); - write_nfloat64 (&mOffsetPV, file_id); - write_nfloat64 (&mScalePV, file_id); - - headersize = lseek (file_id, 0, SEEK_CUR); - lseek (file_id, 0, SEEK_SET); - write_nint16 (&headersize, file_id); + m_pFS->seekp (0); + m_pFS->writeInt16 (headersize); + m_pFS->writeInt16 (signature); + m_pFS->writeInt16 (num_labels); + m_pFS->writeInt16 (mPixelType); + m_pFS->writeInt16 (mPixelSize); + m_pFS->writeInt32 (mNX); + m_pFS->writeInt32 (mNY); + m_pFS->writeInt16 (axis_increment_known); + m_pFS->writeFloat64 (mIncX); + m_pFS->writeFloat64 (mIncY); + m_pFS->writeInt16 (axis_extent_known); + m_pFS->writeFloat64 (mMinX); + m_pFS->writeFloat64 (mMaxX); + m_pFS->writeFloat64 (mMinY); + m_pFS->writeFloat64 (mMaxY); + m_pFS->writeFloat64 (mOffsetPV); + m_pFS->writeFloat64 (mScalePV); + + headersize = m_pFS->tellp(); + m_pFS->writeInt16 (headersize); return (true); } @@ -428,8 +430,8 @@ template bool Array2dFile::arrayDataWrite (void) { - if (file_id < 0) { - sys_error (ERR_WARNING, "Tried to arrayDataWrite with file_id < 0"); + if (! m_pFS) { + sys_error (ERR_WARNING, "Tried to arrayDataWrite with !m_pFS"); return (false); } @@ -437,13 +439,17 @@ Array2dFile::arrayDataWrite (void) if (! da) return (false); - lseek (file_id, headersize, SEEK_SET); - for (unsigned int ix = 0; ix < mNX; ix++) - for (unsigned int iy = 0; iy < mNY; iy++) { - T value = da[ix][iy]; - ConvertReverseNetworkOrder (&value, sizeof(T)); - write (file_id, &value, mPixelSize); - } + m_pFS->seekp (headersize); + for (unsigned int ix = 0; ix < mNX; ix++) { + if (NativeBigEndian()) { + for (unsigned int iy = 0; iy < mNY; iy++) { + T value = da[ix][iy]; + ConvertReverseNetworkOrder (&value, sizeof(T)); + m_pFS->write (&value, mPixelSize); + } + } else + m_pFS->write (da[ix], sizeof(T) * mNY); + } return (true); } @@ -452,8 +458,8 @@ template bool Array2dFile::arrayDataRead (void) { - if (file_id < 0) { - sys_error (ERR_WARNING, "Tried to arrayDataRead with file_id < 0"); + if (! m_pFS) { + sys_error (ERR_WARNING, "Tried to arrayDataRead with !m_pFS"); return (false); } @@ -461,14 +467,19 @@ Array2dFile::arrayDataRead (void) if (! da) return (false); - lseek (file_id, headersize, SEEK_SET); + m_pFS->seekg (headersize); T pixelBuffer; - for (int ix = 0; ix < mNX; ix++) - for (unsigned int iy = 0; iy < mNY; iy++) { - read (file_id, &pixelBuffer, mPixelSize); - ConvertReverseNetworkOrder (&pixelBuffer, sizeof(T)); - da[ix][iy] = pixelBuffer; - } + for (int ix = 0; ix < mNX; ix++) { + if (NativeBigEndian()) { + for (unsigned int iy = 0; iy < mNY; iy++) { + m_pFS->read (&pixelBuffer, mPixelSize); + ConvertReverseNetworkOrder (&pixelBuffer, sizeof(T)); + da[ix][iy] = pixelBuffer; + } + } else + m_pFS->read (da[ix], sizeof(T) * mNY); + } + return (true); } @@ -488,26 +499,21 @@ Array2dFile::labelSeek (unsigned int label_num) } off_t pos = headersize + array.sizeofArray(); - if (lseek (file_id, pos, SEEK_SET) != pos) { - sys_error (ERR_WARNING, "Can't seek to end of data array"); - return (false); - } + m_pFS->seekg (pos); for (int i = 0; i < label_num; i++) { pos += 22; // Skip to string length - if (lseek (file_id, pos, SEEK_SET) != pos) { - sys_error (ERR_WARNING, "Can't seek to string length"); - return (false); - } + m_pFS->seekg (pos); + kuint16 strlength; - read_nint16 (&strlength, file_id); + m_pFS->readInt16 (strlength); pos += 2 + strlength; // Skip label string length + data - if (lseek (file_id, pos, SEEK_SET) != pos) { - sys_error (ERR_WARNING, "Can't seek past label string"); - return (false); - } + m_pFS->seekg (pos); } + + if (! m_pFS) + return false; return (true); } @@ -526,19 +532,19 @@ Array2dFile::labelRead (Array2dFileLabel& label, unsigned int label_num) return (false); } - read_nint16 (&label.m_labelType, file_id); - read_nint16 (&label.m_year, file_id); - read_nint16 (&label.m_month, file_id); - read_nint16 (&label.m_day, file_id); - read_nint16 (&label.m_hour, file_id); - read_nint16 (&label.m_minute, file_id); - read_nint16 (&label.m_second, file_id); - read_nfloat64 (&label.m_calcTime, file_id); + m_pFS->readInt16 (label.m_labelType); + m_pFS->readInt16 (label.m_year); + m_pFS->readInt16 (label.m_month); + m_pFS->readInt16 (label.m_day); + m_pFS->readInt16 (label.m_hour); + m_pFS->readInt16 (label.m_minute); + m_pFS->readInt16 (label.m_second); + m_pFS->readFloat64 (label.m_calcTime); kuint16 strlength; - read_nint16 (&strlength, file_id); + m_pFS->readInt16 (strlength); char str [strlength+1]; - read (file_id, str, strlength); + m_pFS->read (str, strlength); str[strlength] = 0; label.m_strLabel = str; @@ -568,22 +574,21 @@ Array2dFile::labelAdd (const Array2dFileLabel& label) { labelSeek (num_labels); - write_nint16 (&label.m_labelType, file_id); - write_nint16 (&label.m_year, file_id); - write_nint16 (&label.m_month, file_id); - write_nint16 (&label.m_day, file_id); - write_nint16 (&label.m_hour, file_id); - write_nint16 (&label.m_minute, file_id); - write_nint16 (&label.m_second, file_id); - write_nfloat64 (&label.m_calcTime, file_id); + m_pFS->writeInt16 (label.m_labelType); + m_pFS->writeInt16 (label.m_year); + m_pFS->writeInt16 (label.m_month); + m_pFS->writeInt16 (label.m_day); + m_pFS->writeInt16 (label.m_hour); + m_pFS->writeInt16 (label.m_minute); + m_pFS->writeInt16 (label.m_second); + m_pFS->writeFloat64 (label.m_calcTime); kuint16 strlength = label.m_strLabel.length(); - write_nint16 (&strlength, file_id); - write (file_id, static_cast(label.m_strLabel.c_str()), strlength); + m_pFS->writeInt16 (strlength); + m_pFS->write (label.m_strLabel.c_str(), strlength); num_labels++; headerWrite(); - fsync(file_id); } template @@ -674,7 +679,7 @@ class F64Image : public Array2dFile #endif }; -#define IMAGEFILE_64_BITS 1 +#undef IMAGEFILE_64_BITS #ifdef IMAGEFILE_64_BITS typedef F64Image ImageFile; typedef kfloat64 ImageFileValue; diff --git a/include/projections.h b/include/projections.h index c1302a9..104476d 100644 --- a/include/projections.h +++ b/include/projections.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: projections.h,v 1.2 2000/06/19 17:58:20 kevin Exp $ +** $Id: projections.h,v 1.3 2000/06/20 17:54:51 kevin Exp $ ** ** ** This program is free software; you can redistribute it and/or modify @@ -48,8 +48,8 @@ class Projections bool read (const char* fname); bool write (const char* fname); - bool detarrayRead (DetectorArray& darray, const int view_num); - bool detarrayWrite (const DetectorArray& darray, const int view_num); + bool detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int view_num); + bool detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int view_num); void setNView (int nView); // used in MPI to restrict # of views void setRotInc (double rotInc) { m_rotInc = rotInc;} @@ -72,7 +72,6 @@ class Projections { return (*m_projData[iview]); } private: - int m_fd; // Disk file descriptor int m_headerSize; // Size of disk file header int m_geometry; // Geometry of scanner struct DetectorArray **m_projData; // Pointer to array of detarray_st pointers @@ -88,6 +87,8 @@ class Projections bool headerRead (void); bool headerWrite (void); + bool headerRead (fnetorderstream& fs); + bool headerWrite (fnetorderstream& fs); void newProjData (void); void deleteProjData (void); diff --git a/libctgraphics/ezplot.cpp b/libctgraphics/ezplot.cpp index dfd0baf..81af0dd 100644 --- a/libctgraphics/ezplot.cpp +++ b/libctgraphics/ezplot.cpp @@ -6,7 +6,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ezplot.cpp,v 1.2 2000/06/19 19:04:05 kevin Exp $ +** $Id: ezplot.cpp,v 1.3 2000/06/20 17:54:51 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 @@ -123,6 +123,18 @@ static char y_numfmt[10]; /* format to print y tick labels */ *----------------------------------------------------------------------*/ +SGP_ID +ezplot (float x[], double y[], int num) +{ + double dx [num]; + + for (int i = 0; i < num; i++) + dx[i] = x[i]; + + ezplot (dx, y, num); +} + + SGP_ID ezplot (double x[], double y[], int num) { diff --git a/libctgraphics/ezset.cpp b/libctgraphics/ezset.cpp index 1bc08c8..09e85e1 100644 --- a/libctgraphics/ezset.cpp +++ b/libctgraphics/ezset.cpp @@ -6,7 +6,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ezset.cpp,v 1.2 2000/06/19 19:04:05 kevin Exp $ +** $Id: ezset.cpp,v 1.3 2000/06/20 17:54:51 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 @@ -183,7 +183,7 @@ do_cmd (int lx) break; case S_REPLOT: ez.i_plotimmediate = TRUE; - ezplot (NULL, NULL, 0); + ezplot (static_cast(NULL), static_cast(NULL), 0); #if 0 if (modeinteract == TRUE) WAITKEY(); @@ -405,7 +405,7 @@ do_cmd (int lx) ez.o_unknowncurves = FALSE; ez.o_reqcurves = ez.i_numcurves; ez.i_plotimmediate = TRUE; - ezplot (NULL, NULL, 0); + ezplot (static_cast(NULL), static_cast(NULL), 0); ez.i_plotimmediate = FALSE; } } diff --git a/libctsim/Makefile.am b/libctsim/Makefile.am index 3cc76ac..70b8fa5 100644 --- a/libctsim/Makefile.am +++ b/libctsim/Makefile.am @@ -1,5 +1,5 @@ noinst_LIBRARIES = libctsim.a -libctsim_a_SOURCES = filter.cpp scanner.cpp projections.cpp phantom.cpp options.cpp convolve.cpp reconstr.cpp imagefile.cpp reconstr.cpp backprojectors.cpp +libctsim_a_SOURCES = filter.cpp scanner.cpp projections.cpp phantom.cpp options.cpp imagefile.cpp reconstr.cpp backprojectors.cpp INCLUDES=@my_includes@ EXTRA_DIST=Makefile.nt diff --git a/libctsim/convolve.cpp b/libctsim/convolve.cpp deleted file mode 100644 index 52abd94..0000000 --- a/libctsim/convolve.cpp +++ /dev/null @@ -1,84 +0,0 @@ -/***************************************************************************** -** FILE IDENTIFICATION -** -** Name: Convolve.c Convolution functions -** Date Started: 13 Mar 86 -** -** This is part of the CTSim program -** Copyright (C) 1983-2000 Kevin Rosenberg -** -** This program is free software; you can redistribute it and/or modify -** it under the terms of the GNU General Public License (version 2) as -** published by the Free Software Foundation. -** -** This program is distributed in the hope that it will be useful, -** but WITHOUT ANY WARRANTY; without even the implied warranty of -** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -** GNU General Public License for more details. -** -** You should have received a copy of the GNU General Public License -** along with this program; if not, write to the Free Software -** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -******************************************************************************/ - -#include "ct.h" - - -/* NAME - * convolve Discrete convolution of two functions - * - * SYNOPSIS - * r = convolve (f1, f2, dx, n, np, func_type) - * double r Convolved result - * double f1[], f2[] Functions to be convolved - * double dx Difference between successive x values - * int n Array index to center convolution about - * int np Number of points in f1 array - * int func_type EVEN or ODD or EVEN_AND_ODD function f2 - * - * NOTES - * f1 is the projection data, its indices range from 0 to np - 1. - * The index for f2, the filter, ranges from -(np-1) to (np-1). - * There are 3 ways to handle the negative vertices of f2: - * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n]. - * All filters used in reconstruction are even. - * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n] - * 3. If f2 is both ODD AND EVEN, then we must store the value of f2 - * for negative indices. Since f2 must range from -(np-1) to (np-1), - * if we add (np - 1) to f2's array index, then f2's index will - * range from 0 to 2 * (np - 1), and the origin, x = 0, will be - * stored at f2[np-1]. - */ - -double -convolve (const double f1[], const double f2[], const double dx, const int n, const int np, const FunctionSymmetry func_type) -{ - double sum = 0.0; - - if (func_type == FUNC_BOTH) { -#if UNOPTIMIZED_CONVOLUTION - for (int i = 0; i < np; i++) - sum += f1[i] * f2[n - i + (np - 1)]; -#else - f2 += n + (np - 1); - for (int i = 0; i < np; i++) - sum += *f1++ * *f2--; -#endif - } else if (func_type == FUNC_EVEN) { - for (int i = 0; i < np; i++) { - int k = abs (n - i); - sum += f1[i] * f2[k]; - } - } else if (func_type == FUNC_ODD) { - for (int i = 0; i < np; i++) { - int k = n - i; - if (k < 0) - sum -= f1[i] * f2[k]; - else - sum += f1[i] * f2[k]; - } - } else - sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type); - - return (sum * dx); -} diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp index 30d4718..b2d3062 100644 --- a/libctsim/filter.cpp +++ b/libctsim/filter.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: filter.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $ +** $Id: filter.cpp,v 1.2 2000/06/20 17:54:51 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 @@ -46,40 +46,44 @@ * ANALYTIC solutions, use numint = 0 */ -double * -filter_generate (const FilterType filt_type, double bw, double xmin, double xmax, int n, double param, const DomainType domain, int numint) +SignalFilter::SignalFilter (const FilterType filt_type, double bw, double xmin, double xmax, int n, double param, const DomainType domain, int numint) { - double *f = new double [n]; + m_vecFilter = new double [n]; + m_filterType = filt_type; + m_bw = bw; + double xinc = (xmax - xmin) / (n - 1); - if (filt_type == FILTER_SHEPP) { - double a = 2 * bw; + if (m_filterType == FILTER_SHEPP) { + double a = 2 * m_bw; double c = - 4. / (a * a); int center = (n - 1) / 2; int sidelen = center; - f[center] = 4. / (a * a); + m_vecFilter[center] = 4. / (a * a); for (int i = 1; i <= sidelen; i++ ) - f [center + i] = f [center - i] = c / (4 * (i * i) - 1); + m_vecFilter [center + i] = m_vecFilter [center - i] = c / (4 * (i * i) - 1); } else if (domain == D_FREQ) { double x; int i; for (x = xmin, i = 0; i < n; x += xinc, i++) - f[i] = filter_frequency_response (filt_type, x, bw, param); + m_vecFilter[i] = frequencyResponse (x, param); } else if (domain == D_SPATIAL) { double x; int i; for (x = xmin, i = 0; i < n; x += xinc, i++) if (numint == 0) - f[i] = filter_spatial_response_analytic (filt_type, x, bw, param); + m_vecFilter[i] = spatialResponseAnalytic (x, param); else - f[i] = filter_spatial_response_calc (filt_type, x, bw, param, numint); + m_vecFilter[i] = spatialResponseCalc (x, param, numint); } else { sys_error (ERR_WARNING, "Illegal domain %d [filt_generate]", domain); - return (NULL); } - - return (f); +} + +SignalFilter::~SignalFilter (void) +{ + delete m_vecFilter; } @@ -89,21 +93,27 @@ filter_generate (const FilterType filt_type, double bw, double xmin, double xmax * response * * SYNOPSIS - * y = filter_spatial_response_calc (filt_type, x, bw, param, n) + * y = filter_spatial_response_calc (filt_type, x, m_bw, param, n) * double y Filter's response in spatial domain * int filt_type Type of filter (definitions in ct.h) * double x Spatial position to evaluate filter - * double bw Bandwidth of window + * double m_bw Bandwidth of window * double param General parameter for various filters * int n Number of points to calculate integrations */ double -filter_spatial_response_calc (int filt_type, double x, double bw, double param, int n) +SignalFilter::spatialResponseCalc (double x, double param, int n) const +{ + return (spatialResponseCalc (m_filterType, m_bw, x, param, n)); +} + +double +SignalFilter::spatialResponseCalc (FilterType fType, double bw, double x, double param, int n) { double zmin, zmax; - if (filt_type == FILTER_TRIANGLE) { + if (fType == FILTER_TRIANGLE) { zmin = 0; zmax = bw; } else { @@ -115,7 +125,7 @@ filter_spatial_response_calc (int filt_type, double x, double bw, double param, double z = zmin; double q [n]; for (int i = 0; i < n; i++, z += zinc) - q[i] = filter_frequency_response (filt_type, z, bw, param) * cos (TWOPI * z * x); + q[i] = frequencyResponse (fType, bw, z, param) * cos (TWOPI * z * x); double y = 2 * integrateSimpson (zmin, zmax, q, n); @@ -127,21 +137,28 @@ filter_spatial_response_calc (int filt_type, double x, double bw, double param, * filter_frequency_response Return filter frequency response * * SYNOPSIS - * h = filter_frequency_response (filt_type, u, bw, param) + * h = filter_frequency_response (filt_type, u, m_bw, param) * double h Filters frequency response at u * int filt_type Type of filter * double u Frequency to evaluate filter at - * double bw Bandwidth of filter + * double m_bw Bandwidth of filter * double param General input parameter for various filters */ double -filter_frequency_response (int filt_type, double u, double bw, double param) +SignalFilter::frequencyResponse (double u, double param) const +{ + return frequencyResponse (m_filterType, m_bw, u, param); +} + + +double +SignalFilter::frequencyResponse (FilterType fType, double bw, double u, double param) { double q; double au = fabs (u); - switch (filt_type) { + switch (fType) { case FILTER_BANDLIMIT: if (au >= bw / 2) q = 0.; @@ -192,9 +209,7 @@ filter_frequency_response (int filt_type, double u, double bw, double param) break; default: q = 0; - sys_error (ERR_WARNING, - "Frequency response for filter %d not implemented [filter_frequency_response]", - filt_type); + sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", fType); break; } return (q); @@ -208,16 +223,22 @@ filter_frequency_response (int filt_type, double u, double bw, double param) * response * * SYNOPSIS - * y = filter_spatial_response_analytic (filt_type, x, bw, param) + * y = filter_spatial_response_analytic (filt_type, x, m_bw, param) * double y Filter's response in spatial domain * int filt_type Type of filter (definitions in ct.h) * double x Spatial position to evaluate filter - * double bw Bandwidth of window + * double m_bw Bandwidth of window * double param General parameter for various filters */ double -filter_spatial_response_analytic (int filt_type, double x, double bw, double param) +SignalFilter::spatialResponseAnalytic (double x, double param) const +{ + return spatialResponseAnalytic (m_filterType, m_bw, x, param); +} + +double +SignalFilter::spatialResponseAnalytic (FilterType fType, double bw, double x, double param) { double q, temp; double u = TWOPI * x; @@ -225,7 +246,7 @@ filter_spatial_response_analytic (int filt_type, double x, double bw, double par double b = PI / bw; double b2 = TWOPI / bw; - switch (filt_type) { + switch (fType) { case FILTER_BANDLIMIT: q = bw * sinc(u * w, 1.0); break; @@ -237,8 +258,7 @@ filter_spatial_response_analytic (int filt_type, double x, double bw, double par q = sinc(b-u,w) + sinc(b+u,w); break; case FILTER_G_HAMMING: - q = 2 * param * sin(u*w)/u + (1-param) * - (sinc(b2-u, w) + sinc(b2+u, w)); + q = 2 * param * sin(u*w)/u + (1-param) * (sinc(b2-u, w) + sinc(b2+u, w)); break; case FILTER_ABS_BANDLIMIT: q = 2 * integral_abscos (u, w); @@ -264,9 +284,7 @@ filter_spatial_response_analytic (int filt_type, double x, double bw, double par break; case FILTER_ABS_SINC: default: - sys_error (ERR_WARNING, - "Analytic filter type %d not implemented [filter_spatial_response_analytic]", - filt_type); + sys_error (ERR_WARNING, "Analytic filter type %d not implemented [filter_spatial_response_analytic]", fType); q = 0; break; } @@ -287,12 +305,6 @@ filter_spatial_response_analytic (int filt_type, double x, double bw, double par * v = sin(x * mult) / x; */ -double -sinc (double x, double mult) -{ - return (fabs(x) > F_EPSILON ? (sin (x * mult) / x) : 1.0); -} - /* NAME * integral_abscos Returns integral of u*cos(u) @@ -308,7 +320,7 @@ sinc (double x, double mult) */ double -integral_abscos (double u, double w) +SignalFilter::integral_abscos (double u, double w) { if (fabs (u) > F_EPSILON) return (cos(u * w) - 1) / (u * u) + w / u * sin (u * w); @@ -317,3 +329,96 @@ integral_abscos (double u, double w) } +/* NAME + * convolve Discrete convolution of two functions + * + * SYNOPSIS + * r = convolve (f1, f2, dx, n, np, func_type) + * double r Convolved result + * double f1[], f2[] Functions to be convolved + * double dx Difference between successive x values + * int n Array index to center convolution about + * int np Number of points in f1 array + * int func_type EVEN or ODD or EVEN_AND_ODD function f2 + * + * NOTES + * f1 is the projection data, its indices range from 0 to np - 1. + * The index for f2, the filter, ranges from -(np-1) to (np-1). + * There are 3 ways to handle the negative vertices of f2: + * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n]. + * All filters used in reconstruction are even. + * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n] + * 3. If f2 is both ODD AND EVEN, then we must store the value of f2 + * for negative indices. Since f2 must range from -(np-1) to (np-1), + * if we add (np - 1) to f2's array index, then f2's index will + * range from 0 to 2 * (np - 1), and the origin, x = 0, will be + * stored at f2[np-1]. + */ + +double +SignalFilter::convolve (const double func[], const double dx, const int n, const int np, const FunctionSymmetry func_type) const +{ + double sum = 0.0; + + if (func_type == FUNC_BOTH) { +#if UNOPTIMIZED_CONVOLUTION + for (int i = 0; i < np; i++) + sum += func[i] * m_vecFilter[n - i + (np - 1)]; +#else + double* f2 = m_vecFilter + n + (np - 1); + for (int i = 0; i < np; i++) + sum += *func++ * *f2--; +#endif + } else if (func_type == FUNC_EVEN) { + for (int i = 0; i < np; i++) { + int k = abs (n - i); + sum += func[i] * m_vecFilter[k]; + } + } else if (func_type == FUNC_ODD) { + for (int i = 0; i < np; i++) { + int k = n - i; + if (k < 0) + sum -= func[i] * m_vecFilter[k]; + else + sum += func[i] * m_vecFilter[k]; + } + } else + sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type); + + return (sum * dx); +} + + +double +SignalFilter::convolve (const float func[], const double dx, const int n, const int np, const FunctionSymmetry func_type) const +{ + double sum = 0.0; + + if (func_type == FUNC_BOTH) { +#if UNOPTIMIZED_CONVOLUTION + for (int i = 0; i < np; i++) + sum += func[i] * m_vecFilter[n - i + (np - 1)]; +#else + double* f2 = m_vecFilter + n + (np - 1); + for (int i = 0; i < np; i++) + sum += *func++ * *f2--; +#endif + } else if (func_type == FUNC_EVEN) { + for (int i = 0; i < np; i++) { + int k = abs (n - i); + sum += func[i] * m_vecFilter[k]; + } + } else if (func_type == FUNC_ODD) { + for (int i = 0; i < np; i++) { + int k = n - i; + if (k < 0) + sum -= func[i] * m_vecFilter[k]; + else + sum += func[i] * m_vecFilter[k]; + } + } else + sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type); + + return (sum * dx); +} + diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp index 8756bc8..4cae5d2 100644 --- a/libctsim/imagefile.cpp +++ b/libctsim/imagefile.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: imagefile.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $ +** $Id: imagefile.cpp,v 1.2 2000/06/20 17:54:51 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 @@ -110,9 +110,9 @@ image_filter_response (ImageFile& im, const DomainType domain, double bw, const double r = sqrt(i * i + j * j); if (domain == D_SPATIAL) - v[i+hx][j+hy] = filter_spatial_response_analytic (filt_type, r, bw, filt_param); + v[i+hx][j+hy] = SignalFilter::spatialResponseAnalytic (filt_type, bw, r, filt_param); else - v[i+hx][j+hy] = filter_frequency_response (filt_type, r, bw, filt_param); + v[i+hx][j+hy] = SignalFilter::frequencyResponse (filt_type, bw, r, filt_param); if (opt_trace >= TRACE_PHM) printf ("r=%8.4f, v=%8.4f\n", r, v[i+hx][j+hy]); } diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 8010912..54f8f7d 100644 --- a/libctsim/projections.cpp +++ b/libctsim/projections.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: projections.cpp,v 1.2 2000/06/19 17:58:20 kevin Exp $ +** $Id: projections.cpp,v 1.3 2000/06/20 17:54:51 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 @@ -137,7 +137,7 @@ Projections::deleteProjData (void) */ bool -Projections::headerWrite (void) +Projections::headerWrite (fnetorderstream& fs) { kuint32 _hsize = m_headerSize; kuint32 _nView = m_nView; @@ -151,48 +151,39 @@ Projections::headerWrite (void) kfloat64 _detInc = m_detInc; kfloat64 _phmLen = m_phmLen; - if (m_fd < 0) { - sys_error (ERR_SEVERE, "m_fd < 0 [projections_write_header]"); + fs.seekp(0); + if (! fs) return false; - } - - if (lseek (m_fd, (off_t) 0, SEEK_SET) != (off_t) 0) { - sys_error (ERR_WARNING, "Error seeking to beginning of fiel [projections_write_header]"); - return false; - } - if (! write_nint32 (&_hsize, m_fd) || - ! write_nint32 (&_nView, m_fd) || - ! write_nint32 (&_nDet, m_fd) || - ! write_nint32 (&_geom, m_fd) || - ! write_nfloat64 (&_calcTime, m_fd) || - ! write_nfloat64 (&_rotStart, m_fd) || - ! write_nfloat64 (&_rotInc, m_fd) || - ! write_nfloat64 (&_detStart, m_fd) || - ! write_nfloat64 (&_detInc, m_fd) || - ! write_nfloat64 (&_phmLen, m_fd) || - ! write_nint32 (&_remarksize, m_fd) || - (::write (m_fd, m_remark.c_str(), _remarksize) != _remarksize)) { - sys_error (ERR_SEVERE, "Error writing header information [projections_write_header] %ld", _remarksize); - return false; - } - m_headerSize = _hsize = lseek(m_fd, (off_t) 0, SEEK_CUR); - if ((lseek(m_fd, 0, SEEK_SET) != 0) || ! write_nint32 (&_hsize, m_fd)) { - sys_error (ERR_SEVERE, "Error writing header information [projections_write_header]"); - return false; - } + fs.writeInt32 (_hsize); + fs.writeInt32 (_nView); + fs.writeInt32 (_nDet); + fs.writeInt32 (_geom); + fs.writeFloat64 (_calcTime); + fs.writeFloat64 (_rotStart); + fs.writeFloat64 (_rotInc); + fs.writeFloat64 (_detStart); + fs.writeFloat64 (_detInc); + fs.writeFloat64 (_phmLen); + fs.writeInt32 (_remarksize); + fs.write (m_remark.c_str(), _remarksize); + + m_headerSize = fs.tellp(); + _hsize = m_headerSize; + fs.seekp(0); + fs.writeInt32 (_hsize); + if (! fs) + return false; return true; } - /* NAME * projections_read_header Read data header for projections file * */ - bool -Projections::headerRead (void) +Projections::headerRead (fnetorderstream& fs) { kuint32 _hsize; kuint32 _nView; @@ -206,46 +197,43 @@ Projections::headerRead (void) kfloat64 _detInc; kfloat64 _phmLen; - if (m_fd < 0) { - sys_error (ERR_SEVERE, "m_fd < 0 [projections_read_header]"); - return false; - } - - if (lseek (m_fd, (off_t) 0, SEEK_SET) != (off_t) 0) { - sys_error (ERR_WARNING, "Error seeking to beginning of field [projections_read_header]"); - return false; - } + fs.seekg(0); + if (! fs) + return false; off_t testPos; - if (! read_nint32 (&_hsize, m_fd)) { - sys_error (ERR_SEVERE, "Error reading header size , _hsize=%d [projections_read_header]", _hsize); + fs.readInt32 (_hsize); + fs.readInt32 (_nView); + fs.readInt32 (_nDet); + fs.readInt32 (_geom); + fs.readFloat64 (_calcTime); + fs.readFloat64 (_rotStart); + fs.readFloat64 (_rotInc); + fs.readFloat64 (_detStart); + fs.readFloat64 (_detInc); + fs.readFloat64 (_phmLen); + fs.readInt32 (_remarksize); + if (_remarksize > 100000) { + sys_error (ERR_SEVERE, "Insane _remarksize %d [projections::headerRead]", _remarksize); return false; } - if ( ! read_nint32 (&_nView, m_fd) || - ! read_nint32 (&_nDet, m_fd) || - ! read_nint32 (&_geom, m_fd) || - ! read_nfloat64 (&_calcTime, m_fd) || - ! read_nfloat64 (&_rotStart, m_fd) || - ! read_nfloat64 (&_rotInc, m_fd) || - ! read_nfloat64 (&_detStart, m_fd) || - ! read_nfloat64 (&_detInc, m_fd) || - ! read_nfloat64 (&_phmLen, m_fd) || - ! read_nint32 (&_remarksize, m_fd)) - { - sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize); - return false; - } + + if (! fs) { + sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize); + return false; + } char remarkStorage[_remarksize+1]; - if (::read (m_fd, &remarkStorage, _remarksize) != _remarksize) { + fs.read (remarkStorage, _remarksize); + if (! fs) { sys_error (ERR_SEVERE, "Error reading remark, _remarksize = %d", _remarksize); return false; } remarkStorage[_remarksize] = 0; m_remark = remarkStorage; - off_t _hsizeread = lseek(m_fd, 0, SEEK_CUR); - if (_hsizeread != _hsize) { + off_t _hsizeread = fs.tellg(); + if (!fs || _hsizeread != _hsize) { sys_error (ERR_WARNING, "File header size read %ld != file header size stored %ld [read_projections_header]\n_remarksize=%ld", (long int) _hsizeread, _hsize, _remarksize); return false; } @@ -264,31 +252,26 @@ Projections::headerRead (void) return true; } - bool Projections::read (const char* filename) { - if ((m_fd = open(filename, OPEN_RDONLY | O_BINARY)) < 0) { - sys_error (ERR_SEVERE, "Error opening projections file %s for input [open_projections]", filename); + frnetorderstream fileRead (filename, ios::in | ios::binary); + + if (! fileRead) return false; - } - if (! headerRead()) { - sys_error (ERR_SEVERE, "Error reading projections file %s [open_projections]", filename); + if (! headerRead (fileRead)) return false; - } deleteProjData (); - newProjData (); + newProjData(); for (int i = 0; i < m_nView; i++) { - if (! detarrayRead (*m_projData[i], i)) + if (! detarrayRead (fileRead, *m_projData[i], i)) break; } - close (m_fd); - m_fd = -1; - + fileRead.close(); return true; } @@ -296,30 +279,25 @@ Projections::read (const char* filename) bool Projections::write (const char* filename) { - if ((m_fd = open (filename, O_WRONLY | O_CREAT | O_TRUNC | O_BINARY -#ifndef _WIN32 - , S_IWRITE | S_IREAD | S_IRGRP | S_IROTH -#endif - )) < 0) { - sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename); - return false; - } - if (! headerWrite ()) { - sys_error(ERR_SEVERE, "Error writing to file %s [create_projections]", filename); - return false; - } + frnetorderstream fs (filename, ios::out | ios::binary | ios::trunc | ios::ate); - headerWrite (); + if (! fs) { + sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename); + return false; + } + if (! headerWrite (fs)) + return false; if (m_projData != NULL) { for (int i = 0; i < m_nView; i++) { - if (! detarrayWrite (*m_projData[i], i)) + if (! detarrayWrite (fs, *m_projData[i], i)) break; } } + if (! fs) + return false; - close (m_fd); - m_fd = -1; + fs.close(); return true; } @@ -334,7 +312,7 @@ Projections::write (const char* filename) */ bool -Projections::detarrayRead (DetectorArray& darray, const int iview) +Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int iview) { const int detval_bytes = darray.nDet() * sizeof(kfloat32); const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */; @@ -344,32 +322,20 @@ Projections::detarrayRead (DetectorArray& darray, const int iview) kfloat64 view_angle; kuint32 nDet; - if (m_fd < 0) { - sys_error (ERR_SEVERE, "Tried to read detector array on closed file [detarrayRead]"); - return false; - } + fs.seekg (start_data); - if (lseek(m_fd, start_data, SEEK_SET) != start_data) { - sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayRead]"); - return false; - } - - if (! read_nfloat64 (&view_angle, m_fd) || - ! read_nint32 (&nDet, m_fd)) { - sys_error (ERR_SEVERE, "Error reading detector array [detarrayRead]"); - return false; - } + fs.readFloat64 (view_angle); + fs.readInt32 (nDet); darray.setViewAngle (view_angle); // darray.setNDet ( nDet); for (int i = 0; i < nDet; i++) { kfloat32 detval; - if (! read_nfloat32 (&detval, m_fd)) { - sys_error (ERR_SEVERE, "Error reading detector array [detarrayRead]"); - return false; - } + fs.readFloat32 (detval); detval_ptr[i] = detval; } + if (! fs) + return false; return true; } @@ -390,7 +356,7 @@ Projections::detarrayRead (DetectorArray& darray, const int iview) */ bool -Projections::detarrayWrite (const DetectorArray& darray, const int iview) +Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int iview) { const int detval_bytes = darray.nDet() * sizeof(float); const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */; @@ -400,36 +366,26 @@ Projections::detarrayWrite (const DetectorArray& darray, const int iview) kfloat64 view_angle = darray.viewAngle(); kuint32 nDet = darray.nDet(); - if (m_fd < 0) { - sys_error (ERR_SEVERE, "Tried to write detector array with no file open [detarrayWrite]"); + fs.seekp (start_data); + if (! fs) { + sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]"); return false; } - if (lseek(m_fd, start_data, SEEK_SET) != start_data) { - sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]"); - return false; - } + fs.writeFloat64 (view_angle); + fs.writeInt32 (nDet); - if ( ! write_nfloat64 (&view_angle, m_fd)) { - sys_error (ERR_SEVERE, "Error writing detector array [detarrayWrite]"); - return false; - } - if ( ! write_nint32 (&nDet, m_fd)) { - sys_error (ERR_SEVERE, "Error writing detector array [detarrayWrite]"); - return false; - } for (int i = 0; i < nDet; i++) { - kfloat32 detval = detval_ptr[i]; - if ( ! write_nfloat32 (&detval, m_fd)) { - sys_error (ERR_SEVERE, "Error writing detector array [detarrayWrite]"); - return false; - } + kfloat32 detval = detval_ptr[i]; + fs.writeFloat32 (detval); } + if (! fs) + return (false); + return true; } - /* NAME * prt_projections Print projections data * diff --git a/libctsim/reconstr.cpp b/libctsim/reconstr.cpp index 11b9fc5..9c88780 100644 --- a/libctsim/reconstr.cpp +++ b/libctsim/reconstr.cpp @@ -11,7 +11,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: reconstr.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $ +** $Id: reconstr.cpp,v 1.2 2000/06/20 17:54:51 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 @@ -59,9 +59,8 @@ proj_reconst (ImageFile& im, Projections& proj, const FilterType filt_type, doub int nview = proj.nView(); double det_inc = proj.detInc(); double detlen = (ndet - 1) * det_inc; - double *vec_proj = new double [ndet]; // projection data int n_filtered_proj = ndet; - double* filtered_proj = new double [n_filtered_proj]; // convolved result + double filtered_proj [n_filtered_proj]; // convolved result #ifdef HAVE_BSPLINE_INTERP int spline_order = 0, zoom_factor = 0; @@ -77,12 +76,12 @@ proj_reconst (ImageFile& im, Projections& proj, const FilterType filt_type, doub double filterBW = 1. / det_inc; double filterMin = -detlen; double filterMax = detlen; - double* vec_filter = filter_generate (filt_type, filterBW, filterMin, filterMax, n_vec_filter, filt_param, D_SPATIAL, 0); + + SignalFilter filter (filt_type, filterBW, filterMin, filterMax, n_vec_filter, filt_param, D_SPATIAL, 0); #if HAVE_SGP - double *plot_xaxis = NULL; /* array for plotting */ SGP_ID gid; - plot_xaxis = new double [n_vec_filter]; + double plot_xaxis [n_vec_filter]; // array for plotting if (trace > TRACE_TEXT) { int i; @@ -91,7 +90,7 @@ proj_reconst (ImageFile& im, Projections& proj, const FilterType filt_type, doub for (i = 0, f = filterMin; i < n_vec_filter; i++, f += filterInc) plot_xaxis[i] = f; - gid = ezplot (plot_xaxis, vec_filter, n_vec_filter); + gid = ezplot (plot_xaxis, filter.getFilter(), n_vec_filter); cio_put_str ("Press any key to continue"); cio_kb_getc (); sgp2_close (gid); @@ -111,10 +110,7 @@ proj_reconst (ImageFile& im, Projections& proj, const FilterType filt_type, doub DetectorValue* detval = darray.detValues(); for (int j = 0; j < ndet; j++) - vec_proj[j] = detval[j]; - - for (int j = 0; j < ndet; j++) - filtered_proj[j] = convolve (vec_proj, vec_filter, det_inc, j, ndet, FUNC_BOTH); + filtered_proj[j] = filter.convolve (detval, det_inc, j, ndet, FUNC_BOTH); #ifdef HAVE_SGP if (trace >= TRACE_PLOT) { @@ -126,7 +122,7 @@ proj_reconst (ImageFile& im, Projections& proj, const FilterType filt_type, doub ezset ("box."); ezset ("grid."); ezset ("ufinish yes."); - ezplot (vec_proj, plot_xaxis, proj.nDet()); + ezplot (detval, plot_xaxis, proj.nDet()); ezset ("clear."); ezset ("xticks major 5."); ezset ("xlabel "); @@ -168,13 +164,6 @@ proj_reconst (ImageFile& im, Projections& proj, const FilterType filt_type, doub } delete bj; - delete vec_proj; - delete filtered_proj; - delete vec_filter; - -#ifdef HAVE_SGP - delete plot_xaxis; -#endif return (im); } diff --git a/libctsupport/Makefile.am b/libctsupport/Makefile.am index da467e6..df80a70 100644 --- a/libctsupport/Makefile.am +++ b/libctsupport/Makefile.am @@ -1,6 +1,6 @@ noinst_LIBRARIES = libctsupport.a INCLUDES=@my_includes@ -libctsupport_a_SOURCES= filefuncs.cpp strfuncs.cpp syserror.cpp timedate.cpp byteorder.cpp audio.cpp consoleio.cpp normangle.cpp simpson.cpp xform.cpp clip.cpp +libctsupport_a_SOURCES= filefuncs.cpp strfuncs.cpp syserror.cpp timedate.cpp fnetorderstream.cpp audio.cpp consoleio.cpp normangle.cpp simpson.cpp xform.cpp clip.cpp EXTRA_DIST=Makefile.nt diff --git a/libctsupport/byteorder.cpp b/libctsupport/byteorder.cpp deleted file mode 100644 index 6a2ddc9..0000000 --- a/libctsupport/byteorder.cpp +++ /dev/null @@ -1,343 +0,0 @@ -#if HAVE_CONFIG_H -#include "config.h" -#endif - -#if HAVE_UNISTD_H -#include -#endif - -#include "ctsupport.h" -#include "byteorder.h" - - -inline void -SwapBytes2 (void* buffer) -{ - unsigned char* p = static_cast(buffer); - unsigned char temp = p[0]; - p[0] = p[1]; - p[1] = temp; -} - -// 0<->3 1<->2 = 0123 -> 3210 -inline void -SwapBytes4 (void* buffer) -{ - unsigned char* p = static_cast(buffer); - unsigned char temp = p[0]; - p[0] = p[3]; - p[3] = temp; - temp = p[1]; - p[1] = p[2]; - p[2] = temp; -} - -// 0<->7 1<->6 2<->5 3<->4 = 01234567 -> 76543210 -inline void -SwapBytes8 (void* buffer) -{ - unsigned char* p = static_cast(buffer); - unsigned char temp = p[0]; - p[0] = p[7]; - p[7] = temp; - temp = p[1]; - p[1] = p[6]; - p[6] = temp; - temp = p[2]; - p[2] = p[5]; - p[5] = temp; - temp = p[3]; - p[3] = p[4]; - p[4] = temp; -} - -void -ConvertNetworkOrder (void* buffer, size_t bytes) -{ -#if ! defined (WORDS_BIGENDIAN) - if (bytes < 2) - return; - - char* start = static_cast(buffer); - char* end = start + bytes - 1; // last byte - size_t nSwap = bytes / 2; - - while (nSwap-- > 0) { - unsigned char c = *start; - *start++ = *end; - *end-- = c; - } -#endif -} - -void -ConvertReverseNetworkOrder (void* buffer, size_t bytes) -{ -#if defined (WORDS_BIGENDIAN) - if (bytes < 2) - return; - - char* start = static_cast(buffer); - char* end = start + bytes - 1; // last byte - size_t nSwap = bytes / 2; - - while (nSwap-- > 0) { - unsigned char c = *start; - *start++ = *end; - *end-- = c; - } -#endif -} - -bool -write_nint16 (kuint16 const *n_in, int fd) -{ - kuint16 n = *n_in; - -#ifndef WORDS_BIGENDIAN - SwapBytes2 (&n); -#endif - - if (write (fd, &n, 2) != 2) - return false; - else - return true; -} - -bool -read_nint16 (kuint16 *n_out, int fd) -{ - if (read (fd, n_out, 2) != 2) - return false; - -#ifndef WORDS_BIGENDIAN - SwapBytes2 (n_out); -#endif - - return true; -} - -bool -write_nint32 (kuint32 const *n_in, int fd) -{ - kuint32 n = *n_in; - -#ifndef WORDS_BIGENDIAN - SwapBytes4(&n); -#endif - - if (write (fd, &n, 4) != 4) - return false; - else - return true; -} - -bool -read_nint32 (kuint32 *n_out, int fd) -{ - if (read (fd, n_out, 4) != 4) - return false; - -#ifndef WORDS_BIGENDIAN - SwapBytes4 (n_out); -#endif - - return true; -} - -bool -write_nfloat32 (kfloat32 const *f_in, int fd) -{ - kfloat32 f = *f_in; - -#ifndef WORDS_BIGENDIAN - SwapBytes4 (&f); -#endif - - if (write (fd, &f, 4) != 4) - return false; - else - return true; -} - -bool -read_nfloat32 (kfloat32 *f_out, int fd) -{ - if (read (fd, f_out, 4) != 4) - return false; - -#ifndef WORDS_BIGENDIAN - SwapBytes4(f_out); -#endif - - return true; -} - -bool -write_nfloat64 (kfloat64 const *f_in, int fd) -{ - kfloat64 f = *f_in; - -#ifndef WORDS_BIGENDIAN - SwapBytes8 (&f); -#endif - - if (write (fd, &f, 8) != 8) - return false; - else - return true; -} - -bool -read_nfloat64 (kfloat64 *f_out, int fd) -{ - if (read (fd, f_out, 8) != 8) - return false; - -#ifndef WORDS_BIGENDIAN - SwapBytes8 (f_out); -#endif - - return true; -} - - - -onetorderstream& onetorderstream::writeInt16 (kuint16 n) { -#ifndef WORDS_BIGENDIAN - SwapBytes2 (&n); -#endif - write (&n, 2); - return (*this); -} - -inetorderstream& inetorderstream::readInt16 (kuint16& n) { - read (&n, 2); -#ifndef WORDS_BIGENDIAN - SwapBytes2 (&n); -#endif - return (*this); -} - -onetorderstream& onetorderstream::writeInt32 (kuint32 n) { -#ifndef WORDS_BIGENDIAN - SwapBytes4(&n); -#endif - write (&n, 4); - return (*this); -} - -inetorderstream& inetorderstream::readInt32 (kuint32& n) { - read (&n, 4); -#ifndef WORDS_BIGENDIAN - SwapBytes4 (&n); -#endif - return (*this); -} - -onetorderstream& onetorderstream::writeFloat32 (kfloat32 n) { -#ifndef WORDS_BIGENDIAN - SwapBytes4 (&n); -#endif - write (&n, 4); - return (*this); -} - -inetorderstream& inetorderstream::readFloat32 (kfloat32& n) { - read (&n, 4); -#ifndef WORDS_BIGENDIAN - SwapBytes4 (&n); -#endif - return (*this); -} - -onetorderstream& onetorderstream::writeFloat64 (kfloat64 n) { -#ifndef WORDS_BIGENDIAN - SwapBytes8 (&n); -#endif - write (&n, 8); - return (*this); -} - -inetorderstream& inetorderstream::readFloat64 (kfloat64& n) { - read (&n, 8); -#ifndef WORDS_BIGENDIAN - SwapBytes8 (&n); -#endif - return (*this); -} - - - -void -write_rnint16 (kuint16 n, ostream& ostr) -{ -#ifdef WORDS_BIGENDIAN - SwapBytes2 (&n); -#endif - ostr.write (&n, 2); -} - -void -read_rnint16 (kuint16& n, istream& istr) -{ - istr.read (&n, 2); -#ifdef WORDS_BIGENDIAN - SwapBytes2 (&n); -#endif -} - -void -write_rnint32 (kuint32 n, ostream& ostr) -{ -#ifdef WORDS_BIGENDIAN - SwapBytes4(&n); -#endif - ostr.write (&n, 4); -} - -void -read_rnint32 (kuint32& n, istream istr) -{ - istr.read (&n, 4); -#ifdef WORDS_BIGENDIAN - SwapBytes4 (&n); -#endif -} - -void -write_rnfloat32 (kfloat32 n, ostream ostr) -{ -#ifdef WORDS_BIGENDIAN - SwapBytes4 (&n); -#endif - ostr.write (&n, 4); -} - -void -read_rnfloat32 (kfloat32& n, istream istr) -{ - istr.read (&n, 4); -#ifdef WORDS_BIGENDIAN - SwapBytes4 (&n); -#endif -} - -void -write_rnfloat64 (kfloat64 n, ostream ostr) -{ -#ifdef WORDS_BIGENDIAN - SwapBytes8 (&n); -#endif - ostr.write (&n, 8); -} - -void -read_rnfloat64 (kfloat64& n, istream istr) -{ - istr.read (&n, 8); -#ifdef WORDS_BIGENDIAN - SwapBytes8 (&n); -#endif -} - diff --git a/libctsupport/fnetorderstream.cpp b/libctsupport/fnetorderstream.cpp new file mode 100644 index 0000000..7fff7a3 --- /dev/null +++ b/libctsupport/fnetorderstream.cpp @@ -0,0 +1,180 @@ +#if HAVE_CONFIG_H +#include "config.h" +#endif + +#if HAVE_UNISTD_H +#include +#endif + +#include "ctsupport.h" +#include "fnetorderstream.h" + + +void +ConvertNetworkOrder (void* buffer, size_t bytes) +{ +#ifndef WORDS_BIGENDIAN + if (bytes < 2) + return; + + char* start = static_cast(buffer); + char* end = start + bytes - 1; // last byte + size_t nSwap = bytes / 2; + + while (nSwap-- > 0) { + unsigned char c = *start; + *start++ = *end; + *end-- = c; + } +#endif +} + +void +ConvertReverseNetworkOrder (void* buffer, size_t bytes) +{ +#ifdef WORDS_BIGENDIAN + if (bytes < 2) + return; + + char* start = static_cast(buffer); + char* end = start + bytes - 1; // last byte + size_t nSwap = bytes / 2; + + while (nSwap-- > 0) { + unsigned char c = *start; + *start++ = *end; + *end-- = c; + } +#endif +} + +fnetorderstream& fnetorderstream::writeInt16 (kuint16 n) { +#ifndef WORDS_BIGENDIAN + SwapBytes2 (&n); +#endif + write (&n, 2); + return (*this); +} + +fnetorderstream& fnetorderstream::writeInt32 (kuint32 n) { +#ifndef WORDS_BIGENDIAN + SwapBytes4(&n); +#endif + write (&n, 4); + return (*this); +} + +fnetorderstream& fnetorderstream::writeFloat32 (kfloat32 n) { +#ifndef WORDS_BIGENDIAN + SwapBytes4 (&n); +#endif + write (&n, 4); + return (*this); +} + +fnetorderstream& fnetorderstream::writeFloat64 (kfloat64 n) { +#ifndef WORDS_BIGENDIAN + SwapBytes8 (&n); +#endif + write (&n, 8); + return (*this); +} + +fnetorderstream& fnetorderstream::readInt16 (kuint16& n) { + read (&n, 2); +#ifndef WORDS_BIGENDIAN + SwapBytes2 (&n); +#endif + return (*this); +} + +fnetorderstream& fnetorderstream::readInt32 (kuint32& n) { + read (&n, 4); +#ifndef WORDS_BIGENDIAN + SwapBytes4 (&n); +#endif + return (*this); +} + +fnetorderstream& fnetorderstream::readFloat32 (kfloat32& n) { + read (&n, 4); +#ifndef WORDS_BIGENDIAN + SwapBytes4 (&n); +#endif + return (*this); +} + +fnetorderstream& fnetorderstream::readFloat64 (kfloat64& n) { + read (&n, 8); +#ifndef WORDS_BIGENDIAN + SwapBytes8 (&n); +#endif + return (*this); +} + + + +frnetorderstream& frnetorderstream::writeInt16 (kuint16 n) { +#ifdef WORDS_BIGENDIAN + SwapBytes2 (&n); +#endif + write (&n, 2); + return (*this); +} + +frnetorderstream& frnetorderstream::writeInt32 (kuint32 n) { +#ifdef WORDS_BIGENDIAN + SwapBytes4(&n); +#endif + write (&n, 4); + return (*this); +} + +frnetorderstream& frnetorderstream::writeFloat32 (kfloat32 n) { +#ifdef WORDS_BIGENDIAN + SwapBytes4 (&n); +#endif + write (&n, 4); + return (*this); +} + +frnetorderstream& frnetorderstream::writeFloat64 (kfloat64 n) { +#ifdef WORDS_BIGENDIAN + SwapBytes8 (&n); +#endif + write (&n, 8); + return (*this); +} + +frnetorderstream& frnetorderstream::readInt16 (kuint16& n) { + read (&n, 2); +#ifdef WORDS_BIGENDIAN + SwapBytes2 (&n); +#endif + return (*this); +} + +frnetorderstream& frnetorderstream::readInt32 (kuint32& n) { + read (&n, 4); +#ifdef WORDS_BIGENDIAN + SwapBytes4 (&n); +#endif + return (*this); +} + +frnetorderstream& frnetorderstream::readFloat32 (kfloat32& n) { + read (&n, 4); +#ifdef WORDS_BIGENDIAN + SwapBytes4 (&n); +#endif + return (*this); +} + +frnetorderstream& frnetorderstream::readFloat64 (kfloat64& n) { + read (&n, 8); +#ifdef WORDS_BIGENDIAN + SwapBytes8 (&n); +#endif + return (*this); +} + diff --git a/src/sample-ctrec.sh.in b/src/sample-ctrec.sh.in index 11809a2..1a9f105 100755 --- a/src/sample-ctrec.sh.in +++ b/src/sample-ctrec.sh.in @@ -6,6 +6,11 @@ else bin="@prefix@/bin/" fi +if test "$1" = "clean" ; then + rm -f sample-phm.png sample-phm16.png sample-phm.if sample-pj.pj sample-pj.if sample-pj.png sample-pj16.png sample-rec.if sample-rec.png sample-rec16.png + exit +fi + # Generate phantom image ${bin}phm2if sample-phm.if 256 256 --nsample 2 --phantom herman -- 2.34.1