+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
--with-cgibin-dir=..., --with-cgibin-url=..., --with-webdata-dir=...,
--with-webdata-url=..., and --with-html-dir=... must be set.
+
CTSim Specific Configuration Help
---------------------------------
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
-------------
getopt_long appears broken, configure.in checks for cygwin to use
bundled version of getopt_long.
+
Basic Installation
==================
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
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
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.
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
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
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; }
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.
-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
+++ /dev/null
-#ifndef NETORDER_H
-#define NETORDER_H
-
-#include <iostream>
-
-/* 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
** 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
#endif
#include "ctsupport.h"
-#include "byteorder.h"
+#include "fnetorderstream.h"
#ifdef HAVE_SGP
#include "ezplot.h"
using namespace std;
#include "array2d.h"
+#include "fnetorderstream.h"
#include "imagefile.h"
#include "phantom.h"
#include "projections.h"
#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);
** 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
/* 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);
--- /dev/null
+#ifndef NETORDER_H
+#define NETORDER_H
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <fstream>
+#include <iostream>
+#include <string>
+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<unsigned char*>(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<unsigned char*>(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<unsigned char*>(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
#include <string>
#include <sys/types.h>
#include <unistd.h>
+#include <fstream>
+#include <iostream>
#include "ctsupport.h"
+#include "fnetorderstream.h"
#include "array2d.h"
using namespace std;
kuint16 headersize;
string filename;
int file_id;
- iostream *io;
+ frnetorderstream *m_pFS;
bool bHeaderWritten;
bool bDataWritten;
{
mPixelSize = sizeof(T);
signature = ('I' * 256 + 'F');
- file_id = -1;
mNX = 0;
mNY = 0;
headersize = 0;
mMinX = mMaxX = mMinY = mMaxY = 0;
mOffsetPV = 0;
mScalePV = 1;
- io = NULL;
+ m_pFS = NULL;
#if 0
const type_info& t_id = typeid(T);
void
Array2dFile<T>::fileClose (void)
{
- if (file_id >= 0) {
+ if (m_pFS) {
headerWrite ();
- close (file_id);
- file_id = -1;
+ m_pFS->close ();
+ m_pFS = NULL;
}
}
bool
Array2dFile<T>::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);
}
bool
Array2dFile<T>::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);
bool
Array2dFile<T>::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);
bool
Array2dFile<T>::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);
}
bool
Array2dFile<T>::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);
}
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);
}
bool
Array2dFile<T>::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);
}
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);
}
}
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);
}
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;
{
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<const void*>(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<class T>
#endif
};
-#define IMAGEFILE_64_BITS 1
+#undef IMAGEFILE_64_BITS
#ifdef IMAGEFILE_64_BITS
typedef F64Image ImageFile;
typedef kfloat64 ImageFileValue;
** 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
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;}
{ 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
bool headerRead (void);
bool headerWrite (void);
+ bool headerRead (fnetorderstream& fs);
+ bool headerWrite (fnetorderstream& fs);
void newProjData (void);
void deleteProjData (void);
** 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
*----------------------------------------------------------------------*/
+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)
{
** 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
break;
case S_REPLOT:
ez.i_plotimmediate = TRUE;
- ezplot (NULL, NULL, 0);
+ ezplot (static_cast<double*>(NULL), static_cast<double*>(NULL), 0);
#if 0
if (modeinteract == TRUE)
WAITKEY();
ez.o_unknowncurves = FALSE;
ez.o_reqcurves = ez.i_numcurves;
ez.i_plotimmediate = TRUE;
- ezplot (NULL, NULL, 0);
+ ezplot (static_cast<double*>(NULL), static_cast<double*>(NULL), 0);
ez.i_plotimmediate = FALSE;
}
}
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
+++ /dev/null
-/*****************************************************************************
-** 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);
-}
** 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
* 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;
}
* 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 {
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);
* 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.;
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);
* 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;
double b = PI / bw;
double b2 = TWOPI / bw;
- switch (filt_type) {
+ switch (fType) {
case FILTER_BANDLIMIT:
q = bw * sinc(u * w, 1.0);
break;
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);
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;
}
* 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)
*/
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);
}
+/* 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);
+}
+
** 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
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]);
}
** 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
*/
bool
-Projections::headerWrite (void)
+Projections::headerWrite (fnetorderstream& fs)
{
kuint32 _hsize = m_headerSize;
kuint32 _nView = m_nView;
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;
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;
}
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;
}
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;
}
*/
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 */;
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;
}
*/
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 */;
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
*
** 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
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;
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;
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);
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) {
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 ");
}
delete bj;
- delete vec_proj;
- delete filtered_proj;
- delete vec_filter;
-
-#ifdef HAVE_SGP
- delete plot_xaxis;
-#endif
return (im);
}
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
+++ /dev/null
-#if HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#if HAVE_UNISTD_H
-#include <unistd.h>
-#endif
-
-#include "ctsupport.h"
-#include "byteorder.h"
-
-
-inline void
-SwapBytes2 (void* buffer)
-{
- unsigned char* p = static_cast<unsigned char*>(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<unsigned char*>(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<unsigned char*>(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<char*>(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<char*>(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
-}
-
--- /dev/null
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#if HAVE_UNISTD_H
+#include <unistd.h>
+#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<char*>(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<char*>(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);
+}
+
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
projects / ctsim.git / commitdiff