--- /dev/null
+noinst_LIBRARIES = libctsim.a
+libir_a_SOURCES = filter.cpp scanner.cpp projections.cpp phantom.cpp options.cpp convolve.cpp reconstr.cpp phm2image.cpp imagefile.cpp reconstr.cpp phm2image.cpp backprojectors.cpp projections.cpp
+
+INCLUDES=@my_includes@
+EXTRA_DIST=Makefile.nt
+
--- /dev/null
+# Makefile for libir
+
+!include <ntwin32.mak>
+
+CC=cl
+LD=link
+CFLAGS=-O -nologo -I..\include
+LDFLAGS=
+O=.obj
+
+# variables
+OBJ1 = backproj$(O) bproj-trig$(O) mpi_ct$(O) raycollect$(O) imagefile$(O) bproj-d$(O) bproj-d2$(O) bproj-id2$(O) bproj-table$(O) convolve$(O) filter$(O) image$(O) options$(O) phm$(O) phm2image$(O) phmstd$(O) rayio$(O) reconstr$(O) scanner$(O)
+
+all: libir.lib
+
+.obj: .c
+ $(CC) -c $(cvarsdll) $(CFLAGS) $*.c
+
+.obj: .cpp
+ $(CC) -c $(cvarsdll) $(CFLAGS) $*.cpp
+
+
+libir.lib: $(OBJ1)
+ echo something to del > libir.lib
+ del libir.lib
+ lib /out:libir.lib $(OBJ1)
+
+
+clean:
+ echo dummy > a.obj
+ echo dummy > a.lib
+ echo dummy > a.exe
+ del *.obj
+ del *.exe
+ del *.lib
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: backprojectors.cpp Classes for backprojection
+** Programmer: Kevin Rosenberg
+** Date Started: June 2000
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: backprojectors.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+#include "ct.h"
+
+
+// FUNCTION IDENTIFICATION
+// Backproject* projector = selectBackprojector (...)
+//
+// PURPOSE
+// Selects a backprojector based on BackprojType
+// and initializes the backprojector
+
+Backproject* selectBackprojector (BackprojType bjType, const Projections& proj, ImageFile& im, InterpolationType interpType)
+{
+ Backproject* bj = NULL;
+
+ if (bjType == O_BPROJ_TRIG)
+ bj = static_cast<Backproject*>(new BackprojectTrig (proj, im, interpType));
+ else if (bjType == O_BPROJ_TABLE)
+ bj = static_cast<Backproject*>(new BackprojectTable (proj, im, interpType));
+ else if (bjType == O_BPROJ_DIFF)
+ bj = static_cast<Backproject*>(new BackprojectDiff (proj, im, interpType));
+ else if (bjType == O_BPROJ_DIFF2)
+ bj = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, interpType));
+ else if (bjType == O_BPROJ_IDIFF2)
+ bj = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, interpType));
+ else
+ sys_error (ERR_WARNING, "Illegal backproject type %d [selectBackprojector]");
+
+ return (bj);
+}
+
+
+// CLASS IDENTICATION
+// Backproject
+//
+// PURPOSE
+// Pure virtual base class for all backprojectors.
+
+Backproject::Backproject (const Projections& proj, ImageFile& im, const InterpolationType interpType)
+ : proj(proj), im(im), interpType(interpType)
+{
+ detInc = proj.detInc();
+ nDet = proj.nDet();
+ iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection
+ rotInc = proj.rotInc();
+
+ v = im.getArray();
+ nx = im.nx();
+ ny = im.ny();
+ im.arrayDataClear();
+
+ xMin = -proj.phmLen() / 2; // Retangular coords of phantom
+ xMax = xMin + proj.phmLen();
+ yMin = -proj.phmLen() / 2;
+ yMax = yMin + proj.phmLen();
+
+ xInc = (xMax - xMin) / nx; // size of cells
+ yInc = (yMax - yMin) / ny;
+
+ if (interpType != I_NEAREST && interpType != I_LINEAR)
+ sys_error (ERR_WARNING, "Illegal interpType %d [selectBackprojector]", interpType);
+}
+
+Backproject::~Backproject (void)
+{}
+
+void
+Backproject::ScaleImageByRotIncrement (void)
+{
+ for (int ix = 0; ix < nx; ix++)
+ for (int iy = 0; iy < ny; iy++)
+ v[ix][iy] *= rotInc;
+}
+
+void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
+{
+ printf ("r=%f, phi=%f\n", r, phi);
+ errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
+}
+
+void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
+{
+ printf ("ix=%d, iy=%d\n", ix, iy);
+ printf ("theta=%f, L=%f, detInc=%f\n", theta, L, detInc);
+ printf ("proj.ndet=%d, proj.detInc=%.4f, iDetCenter=%d\n", nDet, detInc, iDetCenter);
+ printf ("xMin=%15.8f, xMax=%15.8f, xInc=%15.8f\n", xMin, xMax, xInc);
+ printf ("yMin=%15.8f, yMax=%15.8f, yInc=%15.8f\n", yMin, yMax, yInc);
+ sys_error (ERR_WARNING, "iDetPos index outside bounds: %d [backprojector]", iDetPos);
+}
+
+
+// CLASS IDENTICATION
+// BackprojectTrig
+//
+// PURPOSE
+// Uses trigometric functions at each point in image for backprojection.
+
+void
+BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double theta = HALFPI + view_angle; // Add PI/2 to get perpendicular angle to detector
+ int ix, iy;
+ double x, y; // Rectang coords of center of pixel
+
+ for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
+ for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
+ double r = sqrt (x * x + y * y); // distance of cell from center
+ double phi = atan2 (y, x); // angle of cell from center
+ double L = r * cos (theta - phi); // position on detector
+
+ if (interpType == I_NEAREST) {
+ int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
+
+ if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
+ errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
+ else
+ v[ix][iy] += rotInc * filteredProj[iDetPos];
+ } else if (interpType == I_LINEAR) {
+ double p = L / detInc; // position along detector
+ double pFloor = floor (p);
+ int iDetPos = iDetCenter + static_cast<int>(pFloor);
+ double frac = p - pFloor; // fraction distance from det
+ if (iDetPos < 0 || iDetPos >= nDet - 1) // check for impossible: index outside of raysum pos
+ errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
+ else
+ v[ix][iy] += rotInc * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+ }
+ }
+}
+
+
+// CLASS IDENTICATION
+// BackprojectTable
+//
+// PURPOSE
+// Precalculates trigometric function value for each point in image for backprojection.
+
+BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, InterpolationType interpType)
+ : Backproject::Backproject (proj, im, interpType)
+{
+ arrayR.initSetSize (nx, ny);
+ arrayPhi.initSetSize (nx, ny);
+ r = arrayR.getArray();
+ phi = arrayPhi.getArray();
+
+ double x, y; // Rectang coords of center of pixel
+ int ix, iy;
+ for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
+ for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
+ r[ix][iy] = sqrt (x * x + y * y);
+ phi[ix][iy] = atan2 (y, x);
+ }
+}
+
+BackprojectTable::~BackprojectTable (void)
+{
+ ScaleImageByRotIncrement();
+}
+
+void
+BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double theta = HALFPI + view_angle; // add half PI to view angle to get perpendicular theta angle
+
+ for (int ix = 0; ix < nx; ix++) {
+ ImageFileColumn pImCol = v[ix];
+
+ for (int iy = 0; iy < ny; iy++) {
+ double L = r[ix][iy] * cos (theta - phi[ix][iy]);
+
+ if (interpType == I_NEAREST) {
+ int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
+
+ if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
+ errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
+ else
+ pImCol[iy] += filteredProj[iDetPos];
+ } else if (interpType == I_LINEAR) {
+ double dPos = L / detInc; // position along detector
+ double dPosFloor = floor (dPos);
+ int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
+ double frac = dPos - dPosFloor; // fraction distance from det
+ if (iDetPos < 0 || iDetPos >= nDet - 1)
+ errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
+ else
+ pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+ }
+ } // end for y
+ } // end for x
+}
+
+
+// CLASS IDENTICATION
+// BackprojectDiff
+//
+// PURPOSE
+// Backprojects by precalculating the change in L position for each x & y step in the image.
+// Iterates in x & y direction by adding difference in L position
+
+BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, InterpolationType interpType)
+ : Backproject::Backproject (proj, im, interpType)
+{
+ // calculate center of first pixel v[0][0]
+ double x = xMin + xInc / 2;
+ double y = yMin + yInc / 2;
+ start_r = sqrt (x * x + y * y);
+ start_phi = atan2 (y, x);
+
+ im.arrayDataClear();
+}
+
+BackprojectDiff::~BackprojectDiff()
+{
+ ScaleImageByRotIncrement();
+}
+
+void
+BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
+ double det_dx = xInc * sin (theta);
+ double det_dy = yInc * cos (theta);
+ double lColStart = start_r * cos (theta - start_phi); // calculate L for first point in image
+
+ for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
+ double curDetPos = lColStart;
+ ImageFileColumn pImCol = v[ix];
+
+ for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
+#ifdef DEBUG
+ printf ("[%2d,%2d]: %8.5lf ", ix, iy, curDetPos);
+#endif
+ if (interpType == I_NEAREST) {
+ int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc); // calc index in the filtered raysum vector
+
+ if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
+ errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
+ else
+ pImCol[iy] += filteredProj[iDetPos];
+ } else if (interpType == I_LINEAR) {
+ double detPos = curDetPos / detInc; // position along detector
+ double detPosFloor = floor (detPos);
+ int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
+ double frac = detPos - detPosFloor; // fraction distance from det
+ if (iDetPos < 0 || iDetPos >= nDet - 1)
+ errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
+ else
+ pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+ }
+ } // end for y
+ } // end for x
+}
+
+
+// CLASS IDENTICATION
+// BackprojectDiff2
+//
+// PURPOSE
+// Optimized version of BackprojectDiff
+
+void
+BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
+
+ // Distance betw. detectors for an angle given in units of detectors
+ double det_dx = xInc * sin (theta) / detInc;
+ double det_dy = yInc * cos (theta) / detInc;
+
+ // calculate detPosition for first point in image (ix=0, iy=0)
+ double detPosColStart = start_r * cos (theta - start_phi) / detInc;
+
+#ifdef DEBUG
+ printf ("start_r=%8.5f, start_phi=%8.5f, rotInc=%8.5f\n", start_r, start_phi, rotInc);
+#endif
+ for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
+ double curDetPos = detPosColStart;
+ ImageFileColumn pImCol = v[ix];
+
+ for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
+#ifdef DEBUG
+ printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(L))]);
+#endif
+ if (interpType == I_NEAREST) {
+ int iDetPos = iDetCenter + nearest<int> (curDetPos); // calc index in the filtered raysum vector
+
+ if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
+ errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
+ else
+ *pImCol++ += filteredProj[iDetPos];
+ } else if (interpType == I_LINEAR) {
+ double detPosFloor = floor (curDetPos);
+ int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
+ double frac = curDetPos - detPosFloor; // fraction distance from det
+ if (iDetPos < 0 || iDetPos >= nDet - 1)
+ errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
+ else
+ *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+ }
+ } // end for y
+ } // end for x
+}
+
+// CLASS IDENTICATION
+// BackprojectIntDiff2
+//
+// PURPOSE
+// Integer version of BackprojectDiff2
+
+void
+BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
+
+#if SIZEOF_LONG == 8
+ long int scale = 1 << 32;
+#else
+ long int scale = 1 << 16;
+#endif
+ double dScale = scale;
+ long int halfScale = scale / 2;
+
+ long int det_dx = nearest<long int> (xInc * sin (theta) / detInc * scale);
+ long int det_dy = nearest<long int> (yInc * cos (theta) / detInc * scale);
+
+ // calculate L for first point in image (0, 0)
+ long int detPosColStart = nearest<long int> (start_r * cos (theta - start_phi) / detInc * scale);
+
+ for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
+ long int curDetPos = detPosColStart;
+ ImageFileColumn pImCol = v[ix];
+
+ for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
+ if (interpType == I_NEAREST) {
+ int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
+ int iDetPos = iDetCenter + detPosNearest; // calc index in the filtered raysum vector
+
+ if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
+ errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
+ else
+ *pImCol++ += filteredProj[iDetPos];
+ } else if (interpType == I_LINEAR) {
+ long int detPosFloor = curDetPos / scale;
+ long int detPosRemainder = curDetPos % scale;
+ if (detPosRemainder < 0) {
+ detPosFloor--;
+ detPosRemainder += scale;
+ }
+ int iDetPos = iDetCenter + detPosFloor;
+ double frac = detPosRemainder / dScale;
+ if (iDetPos < 0 || iDetPos >= nDet - 1)
+ errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
+ else
+ *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+ }
+ } // end for y
+ } // end for x
+}
--- /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);
+}
--- /dev/null
+/*****************************************************************************
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: dialogs.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+#include "ct.h"
+
+
+/* NAME
+ * phm_add_pelem_kb Let user specify pelem, and add it to the pic
+ *
+ * SYNOPSIS
+ * phm_add_pelem_kb (phm)
+ * Phantom& pic PHANTOM that we are to add pelem to
+ *
+ * RETURNS
+ * true if user added an pelem to PHANTOM
+ * false if user decided not to add an pelem to the PHANTOM
+ */
+
+bool
+phm_add_pelem_kb (Phantom& phm)
+{
+ int retval = false;
+
+ do {
+ int pelemtype;
+
+ crt_clrline (1);
+ crt_clrline (2);
+ crt_set_cpos (1, 2);
+ crt_put_str ("1 = Rectangle, 2 = Triangle, 3 = Ellipse, 4 = Sector, 5 = Segment");
+ crt_set_cpos (1, 1);
+ crt_put_str ("Enter pelem type (1-5, 0 = no pelem) -- ");
+ scanf ("%d", &pelemtype);
+ crt_clrline (1);
+ crt_clrline (2);
+
+ if (pelemtype < 1) {
+ retval = false;
+ } else {
+ double cx, cy, u, v, rot, dens;
+
+ retval = true;
+ crt_set_cpos (1, 1);
+ crt_put_str ("Enter pelem specs (cx, cy, u, v, rot, density) -- ");
+ scanf ("%lf %lf %lf %lf %lf %lf %*c", &cx, &cy, &u, &v, &rot, &dens);
+ crt_clrline (1);
+ phm.addPelem (phm, pelemtype, cx, cy, u, v, rot, dens);
+ }
+ } while (retval == true);
+
+ return (true);
+}
+
+const Phantom& phm_select (Phantom& phm)
+{
+ string fname;
+ int phmnum;
+
+ printf ("Which phantom do you want to compile into pixels:\n");
+ printf (" 1 - Herman head phantom\n");
+ printf (" 2 - Rowland head phantom\n");
+ printf (" 3 - \n");
+ printf (" 4 - Rowland head phantom (Bordered)\n");
+ printf (" 6 - A Filter\n");
+ printf (" 7 - Unit pulse\n");
+ printf (" 8 - Enter PHANTOM from file\n");
+ printf (" 9 - Enter PHANTOM from keyboard\n");
+ printf ("Enter the number corresponding to your choice: ");
+ scanf ("%d%*c", &phmnum);
+
+ switch (phmnum) {
+ case 1:
+ phm.std_herman ();
+ break;
+ case 2:
+ phm.std_rowland ();
+ break;
+ case 4:
+ phm.std_rowland_bordered ();
+ break;
+ case 6:
+ phm.setComposition (P_FILTER);
+ break;
+ case 7:
+ phm.setComposition (P_UNIT_PULSE);
+ phm.addPelem (1, 0., 0., 100., 100., 0., 0.); /* outline */
+ phm.addPelem (3, 0., 0., 1., 1., 0., 1.); /* pulse */
+ break;
+ case 8:
+ printf ("Enter name of file: ");
+ scanf ("%s %*c", fname);
+ phm.setComposition (P_PELEMS);
+ if (phm.createFromFile (fname) == false)
+ cerr << "File " << fname << " doesn't contain valid Phantom declaration" << endl;
+ break;
+ case 9:
+ crt_clrscrn ();
+ phm_add_pelem_kb (phm);
+ crt_clrscrn ();
+ break;
+ default:
+ sys_error (ERR_FATAL, "Illegal Phantom number %d\n", phmnum);
+ break;
+ }
+
+ return (phm);
+}
+
+/* NAME
+ * interpolation_select Let user select an interpolation method
+ *
+ * SYNOPSIS
+ * interpolation_type = interpolation_select()
+ * int interpolation_type Method of interpolation to use
+ */
+
+int interpolation_select (void)
+{
+ bool got_it = false;
+
+ do {
+ int interp_type;
+
+ printf ("What interpolation method do you want to use:\n");
+ printf (" %2d - Nearest neighbor\n", I_NEAREST);
+ printf (" %2d - Linear\n", I_LINEAR);
+ printf (" %2d - B-Spline\n", I_BSPLINE);
+ printf ("Enter number corresponding to desired method: ");
+ scanf ("%d", &interp_type);
+ printf ("\n");
+
+ if (interp_name_of (interp_type) == NULL) {
+ cout << endl;
+ cio_beep ();
+ } else
+ got_it = true;
+
+ } while (! got_it);
+
+ return (interp_type);
+}
+
+
+/* NAME
+ * filter_select Let user select a filter
+ *
+ * SYNOPSIS
+ * filter_type = filt_select (filt_param)
+ * int filt_type Type of filter to use
+ * double *filt_param Returns parameter to filter
+ * Currently, only used with Hamming filters
+ */
+
+int
+filter_select (double *filt_param)
+{
+ bool got_it = false;
+
+ do {
+ int filt_type;
+
+ printf ("Which filter would you like to use:\n");
+ printf (" %2d - Bandlimiting\n", FILTER_BANDLIMIT);
+ printf (" %2d - Sinc\n", FILTER_SINC);
+ printf (" %2d - Hamming\n", FILTER_G_HAMMING);
+ printf (" %2d - Cosine\n", FILTER_COSINE);
+ printf (" %2d - Triangle\n", FILTER_TRIANGLE);
+ printf (" %2d - Abs * Bandlimit\n", FILTER_ABS_BANDLIMIT);
+ printf (" %2d - Abs * Sinc\n", FILTER_ABS_SINC);
+ printf (" %2d - Abs * Hamming \n", FILTER_ABS_G_HAMMING);
+ printf (" %2d - Abs * Cosine\n", FILTER_ABS_COSINE);
+ printf (" %2d - Shepp-Logan\n", FILTER_SHEPP);
+ printf ("Enter number corresponding to desired filter: ");
+ scanf ("%d", &filt_type);
+
+ if (filter_name_of (filt_type) == NULL) {
+ printf ("\n");
+ cio_beep ();
+ } else
+ got_it = true;
+
+ } while (! got_it);
+
+ if (filt_type == FILTER_G_HAMMING || filt_type == FILTER_ABS_G_HAMMING) {
+ cout << "Enter alpha (0-1): " << flush;
+ cin >> *filt_param;
+ } else
+ *filt_param = 0.0;
+
+ return (filt_type);
+}
+
+
+
+
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: filter.cpp
+** Purpose: Routines for signal-procesing filters
+** Progammer: Kevin Rosenberg
+** Date Started: Aug 1984
+**
+** 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 $
+**
+** 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
+ * filter_generate Generate a filter
+ *
+ * SYNOPSIS
+ * f = filter_generate (filt_type, bw, xmin, xmax, n, param, domain, analytic)
+ * double f Generated filter vector
+ * int filt_type Type of filter wanted
+ * double bw Bandwidth of filter
+ * double xmin, xmax Filter limits
+ * int n Number of points in filter
+ * double param General input parameter to filters
+ * int domain FREQ or SPATIAL domain wanted
+ * int numint Number if intervals for calculating
+ * discrete inverse fourier xform
+ * for spatial domain filters. For
+ * 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)
+{
+ double *f = new double [n];
+ double xinc = (xmax - xmin) / (n - 1);
+
+ if (filt_type == FILTER_SHEPP) {
+ double a = 2 * bw;
+ double c = - 4. / (a * a);
+ int center = (n - 1) / 2;
+ int sidelen = center;
+ f[center] = 4. / (a * a);
+
+ for (int i = 1; i <= sidelen; i++ )
+ f [center + i] = f [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);
+ } 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);
+ else
+ f[i] = filter_spatial_response_calc (filt_type, x, bw, param, numint);
+ } else {
+ sys_error (ERR_WARNING, "Illegal domain %d [filt_generate]", domain);
+ return (NULL);
+ }
+
+ return (f);
+}
+
+
+/* NAME
+ * filter_spatial_response_calc Calculate filter by discrete inverse fourier
+ * transform of filters's frequency
+ * response
+ *
+ * SYNOPSIS
+ * y = filter_spatial_response_calc (filt_type, x, 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 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)
+{
+ double zmin, zmax;
+
+ if (filt_type == FILTER_TRIANGLE) {
+ zmin = 0;
+ zmax = bw;
+ } else {
+ zmin = 0;
+ zmax = bw / 2;
+ }
+ double zinc = (zmax - zmin) / (n - 1);
+
+ 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);
+
+ double y = 2 * integrateSimpson (zmin, zmax, q, n);
+
+ return (y);
+}
+
+
+/* NAME
+ * filter_frequency_response Return filter frequency response
+ *
+ * SYNOPSIS
+ * h = filter_frequency_response (filt_type, u, 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 param General input parameter for various filters
+ */
+
+double
+filter_frequency_response (int filt_type, double u, double bw, double param)
+{
+ double q;
+ double au = fabs (u);
+
+ switch (filt_type) {
+ case FILTER_BANDLIMIT:
+ if (au >= bw / 2)
+ q = 0.;
+ else
+ q = 1;
+ break;
+ case FILTER_ABS_BANDLIMIT:
+ if (au >= bw / 2)
+ q = 0.;
+ else
+ q = au;
+ break;
+ case FILTER_TRIANGLE:
+ if (au >= bw)
+ q = 0;
+ else
+ q = 1 - au / bw;
+ break;
+ case FILTER_COSINE:
+ if (au >= bw / 2)
+ q = 0;
+ else
+ q = cos(PI * u / bw);
+ break;
+ case FILTER_ABS_COSINE:
+ if (au >= bw / 2)
+ q = 0;
+ else
+ q = au * cos(PI * u / bw);
+ break;
+ case FILTER_SINC:
+ q = bw * sinc (PI * bw * u, 1.);
+ break;
+ case FILTER_ABS_SINC:
+ q = au * bw * sinc (PI * bw * u, 1.);
+ break;
+ case FILTER_G_HAMMING:
+ if (au >= bw / 2)
+ q = 0;
+ else
+ q = param + (1 - param) * cos (TWOPI * u / bw);
+ break;
+ case FILTER_ABS_G_HAMMING:
+ if (au >= bw / 2)
+ q = 0;
+ else
+ q = au * (param + (1 - param) * cos(TWOPI * u / bw));
+ break;
+ default:
+ q = 0;
+ sys_error (ERR_WARNING,
+ "Frequency response for filter %d not implemented [filter_frequency_response]",
+ filt_type);
+ break;
+ }
+ return (q);
+}
+
+
+
+/* NAME
+ * filter_spatial_response_analytic Calculate filter by analytic inverse fourier
+ * transform of filters's frequency
+ * response
+ *
+ * SYNOPSIS
+ * y = filter_spatial_response_analytic (filt_type, x, 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 param General parameter for various filters
+ */
+
+double
+filter_spatial_response_analytic (int filt_type, double x, double bw, double param)
+{
+ double q, temp;
+ double u = TWOPI * x;
+ double w = bw / 2;
+ double b = PI / bw;
+ double b2 = TWOPI / bw;
+
+ switch (filt_type) {
+ case FILTER_BANDLIMIT:
+ q = bw * sinc(u * w, 1.0);
+ break;
+ case FILTER_TRIANGLE:
+ temp = sinc (u * w, 1.0);
+ q = bw * temp * temp;
+ break;
+ case FILTER_COSINE:
+ 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));
+ break;
+ case FILTER_ABS_BANDLIMIT:
+ q = 2 * integral_abscos (u, w);
+ break;
+ case FILTER_ABS_COSINE:
+ q = integral_abscos(b-u,w) + integral_abscos(b+u,w);
+ break;
+ case FILTER_ABS_G_HAMMING:
+ q = 2 * param * integral_abscos(u,w) +
+ (1-param)*(integral_abscos(u-b2,w)+integral_abscos(u+b2,w));
+ break;
+ case FILTER_SHEPP:
+ if (fabs (u) < 1E-7)
+ q = 4. / (PI * bw * bw);
+ else
+ q = fabs ((2 / bw) * sin (u * w)) * sinc (u * w, 1.) * sinc (u * w, 1.);
+ break;
+ case FILTER_SINC:
+ if (fabs (x) < bw / 2)
+ q = 1.;
+ else
+ q = 0.;
+ break;
+ case FILTER_ABS_SINC:
+ default:
+ sys_error (ERR_WARNING,
+ "Analytic filter type %d not implemented [filter_spatial_response_analytic]",
+ filt_type);
+ q = 0;
+ break;
+ }
+
+ return (q);
+}
+
+
+/* NAME
+ * sinc Return sin(x)/x function
+ *
+ * SYNOPSIS
+ * v = sinc (x, mult)
+ * double v sinc value
+ * double x, mult
+ *
+ * DESCRIPTION
+ * 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)
+ *
+ * SYNOPSIS
+ * q = integral_abscos (u, w)
+ * double q Integral value
+ * double u Integration variable
+ * double w Upper integration boundary
+ *
+ * DESCRIPTION
+ * Returns the value of integral of u*cos(u)*dV for V = 0 to w
+ */
+
+double
+integral_abscos (double u, double w)
+{
+ if (fabs (u) > F_EPSILON)
+ return (cos(u * w) - 1) / (u * u) + w / u * sin (u * w);
+ else
+ return (w * w / 2);
+}
+
+
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: imagefile.cpp
+** Purpose: Imagefile classes
+** Programmer: Kevin Rosenberg
+** Date Started: June 2000
+**
+** 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 $
+**
+** 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"
+
+
+///////////////////////////////////////////////////////////////////////////
+// CLASS IMPLEMENTATION
+//
+// Name: Array2dFileLabel
+// Purpose: Labels for Array2dFiles
+///////////////////////////////////////////////////////////////////////////
+
+void
+Array2dFileLabel::init (void)
+{
+ m_calcTime = 0;
+ m_labelType = L_EMPTY;
+ TIMEDATE td;
+ td_get_tmdt (&td);
+ m_year = td.d.year;
+ m_month = td.d.month;
+ m_day = td.d.date;
+ m_hour = td.t.hour;
+ m_minute = td.t.minute;
+ m_second = td.t.second;
+}
+
+Array2dFileLabel::Array2dFileLabel()
+{
+ init();
+}
+
+Array2dFileLabel::Array2dFileLabel(const char* const str, double ctime = 0.)
+ : m_strLabel (str)
+{
+ init();
+
+ m_labelType = L_USER;
+ m_calcTime = ctime;
+}
+
+Array2dFileLabel::Array2dFileLabel(const int type, const char* const str, double ctime = 0.)
+ : m_strLabel (str)
+{
+ init();
+
+ m_labelType = type;
+ m_calcTime = ctime;
+}
+
+Array2dFileLabel::~Array2dFileLabel()
+{
+}
+
+const string&
+Array2dFileLabel::getDateString (void) const
+{
+ ostringstream oss;
+ oss << m_month <<"/"<< m_day <<"/"<< m_year << " " << m_hour <<":"<< m_minute <<":"<< m_second;
+ m_strDate = oss.str();
+ return m_strDate;
+}
+
+
+/* FILE
+ * image.c Routines for managing images
+ *
+ * PROGRAMMER: Kevin Rosenberg
+ * DATE: Aug 1984
+ *
+ * FUNCTION
+ * Provides a set of routines for handling image files
+ */
+
+
+void
+image_filter_response (ImageFile& im, const DomainType domain, double bw, const FilterType filt_type, double filt_param, const int opt_trace)
+{
+ int hx = im.nx() / 2;
+ int hy = im.ny() / 2;
+ ImageFileArray v = im.getArray();
+
+ for (int i = -hx; i <= hx; i++) {
+ for (int j = -hy; j <= hy; j++) {
+ 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);
+ else
+ v[i+hx][j+hy] = filter_frequency_response (filt_type, r, bw, filt_param);
+ if (opt_trace >= TRACE_PHM)
+ printf ("r=%8.4f, v=%8.4f\n", r, v[i+hx][j+hy]);
+ }
+ }
+}
+
+int
+image_display (const ImageFile& im)
+{
+ ImageFileValue pmin, pmax;
+
+ im.getPixelValueRange (pmin, pmax);
+
+ return (image_display_scale (im, 1, pmin, pmax));
+}
+
+int image_display_scale (const ImageFile& im, const int scale, const double pmin, const double pmax)
+{
+ int grayscale[256];
+ int nx = im.nx();
+ int ny = im.ny();
+ ImageFileArray v = im.getArray();
+
+#if HAVE_G2_H
+ int* pens = new int [nx * ny * scale * scale ];
+
+ double view_scale = 255 / (pmax - pmin);
+ int id_X11 = g2_open_X11 (nx * scale, ny * scale);
+
+ for (int i = 0; i < 256; i++) {
+ double cval = i / 255.;
+ grayscale[i] = g2_ink (id_X11, cval, cval, cval);
+ }
+
+ for (int i= 0, iy = ny - 1; iy >= 0; iy--) {
+ for (int ix = 0; ix < nx; ix++) {
+ int cval = (int) ((v[ix][iy] - pmin) * view_scale);
+ if (cval < 0)
+ cval = 0;
+ else if (cval > 255)
+ cval = 255;
+ pens[i++] = grayscale[cval];
+ }
+ }
+
+ g2_image (id_X11, 0., 0., nx, ny, pens);
+
+ delete [] pens;
+ return (id_X11);
+
+#endif
+}
+
+
+
+
--- /dev/null
+/*****************************************************************************
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: options.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+#include "ct.h"
+
+int
+opt_set_trace (const char *optarg)
+{
+ int opt;
+
+ if (strcmp(optarg, O_TRACE_NONE_STR) == 0)
+ opt = TRACE_NONE;
+ else if (strcmp(optarg, O_TRACE_TEXT_STR) == 0)
+ opt = TRACE_TEXT;
+ else if (strcmp(optarg, O_TRACE_PHM_STR) == 0)
+ opt = TRACE_PHM;
+ else if (strcmp(optarg, O_TRACE_PLOT_STR) == 0)
+ opt = TRACE_PLOT;
+ else if (strcmp(optarg, O_TRACE_CLIPPING_STR) == 0)
+ opt = TRACE_CLIPPING;
+ else if (strcmp(optarg, O_TRACE_RAYS_STR) == 0)
+ opt = TRACE_RAYS;
+ else {
+ sys_error(ERR_WARNING,"Invalid trace option %s\n", optarg);
+ opt = -1;
+ }
+
+ return (opt);
+}
+
+const char *
+name_of_phantom (const int phmnum)
+{
+ const char *name = "Unknown phantom";
+
+ if (phmnum == O_PHM_HERMAN)
+ name = O_PHM_HERMAN_STR;
+ else if (phmnum == O_PHM_ROWLAND)
+ name = O_PHM_ROWLAND_STR;
+ else if (phmnum == O_PHM_BROWLAND)
+ name = O_PHM_BROWLAND_STR;
+ else if (phmnum == O_PHM_UNITPULSE)
+ name = O_PHM_UNITPULSE_STR;
+
+ return (name);
+}
+
+int
+opt_set_phantom (const char *optarg)
+{
+ int opt;
+
+ if (strcmp(optarg, O_PHM_HERMAN_STR) == 0)
+ opt = O_PHM_HERMAN;
+ else if (strcmp(optarg, O_PHM_ROWLAND_STR) == 0)
+ opt = O_PHM_ROWLAND;
+ else if (strcmp(optarg, O_PHM_BROWLAND_STR) == 0)
+ opt = O_PHM_BROWLAND;
+ else if (strcmp(optarg, O_PHM_UNITPULSE_STR) == 0)
+ opt = O_PHM_UNITPULSE;
+ else {
+ sys_error(ERR_WARNING,"Invalid phantom option %s\n", optarg);
+ opt = -1;
+ }
+
+ return (opt);
+
+}
+
+InterpolationType
+opt_set_interpolation (const char *optarg)
+{
+ InterpolationType opt;
+
+ if (strcmp(optarg, O_INTERP_NEAREST_STR) == 0)
+ opt = I_NEAREST;
+ else if (strcmp(optarg, O_INTERP_LINEAR_STR) == 0)
+ opt = I_LINEAR;
+#if HAVE_BSPLINE_INTERP
+ else if (strcmp(optarg, O_INTERP_BSPLINE_STR) == 0)
+ opt = I_BSPLINE;
+#endif
+ else {
+ sys_error(ERR_WARNING, "Invalid interpolation type %s\n", optarg);
+ opt = static_cast<InterpolationType>(-1);
+ }
+
+ return (opt);
+}
+
+FilterType
+opt_set_filter (const char *optarg)
+{
+ FilterType opt;
+
+ if (strcmp(optarg, O_FILTER_BANDLIMIT_STR) == 0)
+ opt = FILTER_BANDLIMIT;
+ else if (strcmp(optarg, O_FILTER_HAMMING_STR) == 0)
+ opt = FILTER_G_HAMMING;
+ else if (strcmp(optarg, O_FILTER_SINC_STR) == 0)
+ opt = FILTER_SINC;
+ else if (strcmp(optarg, O_FILTER_COS_STR) == 0)
+ opt = FILTER_COSINE;
+ else if (strcmp(optarg, O_FILTER_TRIANGLE_STR) == 0)
+ opt = FILTER_TRIANGLE;
+ else if (strcmp(optarg, O_FILTER_ABS_BANDLIMIT_STR) == 0)
+ opt = FILTER_ABS_BANDLIMIT;
+ else if (strcmp(optarg, O_FILTER_ABS_HAMMING_STR) == 0)
+ opt = FILTER_ABS_G_HAMMING;
+ else if (strcmp(optarg, O_FILTER_ABS_SINC_STR) == 0)
+ opt = FILTER_ABS_SINC;
+ else if (strcmp(optarg, O_FILTER_ABS_COS_STR) == 0)
+ opt = FILTER_ABS_COSINE;
+ else if (strcmp(optarg, O_FILTER_SHEPP_STR) == 0)
+ opt = FILTER_SHEPP;
+ else {
+ sys_error(ERR_WARNING, "Invalid filter type %s\n", optarg);
+ opt = static_cast<FilterType>(-1);
+ }
+
+ return (opt);
+}
+
+const char *
+name_of_filter (const int filter)
+{
+ const char *name = "Unknown filter";
+
+ if (filter == FILTER_SHEPP)
+ name = O_FILTER_SHEPP_STR;
+ else if (filter == FILTER_ABS_COSINE)
+ name = O_FILTER_ABS_COS_STR;
+ else if (filter == FILTER_ABS_SINC)
+ name = O_FILTER_ABS_SINC_STR;
+ else if (filter == FILTER_ABS_G_HAMMING)
+ name = O_FILTER_ABS_HAMMING_STR;
+ else if (filter == FILTER_ABS_BANDLIMIT)
+ name = O_FILTER_ABS_BANDLIMIT_STR;
+ else if (filter == FILTER_COSINE)
+ name = O_FILTER_COS_STR;
+ else if (filter == FILTER_SINC)
+ name = O_FILTER_SINC_STR;
+ else if (filter == FILTER_G_HAMMING)
+ name = O_FILTER_HAMMING_STR;
+ else if (filter == FILTER_BANDLIMIT)
+ name = O_FILTER_BANDLIMIT_STR;
+ else if (filter == FILTER_TRIANGLE)
+ name = O_FILTER_TRIANGLE_STR;
+
+ return (name);
+}
+
+DomainType
+opt_set_filter_domain (const char *optarg)
+{
+ DomainType opt;
+
+ if (strcmp(optarg, D_SPATIAL_STR) == 0)
+ opt = D_SPATIAL;
+ else if (strcmp(optarg, D_FREQ_STR) == 0)
+ opt = D_FREQ;
+ else {
+ sys_error(ERR_WARNING, "Invalid filter domain %s\n", optarg);
+ opt = static_cast<DomainType>(-1);
+ }
+
+ return (opt);
+}
+
+const char *
+name_of_filter_domain (const DomainType domain)
+{
+ const char *name = "Unknown domain";
+
+ if (domain == D_SPATIAL)
+ return(D_SPATIAL_STR);
+ else if (domain == D_FREQ)
+ return(D_FREQ_STR);
+
+ return (name);
+}
+
+
+BackprojType
+opt_set_backproj (const char *optarg)
+{
+ BackprojType opt;
+
+ if (strcmp(optarg, O_BPROJ_TRIG_STR) == 0)
+ opt = O_BPROJ_TRIG;
+ else if (strcmp(optarg, O_BPROJ_TABLE_STR) == 0)
+ opt = O_BPROJ_TABLE;
+ else if (strcmp(optarg, O_BPROJ_DIFF_STR) == 0)
+ opt = O_BPROJ_DIFF;
+ else if (strcmp(optarg, O_BPROJ_DIFF2_STR) == 0)
+ opt = O_BPROJ_DIFF2;
+ else if (strcmp(optarg, O_BPROJ_IDIFF2_STR) == 0)
+ opt = O_BPROJ_IDIFF2;
+ else {
+ sys_error(ERR_WARNING, "Invalid backprojection method %s\n", optarg);
+ opt = static_cast<BackprojType>(-1);
+ }
+
+ return (opt);
+}
+
+const char *
+name_of_backproj(const BackprojType bproj)
+{
+ const char *name = "Unknown backprojection method";
+
+ if (bproj == O_BPROJ_TRIG)
+ name = O_BPROJ_TRIG_STR;
+ else if (bproj == O_BPROJ_TABLE)
+ name = O_BPROJ_TABLE_STR;
+ else if (bproj == O_BPROJ_DIFF)
+ name = O_BPROJ_DIFF_STR;
+ else if (bproj == O_BPROJ_DIFF2)
+ name = O_BPROJ_DIFF2_STR;
+ else if (bproj == O_BPROJ_IDIFF2)
+ name = O_BPROJ_IDIFF2_STR;
+
+ return (name);
+}
+
+
+
+/* NAME
+ * name_of_interp Return name of interpolation method
+ *
+ * SYNOPSIS
+ * name = name_of_interp (interp_type)
+ * char *name Name of interpolation method
+ * int interp_type Method of interpolation
+ *
+ * NOTES
+ * Returns NULL if interp_type is invalid
+ */
+
+const char *
+name_of_interpolation (int interp_type)
+{
+ if (interp_type == I_NEAREST)
+ return (O_INTERP_NEAREST_STR);
+ else if (interp_type == I_LINEAR)
+ return (O_INTERP_LINEAR_STR);
+#if HAVE_BSPLINE_INTERP
+ else if (interp_type == I_BSPLINE)
+ return (O_INTERP_BSPLINE_STR);
+#endif
+ else
+ return ("Unknown interpolation method");
+}
+
+
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: phm.cpp
+** Purpose: Routines for phantom objects
+** Progammer: Kevin Rosenberg
+** Date Started: Aug 1984
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: phantom.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+#include "ct.h"
+
+
+// CLASS IDENTIFICATION
+// Phanton
+//
+
+Phantom::Phantom (void)
+{
+ m_nPElem = 0;
+ m_xmin = 1E30;
+ m_xmax = -1E30;
+ m_ymin = 1E30;
+ m_ymax = -1E30;
+ m_diameter = 0;
+ m_composition = P_PELEMS;
+}
+
+
+Phantom::~Phantom (void)
+{
+ for (PElemIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {
+ delete *i;
+ }
+}
+
+
+void
+Phantom::create (const int phmid)
+{
+ switch (phmid)
+ {
+ case O_PHM_HERMAN:
+ std_herman();
+ break;
+ case O_PHM_ROWLAND:
+ std_rowland();
+ break;
+ case O_PHM_BROWLAND:
+ std_rowland_bordered ();
+ break;
+ case O_PHM_UNITPULSE:
+ m_composition = P_UNIT_PULSE;
+ addPElem ("rectangle", 0., 0., 100., 100., 0., 0.); // outline
+ addPElem ("ellipse", 0., 0., 1., 1., 0., 1.); // pulse
+ break;
+ default:
+ sys_error (ERR_WARNING, "Illegal phantom id %d\n", phmid);
+ break;
+ }
+}
+
+
+/* METHOD IDENTIFICATION
+ * createFromFile Add PhantomElements from file
+ *
+ * SYNOPSIS
+ * createFromFile (filename)
+ *
+ * RETURNS
+ * true if pelem were added
+ * false if an pelem not added
+ */
+
+bool
+Phantom::createFromFile (const char* const fname)
+{
+ bool stoploop = false;
+ bool retval = false;
+ FILE *fp;
+
+ if ((fp = fopen (fname, "r")) == NULL)
+ return (false);
+
+ do {
+ double cx, cy, u, v, rot, dens;
+ char pelemtype[80];
+ int n = fscanf (fp, "%79s %lf %lf %lf %lf %lf %lf",
+ pelemtype, &cx, &cy, &u, &v, &rot, &dens);
+
+ if (n == EOF || n == 0) { /* end of file */
+ stoploop = true;
+ retval = false;
+ } else if (n != 7) {
+ stoploop = true;
+ retval = false;
+ } else {
+ addPElem (pelemtype, cx, cy, u, v, rot, dens);
+ retval = true;
+ }
+ } while (stoploop == false);
+
+ fclose (fp);
+
+ return (retval);
+}
+
+
+/* NAME
+ * addPElem Add pelem
+ *
+ * SYNOPSIS
+ * addPElem (type, cx, cy, u, v, rot, atten)
+ * char *type type of pelem (box, ellipse, etc)
+ * double cx, cy pelem center
+ * double u,v pelem size
+ * double rot rotation angle of pelem (in degrees)
+ * double atten x-ray attenuation cooefficient
+ */
+
+void
+Phantom::addPElem (const char *type, const double cx, const double cy, const double u, const double v, const double rot, const double atten)
+{
+ PhantomElement *pelem = new PhantomElement (type, cx, cy, u, v, rot, atten);
+
+ m_listPElem.push_front (pelem);
+
+ // update phantom limits
+ if (m_xmin > pelem->xmin()) m_xmin = pelem->xmin();
+ if (m_xmax < pelem->xmax()) m_xmax = pelem->xmax();
+ if (m_ymin > pelem->ymin()) m_ymin = pelem->ymin();
+ if (m_ymax < pelem->ymax()) m_ymax = pelem->ymax();
+
+ if (m_diameter < pelem->diameter())
+ m_diameter = pelem->diameter();
+
+ // m_diameter = lineLength(m_xmin, m_ymin, m_xmax, m_ymax);
+
+ m_nPElem++;
+}
+
+
+/*----------------------------------------------------------------------*/
+/* Input-Output Routines */
+/*----------------------------------------------------------------------*/
+
+
+/* NAME
+ * print Print vertices of Phantom pelems
+ *
+ * SYNOPSIS
+ * print (phm)
+ */
+
+void
+Phantom::print (void) const
+{
+ printf("PRINTING Phantom\n\n");
+ printf("number of pelems in Phantom = %d\n", m_nPElem);
+ printf("limits: xmin=%8.2g ymin=%8.2g xmax=%8.2g ymax=%8.2g\n",
+ m_xmin, m_ymin, m_xmax, m_ymax);
+
+ for (PElemIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {
+ printf("PELEM:\n");
+ printf("# pts=%3d atten = %7.4f rot = %7.2f (deg)\n",
+ (*i)->nOutlinePoints(), (*i)->atten(), convertRadiansToDegrees ((*i)->rot()));
+
+ printf("xmin=%7.3g ymin=%7.3g xmax=%7.3g ymax=%7.3g\n",
+ (*i)->xmin(), (*i)->ymin(), (*i)->xmax(), (*i)->ymax());
+
+ // for (int i = 0; i < m_nPoints; i++)
+ // printf("\t%8.3g %8.3g\n", i->xOutline()[i], i->yOutline()[i]);
+ }
+}
+
+
+/* NAME
+ * show Show vector outline of Phantom to user
+ *
+ * SYNOPSIS
+ * show (pic)
+ */
+
+#ifdef HAVE_SGP
+void
+Phantom::show (void) const
+{
+ double wsize = m_xmax - m_xmin;
+ double xmin = m_xmin;
+ double ymin = m_ymin;
+ double xmax, ymax;
+ SGP_ID gid;
+
+ if ((m_ymax - m_ymin) > wsize)
+ wsize = m_ymax - m_ymin;
+ wsize *= 1.1;
+
+ xmax = xmin + wsize;
+ ymax = ymin + wsize;
+
+ printf("Drawing Phantom:\n\n");
+ printf(" data limits: %9.3g, %9.3g, %9.3g, %9.3g\n",
+ m_xmin, m_ymin, m_xmax, m_ymax);
+ printf(" window size: %9.3g, %9.3g, %9.3g, %9.3g\n",
+ xmin, ymin, xmax, ymax);
+
+ gid = sgp2_init(0, 0, "Phantom Show");
+ sgp2_window (xmin, ymin, xmax, ymax);
+
+ draw();
+
+ termgrf2();
+}
+#endif
+
+
+/* NAME
+ * draw Draw vector outline of Phantom
+ *
+ * SYNOPSIS
+ * draw ()
+ */
+
+#ifdef HAVE_SGP
+void
+Phantom::draw (void) const
+{
+ for (PElemIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++)
+ sgp2_polyline_abs ((*i)->xOutline(), (*i)->yOutline(), (*i)->nOutlinePoints());
+}
+#endif
+
+
+/* NAME
+ * std_rowland Make head phantom of S.W. Rowland
+ *
+ * SYNOPSIS
+ * std_rowland ()
+ *
+ * REFERENCES
+ * S. W. Rowland, "Computer Implementation of Image Reconstruction
+ * Formulas", in "Image Reconstruction from Projections: Implementation
+ * and Applications", edited by G. T. Herman, 1978.
+ */
+
+void
+Phantom::std_rowland (void)
+{
+ addPElem("ellipse", 0.0000, 0.0000, 0.6900, 0.9200, 0.0, 1.00);
+ addPElem("ellipse", 0.0000, -0.0184, 0.6624, 0.8740, 0.0, -0.98);
+ addPElem("ellipse", 0.2200, 0.0000, 0.1100, 0.3100, -18.0, -0.02);
+ addPElem("ellipse", -0.2200, 0.0000, 0.1600, 0.4100, 18.0, -0.02);
+ addPElem("ellipse", 0.0000, 0.3500, 0.2100, 0.2500, 0.0, 0.01);
+ addPElem("ellipse", 0.0000, 0.1000, 0.0460, 0.0460, 0.0, 0.01);
+ addPElem("ellipse", 0.0000, -0.1000, 0.0460, 0.0460, 0.0, 0.01);
+ addPElem("ellipse", -0.0800, -0.6050, 0.0460, 0.0230, 0.0, 0.01);
+ addPElem("ellipse", 0.0000, -0.6050, 0.0230, 0.0230, 0.0, 0.01);
+ addPElem("ellipse", 0.0600, -0.6050, 0.0230, 0.0230, 0.0, 0.01);
+ addPElem("ellipse", 0.5538, -0.3858, 0.0330, 0.2060, -18.0, 0.03);
+}
+
+void
+Phantom::std_rowland_bordered (void)
+{
+ std_rowland ();
+ addPElem ("ellipse", 0.000, 0.0000, 0.7500, 1.000, 0.0, 0.00);
+}
+
+/* NAME
+ * std_herman Standard head phantom of G. T. Herman
+ *
+ * SYNOPSIS
+ * std_herman ()
+ *
+ * REFERENCES
+ * G. T. Herman, "Image Reconstructions from Projections: The Fundementals
+ * of Computed Tomography", 1979.
+ */
+
+void
+Phantom::std_herman (void)
+{
+ addPElem("ellipse", 0.000, 1.50, 0.375, 0.3000, 90.00, -0.003);
+ addPElem("ellipse", 0.675, -0.75, 0.225, 0.1500, 140.00, 0.010);
+ addPElem("ellipse", 0.750, 1.50, 0.375, 0.2250, 50.00, 0.003);
+ addPElem("segment", 1.375, -7.50, 1.100, 0.6250, 19.20, -0.204);
+ addPElem("segment", 1.375, -7.50, 1.100, 4.3200, 19.21, 0.204);
+ addPElem("segment", 0.000, -2.25, 1.125, 0.3750, 0.00, -0.003);
+ addPElem("segment", 0.000, -2.25, 1.125, 3.0000, 0.00, 0.003);
+ addPElem("segment", -1.000, 3.75, 1.000, 0.5000, 135.00, -0.003);
+ addPElem("segment", -1.000, 3.75, 1.000, 3.0000, 135.00, 0.003);
+ addPElem("segment", 1.000, 3.75, 1.000, 0.5000, 225.00, -0.003);
+ addPElem("segment", 1.000, 3.75, 1.000, 3.0000, 225.00, 0.003);
+ addPElem("triangle", 5.025, 3.75, 1.125, 0.5000, 110.75, 0.206);
+ addPElem("triangle",-5.025, 3.75, 1.125, 0.9000,-110.75, 0.206);
+ addPElem("ellipse", 0.000, 0.00, 8.625, 6.4687, 90.00, 0.416);
+ addPElem("ellipse", 0.000, 0.00, 7.875, 5.7187, 90.00, -0.206);
+}
+
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////
+// CLASS IDENTIFICATION
+//
+// PhantomElement
+//
+// PURPOSE
+//
+////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+PhantomElement::PhantomElement (const char *type, const double cx, const double cy, const double u, const double v, const double rot, const double atten)
+ : m_cx(cx), m_cy(cy), m_u(u), m_v(v), m_atten(atten), m_nPoints(0), m_xOutline(0), m_yOutline(0)
+{
+ m_rot = convertDegreesToRadians (rot); // convert angle to radians
+
+ if (strcasecmp (type, "rectangle") == 0)
+ m_type = PELEM_RECTANGLE;
+ else if (strcasecmp (type, "triangle") == 0)
+ m_type = PELEM_TRIANGLE;
+ else if (strcasecmp (type, "ellipse") == 0)
+ m_type = PELEM_ELLIPSE;
+ else if (strcasecmp (type, "sector") == 0)
+ m_type = PELEM_SECTOR;
+ else if (strcasecmp (type, "segment") == 0)
+ m_type = PELEM_SEGMENT;
+ else {
+ sys_error (ERR_WARNING, "Unknown PhantomElement type %s [PhantomElement::PhantomElement]", type);
+ m_type = PELEM_INVALID;
+ }
+
+ makeTransformMatrices (); // calc transform matrices between phantom and normalized phantomelement
+ makeVectorOutline (); // calculate vector outline of pelem
+
+ // Find maximum diameter of Object
+ double r2Max = 0;
+ for (int i = 0; i < m_nPoints; i++) {
+ double r2 = (m_xOutline[i] * m_xOutline[i]) + (m_yOutline[i] * m_yOutline[i]);
+ if (r2 > r2Max)
+ r2Max = r2;
+ }
+ m_diameter = 2 * sqrt( r2Max );
+
+ m_rectLimits[0] = m_xmin; m_rectLimits[1] = m_ymin;
+ m_rectLimits[2] = m_xmax; m_rectLimits[3] = m_ymax;
+}
+
+
+PhantomElement::~PhantomElement (void)
+{
+ delete m_xOutline;
+ delete m_yOutline;
+}
+
+void
+PhantomElement::makeTransformMatrices (void)
+{
+ GRFMTX_2D temp;
+
+ // To map normalized Pelem coords to world Phantom
+ // scale by (u, v)
+ // rotate by rot
+ // translate by (cx, cy)
+
+ scale_mtx2 (m_xformObjToPhm, m_u, m_v);
+ rot_mtx2 (temp, m_rot);
+ mult_mtx2 (m_xformObjToPhm, temp, m_xformObjToPhm);
+ xlat_mtx2 (temp, m_cx, m_cy);
+ mult_mtx2 (m_xformObjToPhm, temp, m_xformObjToPhm);
+
+ // to map world Phantom coodinates to normalized PElem coords
+ // translate by (-cx, -cy)
+ // rotate by -rot
+ // scale by (1/u, 1/v)
+
+ xlat_mtx2 (m_xformPhmToObj, -m_cx, -m_cy);
+ rot_mtx2 (temp, -m_rot);
+ mult_mtx2 (m_xformPhmToObj, temp, m_xformPhmToObj);
+ scale_mtx2 (temp, 1 / m_u, 1 / m_v);
+ mult_mtx2 (m_xformPhmToObj, temp, m_xformPhmToObj);
+}
+
+
+/* NAME
+ * pelem_make_points INTERNAL routine to calculate point array for an pelem
+ *
+ * SYNOPSIS
+ * makepelempts (pelem)
+ * PELEM *pelem pelem whose points we are calculating
+ *
+ * NOTES
+ * Called by phm_add_pelem()
+ */
+
+static const double SCALE_PELEM_EXTENT=0.005; // increase pelem limits by 0.5%
+
+void
+PhantomElement::makeVectorOutline (void)
+{
+ double radius, theta, start, stop;
+ double xfact, yfact;
+ int cpts;
+
+ m_nPoints = 0;
+ switch (m_type) {
+ case PELEM_RECTANGLE:
+ m_nPoints = 5;
+ m_xOutline = new double [m_nPoints];
+ m_yOutline = new double [m_nPoints];
+ m_xOutline[0] =-m_u; m_yOutline[0] =-m_v;
+ m_xOutline[1] = m_u; m_yOutline[1] =-m_v;
+ m_xOutline[2] = m_u; m_yOutline[2] = m_v;
+ m_xOutline[3] =-m_u; m_yOutline[3] = m_v;
+ m_xOutline[4] =-m_u; m_yOutline[4] =-m_v;
+ break;
+ case PELEM_TRIANGLE:
+ m_nPoints = 4;
+ m_xOutline = new double [m_nPoints];
+ m_yOutline = new double [m_nPoints];
+ m_xOutline[0] =-m_u; m_yOutline[0] = 0.0;
+ m_xOutline[1] = m_u; m_yOutline[1] = 0.0;
+ m_xOutline[2] = 0.0; m_yOutline[2] = m_v;
+ m_xOutline[3] =-m_u; m_yOutline[3] = 0.0;
+ break;
+ case PELEM_ELLIPSE:
+ cpts = numCirclePoints (TWOPI);
+ m_nPoints = cpts;
+ m_xOutline = new double [m_nPoints];
+ m_yOutline = new double [m_nPoints];
+ calcEllipsePoints (m_xOutline, m_yOutline, cpts, m_u, m_v);
+ break;
+ case PELEM_SECTOR:
+ radius = sqrt(m_u * m_u + m_v * m_v);
+ theta = atan(m_u / m_v); // angle with y-axis
+ start = 3.0 * HALFPI - theta;
+ stop = 3.0 * HALFPI + theta;
+ cpts = numCirclePoints (stop - start);
+ m_nPoints = 3 + cpts;
+ m_xOutline = new double [m_nPoints];
+ m_yOutline = new double [m_nPoints];
+
+ m_xOutline[0] = 0.0; m_yOutline[0] = m_v;
+ m_xOutline[1] =-m_u; m_yOutline[1] = 0.0;
+ calcArcPoints (&m_xOutline[2], &m_yOutline[2], cpts, 0.0, m_v, radius, start, stop);
+ m_xOutline[cpts + 2] = 0.0;
+ m_yOutline[cpts + 2] = m_v;
+ break;
+ case PELEM_SEGMENT:
+ radius = sqrt(m_u * m_u + m_v * m_v);
+ theta = atan (m_u / m_v); // angle with y-axis
+ start = 3.0 * HALFPI - theta;
+ stop = 3.0 * HALFPI + theta;
+
+ cpts = numCirclePoints (stop - start);
+ m_nPoints = cpts + 1;
+ m_xOutline = new double [m_nPoints];
+ m_yOutline = new double [m_nPoints];
+
+ calcArcPoints (m_xOutline, m_yOutline, cpts, 0.0, m_v, radius, start, stop);
+ m_xOutline[cpts] = -m_u;
+ m_yOutline[cpts] = 0.0;
+ break;
+ default:
+ sys_error(ERR_WARNING, "illegal pelem type %d [makeVectorOutline]", m_type);
+ return;
+ }
+
+ rotate2d (m_xOutline, m_yOutline, m_nPoints, m_rot);
+ xlat2d (m_xOutline, m_yOutline, m_nPoints, m_cx, m_cy);
+
+ minmax_array (m_xOutline, m_nPoints, m_xmin, m_xmax);
+ minmax_array (m_yOutline, m_nPoints, m_ymin, m_ymax);
+
+ // increase pelem extent by SCALE_PELEM_EXTENT to eliminate chance of
+ // missing actual pelem maximum due to polygonal sampling
+
+ xfact = (m_xmax - m_xmin) * SCALE_PELEM_EXTENT;
+ yfact = (m_ymax - m_ymin) * SCALE_PELEM_EXTENT;
+
+ m_xmin -= xfact;
+ m_ymin -= yfact;
+ m_xmax += xfact;
+ m_ymax += yfact;
+}
+
+
+
+/* NAME
+ * calc_arc Calculate outline of a arc of a circle
+ *
+ * SYNOPSIS
+ * calc_arc (x, y, xcent, ycent, pts, r, start, stop)
+ * double x[], y[]; Array of points
+ * int pts Number of points in array
+ * double xcent, ycent Center of cirlce
+ * double r Radius of circle
+ * double start, stop Beginning & ending angles
+ */
+
+void
+PhantomElement::calcArcPoints (double x[], double y[], const int pts, const double xcent, const double ycent, const double r, const double start, const double stop)
+{
+ if (r <= 0.0)
+ sys_error (ERR_WARNING, "negative or zero radius in calc_arc()");
+
+ double theta = (stop - start) / (pts - 1); // angle incr. between points
+ double c = cos(theta);
+ double s = sin(theta);
+
+ x[0] = r * cos (start) + xcent;
+ y[0] = r * sin (start) + ycent;
+
+ double xp = x[0] - xcent;
+ double yp = y[0] - ycent;
+ for (int i = 1; i < pts; i++) {
+ double xc = c * xp - s * yp;
+ double yc = s * xp + c * yp;
+ x[i] = xc + xcent;
+ y[i] = yc + ycent;
+ xp = xc; yp = yc;
+ }
+}
+
+
+// NAME
+// PhantomElement::calcEllipsePoints Calculate outline of a ellipse
+//
+// SYNOPSIS
+// calcEllipsePoints ()
+//
+
+
+void
+PhantomElement::calcEllipsePoints (double x[], double y[], const int pts, const double u, const double v)
+{
+ calcArcPoints (x, y, m_nPoints, 0.0, 0.0, 1.0, 0.0, TWOPI); // make a unit circle
+ scale2d (x, y, m_nPoints, m_u, m_v); // scale to ellipse
+}
+
+
+/* NAME
+ * circle_pts Calculate number of points to use for circle segment
+ *
+ * SYNOPSIS
+ * n = circle_pts (theta)
+ * int n Number of points to use for arc
+ * double theta Length of arc in radians
+ */
+
+int
+PhantomElement::numCirclePoints (double theta) const
+{
+ if (theta < 0.0 || theta > TWOPI)
+ sys_error(ERR_WARNING, "illegal values sent to circle_pts");
+
+ return static_cast<int> (POINTS_PER_CIRCLE * theta / TWOPI + 1.5);
+}
+
+
+bool
+PhantomElement::clipLineWorldCoords (double& x1, double& y1, double& x2, double &y2) const
+{
+ /* check if ray is outside of pelem extents */
+ double cx1 = x1, cy1 = y1, cx2 = x2, cy2 = y2;
+ if (! clip_rect (cx1, cy1, cx2, cy2, m_rectLimits))
+ return false;
+
+ // convert phantom coordinates to pelem coordinates
+ xform_mtx2 (m_xformPhmToObj, x1, y1);
+ xform_mtx2 (m_xformPhmToObj, x2, y2);
+
+ if (! clipLineNormalizedCoords (x1, y1, x2, y2))
+ return false;
+
+ // convert standard pelem coordinates back to phantom coordinates
+ xform_mtx2 (m_xformObjToPhm, x1, y1);
+ xform_mtx2 (m_xformObjToPhm, x2, y2);
+
+ return true;
+}
+
+
+/* NAME
+ * pelem_clip_line Clip pelem against an arbitrary line
+ *
+ * SYNOPSIS
+ * pelem_clip_line (pelem, x1, y1, x2, y2)
+ * PhantomElement& pelem; Pelem to be clipped
+ * double *x1, *y1, *x2, *y2 Endpoints of line to be clipped
+ *
+ * RETURNS
+ * true if line passes through pelem
+ * (x1, y1, x2, y2 hold coordinates of new line)
+ * false if line do not pass through pelem
+ * (x1, y1, x2, y2 are undefined)
+ */
+
+bool
+PhantomElement::clipLineNormalizedCoords (double& x1, double& y1, double& x2, double& y2) const
+{
+ bool accept = false;
+
+ switch (m_type) {
+ case PELEM_RECTANGLE:
+ double rect[4];
+ rect[0] = -1.0; rect[1] = -1.0;
+ rect[2] = 1.0; rect[3] = 1.0;
+ accept = clip_rect (x1, y1, x2, y2, rect);
+ break;
+ case PELEM_ELLIPSE:
+ accept = clip_circle (x1, y1, x2, y2, 0.0, 0.0, 1.0, 0.0, 0.0);
+ break;
+ case PELEM_TRIANGLE:
+ accept = clip_triangle (x1, y1, x2, y2, 1.0, 1.0, true);
+ break;
+ case PELEM_SEGMENT:
+ accept = clip_segment (x1, y1, x2, y2, m_u, m_v);
+ break;
+ case PELEM_SECTOR:
+ accept = clip_sector (x1, y1, x2, y2, m_u, m_v);
+ break;
+ default:
+ sys_error (ERR_WARNING, "Illegal pelem type %d [pelem_clip_line]", m_type);
+ break;
+ }
+
+ return(accept);
+}
+
+
+// METHOD IDENTIFICATION
+// PhantomElement::isPointInside Check if point is inside pelem
+//
+// SYNOPSIS
+// is_point_inside (pelem, x, y, coord_type)
+// double x, y Point to see if lies in pelem
+// int coord_type Coordinate type (PELEM_COORD or PHM_COORD)
+//
+// RETURNS
+// true if point lies within pelem
+// false if point lies outside of pelem
+
+bool
+PhantomElement::isPointInside (double x, double y, const CoordType coord_type)
+{
+ if (coord_type == PHM_COORD) {
+ xform_mtx2 (m_xformPhmToObj, x, y);
+ } else if (coord_type != PELEM_COORD) {
+ sys_error(ERR_WARNING, "Illegal coordinate type in pelem_is_point_inside");
+ return (false);
+ }
+
+ switch (m_type) {
+ case PELEM_RECTANGLE:
+ if (x > 1. || x < -1. || y > 1. || y < -1.)
+ return (false);
+ else
+ return (true);
+ break;
+ case PELEM_TRIANGLE:
+ if (y < 0. || y > 1. - x || y > 1. + x)
+ return (false);
+ else
+ return (true);
+ break;
+ case PELEM_ELLIPSE:
+ if (x > 1. || x < -1. || y > 1. || y < -1.)
+ return (false);
+ if (x * x + y * y > 1.) // check if inside unit circle
+ return (false);
+ else
+ return (true);
+ break;
+
+ // for clipping segments & sectors, must NOT scale by (1/u, 1/v)
+ // because this destroys information about size of arc component
+
+ case PELEM_SEGMENT:
+ if (x > 1. || x < -1. || y > 0.)
+ return (false); // clip against y > 0
+ x *= m_u; // put back u & v scale
+ y *= m_v;
+ if (x * x + (y-m_v) * (y-m_v) > m_u * m_u + m_v * m_v)
+ return (false); // clip against circle, r = sqrt(@)
+ else
+ return (true);
+ break;
+ case PELEM_SECTOR:
+ if (x > 1. || x < -1. || y > 1.) // extent
+ return (false);
+ if (y > 1. - x || y > 1. + x) // triangle
+ return (false); // clip against triangle
+ x *= m_u; // circle: put back u & v scale
+ y *= m_v;
+ if (x * x + (y-m_v) * (y-m_v) > m_u * m_u + m_v * m_v)
+ return (false); // clip against circle
+ else
+ return (true);
+ break;
+ default:
+ sys_error (ERR_WARNING, "Illegal pelem type in pelem_is_point_inside()");
+ break;
+ }
+
+ return (false);
+}
+
+
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: phm2image.cpp
+** Purpose: Convert phantom objects to rasterized images
+** Programmer: Kevin Rosenberg
+** Date Started: 1984
+**
+** This part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: phm2image.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+#include "ct.h"
+
+
+/* NAME
+ * pic_to_imagefile Make image array from Phantom
+ *
+ * SYNOPSIS
+ * pic_to_imagefile (pic, im, nsample)
+ * Phantom& pic Phantom definitions
+ * ImageFile *im Computed pixel array
+ * int nsample Number of samples along each axis for each pixel
+ * (total samples per pixel = nsample * nsample)
+ */
+
+void
+phm_to_imagefile (const Phantom& phm, ImageFile& im, const int col_start, const int col_count, const int in_nsample, const int trace)
+{
+ int nsample = in_nsample;
+ double dx = phm.xmax() - phm.xmin();
+ double dy = phm.ymax() - phm.ymin();
+ double xcent = phm.xmin() + dx / 2;
+ double ycent = phm.ymin() + dy / 2;
+ double phmlen = (dx > dy ? dx : dy);
+
+ double phmradius = phmlen / 2;
+
+ double xmin = xcent - phmradius;
+ double xmax = xcent + phmradius;
+ double ymin = ycent - phmradius;
+ double ymax = ycent + phmradius;
+
+ im.setAxisExtent (xmin, xmax, ymin, ymax);
+ int nx = im.nx();
+ int ny = im.ny();
+
+/* Each pixel holds the average of the intensity of the cell with (ix,iy)
+ * at the center of the pixel
+ *
+ * Set major increments so that the last cell v[nx-1][ny-1] will start at
+ * (xmax - xinc, ymax - yinc).
+ *
+ * Set minor increments so that sample points are centered in cell
+ */
+
+ if (nx < 2 || ny < 2)
+ return;
+
+ double xinc = (xmax - xmin) / nx;
+ double yinc = (ymax - ymin) / ny;
+ im.setAxisIncrement (xinc, yinc);
+
+ if (nsample < 1)
+ nsample = 1;
+
+ double kxinc = xinc / nsample; /* interval between samples */
+ double kyinc = yinc / nsample;
+ double kxofs = kxinc / 2; /* offset of 1st point */
+ double kyofs = kyinc / 2;
+
+#ifdef HAVE_SGP
+ SGP_ID gid = NULL;
+ if (trace > TRACE_TEXT) {
+ gid = sgp2_init (0, 0, "Phantom to Image");
+ sgp2_window (xmin, ymin, xmax, ymax);
+ sgp2_frame_vpt();
+ }
+#endif
+
+ ImageFileArray v = im.getArray();
+
+ double x_start = xmin + (col_start * xinc);
+ for (PElemConstIterator i = phm.listPElem().begin(); i != phm.listPElem().end(); i++) {
+#ifdef HAVE_SGP
+ if (trace > TRACE_TEXT)
+ sgp2_polyline_abs ((*i)->xOutline(), (*i)->yOutline(), (*i)->nOutlinePoints());
+#endif
+ double x, y, xi, yi;
+ int ix, iy, kx, ky;
+ for (ix = 0, x = x_start; ix < col_count; ix++, x += xinc) {
+ for (iy = 0, y = ymin; iy < ny; iy++, y += yinc) {
+ for (kx = 0, xi = x + kxofs; kx < nsample; kx++, xi += kxinc) {
+ for (ky = 0, yi = y + kyofs; ky < nsample; ky++, yi += kyinc)
+ if ((*i)->isPointInside (xi, yi, PHM_COORD) == TRUE)
+ v[ix][iy] += (*i)->atten();
+ } // for kx
+ } /* for iy */
+ } /* for ix */
+ } /* for pelem */
+
+
+ if (nsample > 0) {
+ double factor = 1.0 / (nsample * nsample);
+
+ for (int ix = 0; ix < col_count; ix++)
+ for (int iy = 0; iy < ny; iy++)
+ v[ix][iy] *= factor;
+ }
+
+#if HAVE_SGP
+ if (trace > TRACE_TEXT)
+ sgp2_close(gid);
+#endif
+
+}
+
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: projections.cpp Projection data classes
+** Programmer: Kevin Rosenberg
+** Date Started: Aug 84
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: projections.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+#include "ct.h"
+
+
+/* NAME
+ * Projections Constructor for projections matrix storage
+ *
+ * SYNOPSIS
+ * proj = projections_create (filename, nView, nDet)
+ * Projections& proj Allocated projections structure & matrix
+ * int nView Number of rotated view
+ * int nDet Number of detectors
+ *
+ */
+
+Projections::Projections (const Scanner& scanner)
+{
+ init (scanner.nView(), scanner.nDet());
+
+ m_phmLen = scanner.phmLen();
+ m_rotInc = scanner.rotInc();
+ m_detInc = scanner.detInc();
+ m_rotStart = 0;
+ m_detStart = -scanner.radius() + (scanner.detInc() / 2);
+ m_phmLen = scanner.phmLen();
+}
+
+
+Projections::Projections (const int nView, const int nDet)
+{
+ init (nView, nDet);
+}
+
+Projections::Projections (void)
+{
+ init (0, 0);
+}
+
+Projections::~Projections (void)
+{
+ deleteProjData();
+}
+
+
+void
+Projections::init (const int nView, const int nDet)
+{
+ m_nView = nView;
+ m_nDet = nDet;
+ m_projData = NULL;
+ newProjData ();
+ m_fd = -1;
+}
+
+
+// NAME
+// newProjData
+
+void
+Projections::newProjData (void)
+{
+ if (m_projData)
+ sys_error(ERR_WARNING, "m_projData != NULL [newProjData]");
+
+ if (m_nView > 0 && m_nDet) {
+ m_projData = new DetectorArray* [m_nView];
+
+ for (int i = 0; i < m_nView; i++)
+ m_projData[i] = new DetectorArray (m_nDet);
+ }
+}
+
+
+/* NAME
+ * projections_free Free memory allocated to projections
+ *
+ * SYNOPSIS
+ * projections_free(proj)
+ * Projections& proj Projectionss to be deallocated
+ */
+
+void
+Projections::deleteProjData (void)
+{
+ if (m_projData != NULL) {
+ for (int i = 0; i < m_nView; i++)
+ delete m_projData[i];
+
+ delete m_projData;
+ m_projData = NULL;
+ }
+}
+
+
+/* NAME
+ * projections_write_header Write data header for projections file
+ *
+ */
+
+bool
+Projections::headerWrite (void)
+{
+ kuint32 _hsize = m_headerSize;
+ kuint32 _nView = m_nView;
+ kuint32 _nDet = m_nDet;
+ kuint32 _geom = m_geometry;
+ kuint32 _remarksize = m_remark.length();
+ kfloat64 _calcTime = m_calcTime;
+ kfloat64 _rotStart = m_rotStart;
+ kfloat64 _rotInc = m_rotInc;
+ kfloat64 _detStart = m_detStart;
+ kfloat64 _detInc = m_detInc;
+ kfloat64 _phmLen = m_phmLen;
+
+ if (m_fd < 0) {
+ sys_error (ERR_SEVERE, "m_fd < 0 [projections_write_header]");
+ 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;
+ }
+
+ return true;
+}
+
+
+/* NAME
+ * projections_read_header Read data header for projections file
+ *
+ */
+
+bool
+Projections::headerRead (void)
+{
+ kuint32 _hsize;
+ kuint32 _nView;
+ kuint32 _nDet;
+ kuint32 _geom;
+ kuint32 _remarksize = 0;
+ kfloat64 _calcTime;
+ kfloat64 _rotStart;
+ kfloat64 _rotInc;
+ kfloat64 _detStart;
+ 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;
+ }
+
+ off_t testPos;
+ if (! read_nint32 (&_hsize, m_fd)) {
+ sys_error (ERR_SEVERE, "Error reading header size , _hsize=%d [projections_read_header]", _hsize);
+ 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;
+ }
+
+ char remarkStorage[_remarksize+1];
+ if (::read (m_fd, &remarkStorage, _remarksize) != _remarksize) {
+ 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) {
+ 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;
+ }
+
+ m_headerSize = _hsize;
+ m_nView = _nView;
+ m_nDet = _nDet;
+ m_geometry = _geom;
+ m_calcTime = _calcTime;
+ m_rotStart = _rotStart;
+ m_rotInc = _rotInc;
+ m_detStart = _detStart;
+ m_detInc = _detInc;
+ m_phmLen = _phmLen;
+
+ 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);
+ return false;
+ }
+
+ if (! headerRead()) {
+ sys_error (ERR_SEVERE, "Error reading projections file %s [open_projections]", filename);
+ return false;
+ }
+
+ deleteProjData ();
+ newProjData ();
+
+ for (int i = 0; i < m_nView; i++) {
+ if (! detarrayRead (*m_projData[i], i))
+ break;
+ }
+
+ close (m_fd);
+ m_fd = -1;
+
+ 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;
+ }
+
+ headerWrite ();
+
+ if (m_projData != NULL) {
+ for (int i = 0; i < m_nView; i++) {
+ if (! detarrayWrite (*m_projData[i], i))
+ break;
+ }
+ }
+
+ close (m_fd);
+ m_fd = -1;
+
+ return true;
+}
+
+/* NAME
+ * detarrayRead Read a Detector Array structure from the disk
+ *
+ * SYNOPSIS
+ * detarrayRead (proj, darray, view_num)
+ * DETARRAY *darray Detector array storage location to be filled
+ * int view_num View number to read
+ */
+
+bool
+Projections::detarrayRead (DetectorArray& darray, const int iview)
+{
+ const int detval_bytes = darray.nDet() * sizeof(kfloat32);
+ const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
+ const int view_bytes = detheader_bytes + detval_bytes;
+ const off_t start_data = m_headerSize + (iview * view_bytes);
+ DetectorValue* detval_ptr = darray.detValues();
+ kfloat64 view_angle;
+ kuint32 nDet;
+
+ if (m_fd < 0) {
+ sys_error (ERR_SEVERE, "Tried to read detector array on closed file [detarrayRead]");
+ return false;
+ }
+
+ 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;
+ }
+ 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;
+ }
+ detval_ptr[i] = detval;
+ }
+
+ return true;
+}
+
+
+/* NAME
+ * detarrayWrite Write detector array data to the disk
+ *
+ * SYNOPSIS
+ * detarrayWrite (darray, view_num)
+ * DETARRAY *darray Detector array data to be written
+ * int view_num View number to write
+ *
+ * DESCRIPTION
+ * This routine writes the detarray data from the disk sequentially to
+ * the file that was opened with open_projections(). Data is written in
+ * binary format.
+ */
+
+bool
+Projections::detarrayWrite (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 */;
+ const int view_bytes = detheader_bytes + detval_bytes;
+ const off_t start_data = m_headerSize + (iview * view_bytes);
+ const DetectorValue* const detval_ptr = darray.detValues();
+ 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]");
+ return false;
+ }
+
+ if (lseek(m_fd, start_data, SEEK_SET) != start_data) {
+ sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]");
+ return false;
+ }
+
+ 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;
+ }
+ }
+
+ return true;
+}
+
+
+/* NAME
+ * prt_projections Print projections data
+ *
+ * SYNOPSIS
+ * prt_projections (proj)
+ * Projections& proj Projection data to be printed
+ */
+
+void
+Projections::printProjectionData (void)
+{
+ printf("Projections Print\n\n");
+ printf("Description: %s\n", m_remark.c_str());
+ printf("nView = %d nDet = %d\n", m_nView, m_nDet);
+ printf("rotStart = %8.4f rotInc = %8.4f\n", m_rotStart, m_rotInc);
+ printf("detStart = %8.4f detInc = %8.4f\n", m_detStart, m_detInc);
+
+ if (m_projData != NULL) {
+ for (int ir = 0; ir < m_nView; ir++) {
+ DetectorValue* detval = m_projData[ir]->detValues();
+ for (int id = 0; id < m_projData[ir]->nDet(); id++)
+ printf("%8.4f ", detval[id]);
+ printf("\n");
+ }
+ }
+}
+
+void
+Projections::printScanInfo (void) const
+{
+ printf ("Number of detectors: %d\n", m_nDet);
+ printf (" Number of views: %d\n", m_nView);
+ printf (" Remark: %s\n", m_remark.c_str());
+ printf (" phmLen: %f\n", m_phmLen);
+ printf (" detStart: %f\n", m_detStart);
+ printf (" detInc: %f\n", m_detInc);
+ printf (" rotStart: %f\n", m_rotStart);
+ printf (" rotInc: %f\n", m_rotInc);
+}
+
+
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: rec_t.c Image Reconstruction Procedures
+** Programmer: Kevin Rosenberg
+** Date Started: Aug 84
+**
+** GLOBAL FUNCTIONS
+** proj_reconst() Reconstruct Image from Projections
+**
+** 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 $
+**
+** 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
+ * proj_reconst Reconstruct Image from Projections
+ *
+ * SYNOPSIS
+ * im = proj_reconst (im, proj, filt_type, filt_param, interp_type)
+ * IMAGE *im Output image
+ * RAYSUM *proj Raysum data if in memory
+ * int filt_type Type of convolution filter to use
+ * double filt_param Filter specific parameter
+ * Currently, used only with Hamming filters
+ * int interp_type Type of interpolation method to use
+ *
+ * ALGORITHM
+ *
+ * Calculate one-dimensional filter in spatial domain
+ * Allocate & clear (zero) the 2d output image array
+ * For each projection view
+ * Convolve raysum array with filter
+ * Backproject raysums and add (summate) to image array
+ * end
+ */
+
+ImageFile&
+proj_reconst (ImageFile& im, Projections& proj, const FilterType filt_type, double filt_param, InterpolationType interp_type, int interp_param, const BackprojType backproj_type, const int trace)
+{
+ int ndet = proj.nDet();
+ 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
+
+#ifdef HAVE_BSPLINE_INTERP
+ int spline_order = 0, zoom_factor = 0;
+ if (interp_type == I_BSPLINE) {
+ zoom_factor = interp_param;
+ spline_order = 3;
+ zoom_factor = 3;
+ n_filtered_proj = (ndet - 1) * (zoom_factor + 1) + 1;
+ }
+#endif
+
+ int n_vec_filter = 2 * ndet - 1;
+ 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);
+
+#if HAVE_SGP
+ double *plot_xaxis = NULL; /* array for plotting */
+ SGP_ID gid;
+ plot_xaxis = new double [n_vec_filter];
+
+ if (trace > TRACE_TEXT) {
+ int i;
+ double f;
+ double filterInc = (filterMax - filterMin) / (n_vec_filter - 1);
+ 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);
+ cio_put_str ("Press any key to continue");
+ cio_kb_getc ();
+ sgp2_close (gid);
+ }
+ if (trace >= TRACE_TEXT) {
+ printf ("nview=%d, ndet=%d, det_start=%.4f, det_inc=%.4f\n", proj.nView(), proj.nDet(), proj.detStart(), proj.detInc());
+ }
+#endif //HAVE_SGP
+
+ Backproject* bj = selectBackprojector (backproj_type, proj, im, interp_type);
+
+ for (int iview = 0; iview < nview; iview++) {
+ if (trace >= TRACE_TEXT)
+ printf ("Reconstructing view %d (last = %d)\n", iview, nview - 1);
+
+ DetectorArray& darray = proj.getDetectorArray (iview);
+ 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);
+
+#ifdef HAVE_SGP
+ if (trace >= TRACE_PLOT) {
+ ezset ("clear.");
+ ezset ("xticks major 5.");
+ ezset ("xlabel ");
+ ezset ("ylabel ");
+ ezset ("xlength .5.");
+ ezset ("box.");
+ ezset ("grid.");
+ ezset ("ufinish yes.");
+ ezplot (vec_proj, plot_xaxis, proj.nDet());
+ ezset ("clear.");
+ ezset ("xticks major 5.");
+ ezset ("xlabel ");
+ ezset ("ylabel ");
+ ezset ("ustart yes.");
+ ezset ("xporigin .5.");
+ ezset ("xlength .5.");
+ ezset ("box");
+ ezset ("grid");
+ gid = ezplot (filtered_proj, plot_xaxis, proj.nDet());
+ }
+#endif //HAVE_SGP
+
+#ifdef HAVE_BSPLINE_INTERP
+ if (interp_type == I_BSPLINE)
+ bspline (proj.nDet(), zoom_factor, spline_order, filtered_proj, filtered_proj);
+
+#ifdef HAVE_SGP
+ if (trace >= TRACE_PLOT && interp_type == I_BSPLINE) {
+ bspline (proj.nDet(), zoom_factor, spline_order, filtered_proj, filtered_proj);
+ ezplot_1d (filtered_proj, n_filtered_proj);
+ }
+#endif
+#endif
+
+ bj->BackprojectView (filtered_proj, darray.viewAngle());
+
+#ifdef HAVE_SGP
+ if (trace >= TRACE_PLOT) {
+ char str[256];
+ printf ("Do you want to exit with current pic (y/n) -- ");
+ fgets(str, sizeof(str), stdin);
+ sgp2_close (sgp2_get_active_win());
+ if (tolower(str[0]) == 'y') {
+ break;
+ }
+ }
+#endif //HAVE_SGP
+ }
+
+ delete bj;
+ delete vec_proj;
+ delete filtered_proj;
+ delete vec_filter;
+
+#ifdef HAVE_SGP
+ delete plot_xaxis;
+#endif
+
+ return (im);
+}
+
--- /dev/null
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: scanner.cpp
+** Purpose: Classes for CT scanner
+** Programmer: Kevin Rosenberg
+** Date Started: 1984
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: scanner.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+#include "ct.h"
+
+
+// NAME
+// DetectorArray Construct a DetectorArray
+
+DetectorArray::DetectorArray (const int nDet)
+{
+ m_nDet = nDet;
+ m_detValues = new DetectorValue [m_nDet];
+}
+
+
+// NAME
+// ~DetectorArray Free memory allocated to a detector array
+
+DetectorArray::~DetectorArray (void)
+{
+ delete [] m_detValues;
+}
+
+
+
+/* NAME
+ * Scanner::Scanner Construct a user specified detector structure
+ *
+ * SYNOPSIS
+ * Scanner (phm, nDet, nView, nSample)
+ * Phantom& phm PHANTOM that we are making detector for
+ * int geomety Geometry of detector
+ * int nDet Number of detector along detector array
+ * int nView Number of rotated views
+ * int nSample Number of rays per detector
+ */
+
+Scanner::Scanner (const Phantom& phm, const ScannerGeometry geometry, int nDet, int nView, int nSample, const double rot_anglen)
+{
+ m_phmLen = phm.maxAxisLength(); // maximal length along an axis
+
+ if (nView < 1)
+ nView = 1;
+ if (nSample < 1)
+ m_nSample = 1;
+ if (nDet < 1)
+ nDet = 1;
+ if ((nDet % 2) == 0)
+ ++nDet; // ensure odd number of detectors
+
+ m_geometry = geometry;
+ m_nDet = nDet;
+ m_nView = nView;
+ m_nSample = nSample;
+ m_detLen = SQRT2 * m_phmLen * ((m_nDet + N_EXTRA_DETECTORS) / static_cast<double>(m_nDet));
+
+ m_rotLen = rot_anglen;
+
+ m_radius = m_detLen / 2;
+ m_detInc = m_detLen / m_nDet;
+ m_rotInc = m_rotLen / m_nView;
+
+ m_initPos.xd1 = m_detLen / 2;
+ m_initPos.yd1 = -m_detLen / 2;
+ m_initPos.xd2 = m_detLen / 2;
+ m_initPos.yd2 = m_detLen / 2;
+ m_initPos.xs1 = -m_detLen / 2;
+ m_initPos.ys1 = -m_detLen / 2;
+ m_initPos.xs2 = -m_detLen / 2;
+ m_initPos.ys2 = m_detLen / 2;
+ m_initPos.angle = 0.0;
+}
+
+Scanner::~Scanner (void)
+{
+}
+
+
+
+/* NAME
+ * raysum_collect Calculate ray sums for a Phantom
+ *
+ * SYNOPSIS
+ * rs = raysum_collect (det, phm, start_view, trace, unit_pulse)
+ * Scanner& det Scanner specifications**
+ * RAYSUM *rs Calculated ray sum matrix
+ * Phantom& phm Phantom we are taking ray sums of
+ * int trace Boolean flag to signal ray sum tracing
+*/
+
+void
+Scanner::collectProjections (Projections& proj, const Phantom& phm, const int start_view, const int trace)
+{
+ GRFMTX_2D rotmtx_initial, temp;
+ GRFMTX_2D rotmtx_incr;
+
+ double start_angle = start_view * proj.rotInc();
+ double xcent = phm.xmin() + (phm.xmax() - phm.xmin()) / 2;
+ double ycent = phm.ymin() + (phm.ymax() - phm.ymin()) / 2;
+
+ double xd1 = xcent + m_initPos.xd1;
+ double yd1 = ycent + m_initPos.yd1;
+ double xd2 = xcent + m_initPos.xd2;
+ double yd2 = ycent + m_initPos.yd2;
+ double xs1 = xcent + m_initPos.xs1;
+ double ys1 = ycent + m_initPos.ys1;
+ double xs2 = xcent + m_initPos.xs2;
+ double ys2 = ycent + m_initPos.ys2;
+
+ m_trace = trace;
+
+#ifdef HAVE_SGP
+ SGP_ID gid;
+ if (m_trace >= TRACE_PHM) {
+ double wsize = 1.42 * m_phmLen / 2; /* sqrt(2) * radius */
+
+ gid = sgp2_init (512, 512, "RayCollect");
+ sgp2_color (C_LTBLUE);
+ sgp2_window (xcent - wsize, ycent - wsize, xcent + wsize, ycent + wsize);
+ sgp2_color (C_BROWN);
+#if RADIUS
+ sgp2_draw_circle (m_phmLen / 2);
+#else
+ sgp2_draw_rect (xcent - m_phmLen / 2, ycent - m_phmLen / 2,
+ xcent + m_phmLen / 2, ycent + m_phmLen / 2);
+#endif
+ sgp2_color (C_BROWN);
+ sgp2_move_abs (0., 0.);
+ sgp2_draw_circle (wsize);
+ // raysum_trace_menu_column = (crt->xsize * crt->asp) / 8 + 3;
+ traceShowParam ("X-Ray Simulator", "%s", RAYSUM_TRACE_ROW_TITLE, 8+C_LTWHITE, " ");
+ traceShowParam ("---------------", "%s", RAYSUM_TRACE_ROW_TITLE2, 8+C_LTWHITE, " ");
+ traceShowParam ("Phantom:", "%s", RAYSUM_TRACE_ROW_PHANT_ID, C_YELLOW, " Herman");
+ traceShowParam ("Chomaticity :", "%s", RAYSUM_TRACE_ROW_CHROMATIC, C_LTGREEN, "Mono");
+ traceShowParam ("Scatter :", "%5.1f", RAYSUM_TRACE_ROW_SCATTER, C_LTGREEN, 0.);
+ traceShowParam ("Photon Uncert:", "%5.1f", RAYSUM_TRACE_ROW_PHOT_STAT, C_LTGREEN, 0.);
+ traceShowParam ("Num Scanners:", "%5d", RAYSUM_TRACE_ROW_NDET, C_LTRED, proj.nDet());
+ traceShowParam ("Num Views :", "%5d", RAYSUM_TRACE_ROW_NVIEW, C_LTRED, proj.nView());
+ traceShowParam ("Samples / Ray:", "%5d", RAYSUM_TRACE_ROW_SAMPLES, C_LTRED, m_nSample);
+
+ sgp2_color (C_LTGREEN);
+ phm.draw();
+
+ initmarker (BDIAMOND, 129);
+ }
+#endif
+
+/* Calculate initial rotation matrix */
+ xlat_mtx2 (rotmtx_initial, -xcent, -ycent);
+ rot_mtx2 (temp, start_angle);
+ mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
+ xlat_mtx2 (temp, xcent, ycent);
+ mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
+
+ xform_mtx2 (rotmtx_initial, xd1, yd1); /* rotate detector endpoints */
+ xform_mtx2 (rotmtx_initial, xd2, yd2); /* to initial view_angle */
+ xform_mtx2 (rotmtx_initial, xs1, ys1);
+ xform_mtx2 (rotmtx_initial, xs2, ys2);
+
+/* Calculate incrementatal rotation matrix */
+ xlat_mtx2 (rotmtx_incr, -xcent, -ycent);
+ rot_mtx2 (temp, proj.rotInc());
+ mult_mtx2 (rotmtx_incr, temp, rotmtx_incr);
+ xlat_mtx2 (temp, xcent, ycent);
+ mult_mtx2 (rotmtx_incr, temp, rotmtx_incr);
+
+ int iview;
+ double viewAngle;
+ for (iview = 0, viewAngle = start_angle; iview < proj.nView(); iview++, viewAngle += proj.rotInc()) {
+ DetectorArray& detArray = proj.getDetectorArray( iview );
+
+#ifdef HAVE_SGP
+ if (m_trace >= TRACE_PHM) {
+ sgp2_move_abs (xd1, yd1);
+ sgp2_line_abs (xd2, yd2);
+ sgp2_move_abs (xs1, ys1);
+ sgp2_line_abs (xs2, ys2);
+ }
+#endif
+ if (m_trace)
+ traceShowParam ("Current View :", "%5d", RAYSUM_TRACE_ROW_CURR_VIEW, C_LTMAGENTA, iview);
+
+ projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2);
+ detArray.setViewAngle (viewAngle);
+
+#ifdef HAVE_SGP
+ if (m_trace >= TRACE_PHM) {
+ // rs_plot (detArray, xd1, yd1, xcent, ycent, theta);
+ sgp2_move_abs (xd1, yd1);
+ sgp2_line_abs (xd2, yd2);
+ sgp2_move_abs (xs1, ys1);
+ sgp2_line_abs (xs2, ys2);
+ }
+#endif
+ xform_mtx2 (rotmtx_incr, xd1, yd1); // rotate detector endpoints
+ xform_mtx2 (rotmtx_incr, xd2, yd2);
+ xform_mtx2 (rotmtx_incr, xs1, ys1);
+ xform_mtx2 (rotmtx_incr, xs2, ys2);
+ } /* for each iview */
+}
+
+
+/* NAME
+ * rayview Calculate raysums for a view at any angle
+ *
+ * SYNOPSIS
+ * rayview (phm, detArray, xd1, nSample, yd1, xd2, yd2, xs1, ys1, xs2, ys2)
+ * Phantom& phm Phantom to scan
+ * DETARRAY *detArray Storage of values for detector array
+ * Scanner& det Scanner parameters
+ * double xd1, yd1, xd2, yd2 Beginning & ending detector positions
+ * double xs1, ys1, xs2, ys2 Beginning & ending source positions
+ *
+ * RAY POSITIONING
+ * For each detector, have there are a variable number of rays traced.
+ * The source of each ray is the center of the source x-ray cell. The
+ * detector positions are equally spaced within the cell
+ *
+ * The increments between rays are calculated so that the cells start
+ * at the beginning of a detector cell and they end on the endpoint
+ * of the cell. Thus, the last cell starts at (xd2-ddx),(yd2-ddy).
+ * The exception to this is if there is only one ray per detector.
+ * In that case, the detector position is the center of the detector cell.
+ */
+
+void
+Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const double xd1, const double yd1, const double xd2, const double yd2, const double xs1, const double ys1, const double xs2, const double ys2)
+{
+ double ddx = (xd2 - xd1) / detArray.nDet(); // change in coords between detectors
+ double ddy = (yd2 - yd1) / detArray.nDet();
+ double sdx = (xs2 - xs1) / detArray.nDet(); // change in coords between source
+ double sdy = (ys2 - ys1) / detArray.nDet();
+
+ double ddx2 = ddx / m_nSample; // Incr. between rays with detector cell
+ double ddy2 = ddy / m_nSample; // Doesn't include detector endpoints
+ double ddx2_ofs = ddx2 / 2; // offset of 1st ray from start of detector cell
+ double ddy2_ofs = ddy2 / 2;
+
+ double xd_maj = xd1 + ddx2_ofs; // Incr. between detector cells
+ double yd_maj = yd1 + ddy2_ofs;
+ double xs_maj = xs1 + (sdx / 2); // put ray source in center of cell
+ double ys_maj = ys1 + (sdy / 2);
+
+ DetectorValue* detval = detArray.detValues();
+
+ if (phm.getComposition() == P_UNIT_PULSE) { // put unit pulse in center of view
+ for (int d = 0; d < detArray.nDet(); d++)
+ if (detArray.nDet() / 2 == d && (d % 2) == 1)
+ detval[d] = 1;
+ else
+ detval[d] = 0;
+ } else {
+ for (int d = 0; d < detArray.nDet(); d++) {
+ double xd = xd_maj;
+ double yd = yd_maj;
+ double xs = xs_maj;
+ double ys = ys_maj;
+ double sum = 0.0;
+ for (int i = 0; i < m_nSample; i++) {
+#ifdef HAVE_SGP
+ if (m_trace >= TRACE_RAYS) {
+ sgp2_move_abs (xs, ys);
+ sgp2_line_abs (xd, yd);
+ }
+#endif
+ sum += projectSingleLine (phm, xd, yd, xs, ys);
+
+ if (m_trace >= TRACE_RAYS)
+ traceShowParam ("Attenuation :", "%5.2f", RAYSUM_TRACE_ROW_ATTEN, C_LTMAGENTA, "sum");
+
+#ifdef HAVE_SGP
+ if (m_trace >= TRACE_RAYS) {
+ sgp2_move_abs (xs, ys);
+ sgp2_line_abs (xd, yd);
+ }
+#endif
+ xd += ddx2;
+ yd += ddy2;
+ }
+
+ detval[d] = sum / m_nSample;
+ xd_maj += ddx;
+ yd_maj += ddy;
+ xs_maj += sdx;
+ ys_maj += sdy;
+ } /* for each detector */
+ } /* if not unit pulse */
+}
+
+
+void
+Scanner::traceShowParam (const char *label, const char *fmt, int row, int color, ...)
+{
+ char s[80];
+ va_list arg;
+
+ va_start(arg, color);
+ // cio_set_cpos (raysum_trace_menu_column, row);
+ snprintf (s, sizeof(s), label, "%s");
+ // cio_set_text_clr (color - 8, 0);
+ cio_put_str (s);
+ vsnprintf (s, sizeof(s), fmt, arg);
+ // cio_set_text_clr (color, 0);
+ cio_put_str (s);
+
+ va_end(arg);
+}
+
+
+/* NAME
+ * projectSingleLine INTERNAL: Calculates raysum along a line for a Phantom
+ *
+ * SYNOPSIS
+ * rsum = phm_ray_attenuation (phm, x1, y1, x2, y2)
+ * double rsum Ray sum of Phantom along given line
+ * Phantom& phm; Phantom from which to calculate raysum
+ * double *x1, *y1, *x2, y2 Endpoints of ray path (in Phantom coords)
+ */
+
+double
+Scanner::projectSingleLine (const Phantom& phm, const double x1, const double y1, const double x2, const double y2)
+{
+ // check ray against each pelem in Phantom
+ double rsum = 0.0;
+ for (PElemConstIterator i = phm.listPElem().begin(); i != phm.listPElem().end(); i++)
+ rsum += projectLineAgainstPElem (**i, x1, y1, x2, y2);
+
+ return (rsum);
+}
+
+
+/* NAME
+ * pelem_ray_attenuation Calculate raysum of an pelem along one line
+ *
+ * SYNOPSIS
+ * rsum = pelem_ray_attenuation (pelem, x1, y1, x2, y2)
+ * double rsum Computed raysum
+ * PhantomElement& pelem Pelem to scan
+ * double x1, y1, x2, y2 Endpoints of raysum line
+ */
+
+double
+Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double y1, double x2, double y2)
+{
+ if (! pelem.clipLineWorldCoords (x1, y1, x2, y2)) {
+ if (m_trace == TRACE_CLIPPING)
+ cio_tone (1000., 0.05);
+ return (0.0);
+ }
+
+#ifdef HAVE_SGP
+ if (m_trace == TRACE_CLIPPING) {
+ sgp2_move_abs (x1, y1);
+ sgp2_line_abs (x2, y2);
+ cio_tone (8000., 0.05);
+ sgp2_move_abs (x1, y1);
+ sgp2_line_abs (x2, y2);
+ }
+#endif
+
+ double len = lineLength (x1, y1, x2, y2);
+ return (len * pelem.atten());
+}
+
--- /dev/null
+noinst_LIBRARIES = libctsupport.a
+INCLUDES=@my_includes@
+libk_a_SOURCES= filefuncs.cpp strfuncs.cpp syserror.cpp timedate.cpp byteorder.cpp audio.cpp consoleio.cpp normangle.cpp simpson.cpp xform.cpp clip.cpp
+EXTRA_DIST=Makefile.nt
+
+
+
--- /dev/null
+# Makefile for libk\r
+\r
+!include <ntwin32.mak>\r
+\r
+CC=cl\r
+LD=link\r
+CFLAGS=-O -nologo -I..\include -DENDIAN_HIGH=0 -DENDIAN_LOW=1\r
+LDFLAGS=\r
+O=.obj\r
+\r
+# variables\r
+OBJ1 = kbasename$(O) allocnum$(O) fexist$(O) iclip$(O) sysalloc$(O) syserror$(O) sysfopen$(O) sysfree$(O) s_head$(O) s_lower$(O) s_rmtail$(O) s_save$(O) timedate$(O) netorder$(O)\r
+\r
+all: libk.lib\r
+\r
+.obj: .c\r
+ $(CC) -c $(cvarsdll) $(CFLAGS) $*.c\r
+\r
+\r
+libk.lib: $(OBJ1)\r
+ echo something to del > libk.lib\r
+ del libk.lib\r
+ lib /out:libk.lib $(OBJ1)\r
+\r
+\r
+clean:\r
+ echo dummy > a.obj\r
+ echo dummy > a.exe\r
+ echo dummy > a.lib\r
+ del *.obj\r
+ del *.exe\r
+ del *.lib\r
--- /dev/null
+/*****************************************************************************
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: audio.cpp,v 1.1 2000/06/19 02:58:08 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+/* NAME
+ * beep sound a beep to user
+ *
+ * SYNOPSIS
+ * beep()
+ */
+
+#include "cio.h"
+
+void cio_beep (void)
+{
+ cio_tone (2000.0, 0.3);
+}
+
+/* NAME
+ * tone play a frequency sound for some duration
+ *
+ * SYNOPSIS
+ * tone (freq, length)
+ * double freq frequency to play in Hertz
+ * double length duration to play note in seconds
+ */
+
+#include <stdio.h>
+#include "cio.h"
+
+void
+cio_tone (double freq, double length)
+{
+#if 1
+ fprintf(stdout, "\007");
+#else
+ cio_spkr_freq (freq); /* Set frequency of tone */
+ cio_spkr_on (); /* Turn on speaker */
+ pause (length); /* Pause for length seconds */
+ cio_spkr_off (); /* Turn off speaker */
+#endif
+}
--- /dev/null
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#if HAVE_UNISTD_H
+#include <unistd.h>
+#endif
+
+#include "kstddef.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
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+** Name: clip.c
+** Purpose: Line clipping routines
+** Programmer: Kevin Rosenberg
+** Date Started: 1984
+**
+** OVERVIEW
+** Routines to clip lines against objects
+** All routines get the endpoints of the line, and
+** the SNARK size of the object (u,v)
+**
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: clip.cpp,v 1.1 2000/06/19 02:58:08 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+#include "kstddef.h"
+#include "kmath.h"
+
+
+/* NAME
+ * clip_segment Clip against a segment of a circle
+ *
+ * SYNOPSIS
+ * clip_segment (x1, y1, x2, y2, u, v)
+ * double& x1,*y1,*x2,*y2 Endpoints of line
+ * double u,v Dimensions of segment
+ */
+
+int
+clip_segment (double& x1, double& y1, double& x2, double& y2, const double u, const double v)
+{
+ double xc1 = x1 * u;
+ double yc1 = y1 * v;
+ double xc2 = x2 * u;
+ double yc2 = y2 * v;
+
+ if (yc1 > 0. && yc2 > 0.) // reject lines above y-axis
+ return false;
+
+ double radius = sqrt (u * u + v * v);
+
+ if (clip_circle (xc1, yc1, xc2, yc2, 0.0, v, radius, 0.0, 0.0) == false)
+ return false;
+
+ if (yc1 > 0. && yc2 > 0.) // trivial reject above y-axis
+ return false;
+
+ // clip above x-axis
+ if (yc1 > 0.) {
+ xc1 = xc1 + (xc2-xc1)*(0.0-yc1)/(yc2-yc1);
+ yc1 = 0.0;
+ } else if (yc2 > 0.) {
+ xc2 = xc1 + (xc2-xc1)*(0.0-yc1)/(yc2-yc1);
+ yc2 = 0.0;
+ }
+
+ x1 = xc1 / u;
+ y1 = yc1 / v;
+ x2 = xc2 / u;
+ y2 = yc2 / v;
+
+ return true;
+}
+
+
+/* NAME
+ * clip_sector Clip a line against a sector of a circle
+ *
+ * SYNOPSIS
+ * clip_sector (x1, y1, x2, y2, u, v)
+ * double& x1,*y1,*x2,*y2 Endpoints of line
+ * double u,v Size of sector
+ */
+
+int
+clip_sector (double& x1, double& y1, double& x2, double& y2, const double u, const double v)
+{
+ double xc1 = x1 * u;
+ double yc1 = y1 * v;
+ double xc2 = x2 * u;
+ double yc2 = y2 * v;
+
+ double radius = sqrt (u * u + v * v);
+
+ if (clip_circle (xc1, yc1, xc2, yc2, 0.0, v, radius, 0.0, 0.0) == false)
+ return false;
+
+ if (clip_triangle (xc1, yc1, xc2, yc2, u, v, false) == false)
+ return false;
+
+ x1 = xc1 / u;
+ y1 = yc1 / v;
+ x2 = xc2 / u;
+ y2 = yc2 / v;
+
+ return true;
+}
+
+
+/* NAME
+ * clip_circle Clip a line against a circle
+ *
+ * SYNOPSIS
+ * clip_circle (x1,y1,x2,y2,cx,cy,radius,t1,t2)
+ * double& x1,*y1,*x2,*y2 Endpoints of line to be clipped
+ * double cx,cy Center of circle
+ * double radius Radius of circle
+ * double t1,t2 Starting & stopping angles of clipping
+ */
+
+int
+clip_circle (double& x1, double& y1, double& x2, double& y2, const double cx, const double cy, const double radius, double t1, double t2)
+{
+ double xc1 = x1;
+ double yc1 = y1;
+ double xc2 = x2;
+ double yc2 = y2;
+ double ccx = cx;
+ double ccy = cy;
+
+ double xtrans = -xc1; // move (x1,y1) to origin
+ double ytrans = -yc1;
+
+ xc1 += xtrans;
+ yc1 += ytrans;
+ xc2 += xtrans;
+ yc2 += ytrans;
+ ccx += xtrans;
+ ccy += ytrans;
+
+ double theta = -atan2 (yc2, xc2); // rotate line to lie on x-axis
+ GRFMTX_2D rotmtx;
+ rot_mtx2 (rotmtx, theta); // xc1,yc1 is at origin, no need to rot
+ xform_mtx2 (rotmtx, xc2, yc2);
+ xform_mtx2 (rotmtx, ccx, ccy);
+ t1 += theta; // rotate start and stop angles
+ t2 += theta;
+ t1 = norm_ang (t1);
+ t2 = norm_ang (t2);
+
+ if (xc2 < -D_EPSILON || fabs(yc2) > F_EPSILON) {
+ sys_error (ERR_SEVERE, "Internal error in clip_circle\n x1=%6.2f, y1=%6.2f, x2=%6.2f, y2=%6.2f, xc2=%6.2f, yc2=%6.2f, theta=%6.2f", x1, y1, x2, y2, xc2, yc2, theta);
+ return false;
+ }
+
+ if (fabs(ccy) > radius) /* check if can reject */
+ return false;
+
+ double temp = sqrt (radius * radius - ccy * ccy);
+ double xcmin = ccx - temp;
+ double xcmax = ccx + temp;
+
+ if (fabs(t2 - t1) < D_EPSILON) {
+ if (xc1 < xcmin)
+ xc1 = xcmin;
+ if (xc2 > xcmax)
+ xc2 = xcmax;
+ } else if (t1 < t2) {
+ if (t1 < PI && t2 > PI)
+ if (xc1 < xcmin)
+ xc1 = xcmin;
+ } else if (t1 > t2) {
+ if (t1 < PI)
+ if (xc1 < xcmin)
+ xc1 = xcmin;
+ if (xc2 > xcmax)
+ xc2 = xcmax;
+ }
+
+ rot_mtx2 (rotmtx, -theta);
+ xform_mtx2 (rotmtx, xc1, yc1);
+ xform_mtx2 (rotmtx, xc2, yc2);
+ xc1 += -xtrans;
+ yc1 += -ytrans;
+ xc2 += -xtrans;
+ yc2 += -ytrans;
+
+ x1 = xc1;
+ y1 = yc1;
+ x2 = xc2;
+ y2 = yc2;
+
+ return true;
+}
+
+
+/* NAME
+ * clip_triangle Clip a line against a triangle
+ *
+ * SYNOPSIS
+ * clip_triangle (x1, y1, x2, y2, u, v, clip_xaxis)
+ * double& x1, *y1, *x2, *y2 Endpoints of line
+ * double u, v Size of 1/2 base len & height
+ * int clip_xaxis Boolean flag whether to clip against x axis
+ * (Use true for all triangles)
+ * (false if used internally by sector clipping routine)
+ *
+ * DESCRIPTION
+ * x
+ * /|\ Note that vertices of triangle are
+ * / | \ (-u, 0)
+ * / | \ (u, 0)
+ * / | \ (0, v)
+ * / | v \
+ * / | \
+ * +------+------+
+ * (0,0) u
+ *
+ * NOTES
+ * 1) Inside of this routine, values of (u,v) are assumed to be (1,1)
+ *
+ * 2) Derivation of clipping equations:
+ * Using parametric equations for the line
+ * xv = x1 + t * (x2 - x1)
+ * yv = y1 + t * (y2 - y1)
+ * so,
+ * t = (xv - x1) / (x2 - x1)
+ * yv = y1 + (xv - x1) * (y2 - y1) / (x2 - x1)
+ * yv = y1 + (xv - x1) * dy / dx
+ *
+ * Now, find the intersections with the following clipping boundries:
+ * yv = v - (v/u) * xv (yv = mx + b)
+ * yv = v + (v/u) * xv (m = v/u, b = v);
+ */
+
+static int tcode (const double x, const double y, const double m, const double b, const int clip_xaxis);
+
+int
+clip_triangle (double& x1, double& y1, double& x2, double& y2, const double u, const double v, const int clip_xaxis)
+{
+ double m = v / u; // slope of triangle lines
+ double b = v; // y-intercept of triangle lines
+
+ int c1 = tcode (x1, y1, m, b, clip_xaxis);
+ int c2 = tcode (x2, y2, m, b, clip_xaxis);
+
+#ifdef DEBUG
+ crt_set_cpos (1,1);
+ printf ("x1:%6.2f y1:%6.2f code1:%2d x2:%6.2f y2:%6.2f code2:%2d",
+ x1, y1, c1, x2, y2, c2);
+#endif
+ while ( c1 || c2 ) {
+ if ( c1 & c2 ) {
+ return false; // trivial reject
+ }
+ int c = c1;
+ if (c1 == 0)
+ c = c2;
+
+ double x = 0, y = 0;
+ if (c & 1) { // below
+ x = x1 + (x2-x1)*(0.0-y1)/(y2-y1);
+ y = 0.0;
+ } else if (c & 2) { // right
+ double dx, dy;
+ dx = x2 - x1;
+ dy = y2 - y1;
+ if (fabs(dx) > D_EPSILON)
+ x = (-y1 + b + x1 * dy/dx) / (m + dy/dx);
+ else
+ x = x1;
+ y = -m * x + b;
+ } else if (c & 4) { /* left */
+ double dx, dy;
+ dx = x2 - x1;
+ dy = y2 - y1;
+ if (fabs(dx) > D_EPSILON) {
+ x = (y1 - b - x1 * dy/dx);
+ x /= (m - dy/dx);
+ } else
+ x = x1;
+ y = m * x + b;
+ }
+
+ if (c == c1) {
+ x1=x; y1=y; c1=tcode (x1,y1,m,b,clip_xaxis);
+ } else {
+ x2=x; y2=y; c2=tcode (x2,y2,m,b,clip_xaxis);
+ }
+#ifdef DEBUG
+ crt_set_cpos (1,1);
+ printf ("x1:%6.2f y1:%6.2f code1:%2d x2:%6.2f y2:%6.2f code2:%2d", x1, y1, c1, x2, y2, c2);
+#endif
+ }
+
+ return true; /* we have clipped the line, and it is good */
+}
+
+
+/* compute region code */
+static int
+tcode (const double x, const double y, const double m, const double b, const int clip_xaxis)
+{
+ int c = 0;
+
+ if (clip_xaxis && y < 0.) // below triange
+ c = 1;
+
+ if (y > -m * x + b + D_EPSILON) // right of triangle
+ c += 2;
+ if (y > m * x + b + D_EPSILON) // left of triangle
+ c += 4;
+
+ return (c);
+}
+
+
+/* NAME
+ * clip_rect Clip a line against a rectangle
+ *
+ * SYNOPSIS
+ * clip_rect (x1, y1, x2, y2, rect)
+ * double& x1, *y1, *x2, *y2 Endpoints of line
+ * double rect[4] Rectangle to clip against
+ * ordered xmin, ymin, xmax, ymax
+ */
+
+static int rectcode (double x, double y, const double rect[4]);
+
+int
+clip_rect (double& x1, double& y1, double& x2, double& y2, const double rect[4])
+{
+ double x = 0, y = 0;
+
+ int c1 = rectcode (x1, y1, rect);
+ int c2 = rectcode (x2, y2, rect);
+
+ while (c1 || c2) {
+ if (c1 & c2)
+ return false; // trivial reject
+ int c = c1;
+ if (c1 == 0)
+ c = c2;
+ if (c & 1) { // left
+ y = y1 + (y2-y1)*(rect[0]-x1)/(x2-x1);
+ x = rect[0];
+ } else if (c & 2) { // right
+ y = y1 + (y2-y1)*(rect[2]-x1)/(x2-x1);
+ x = rect[2];
+ } else if (c & 4) { // bottom
+ x = x1 + (x2-x1)*(rect[1]-y1)/(y2-y1);
+ y = rect[1];
+ } else if (c & 8) { // top
+ x = x1 + (x2-x1)*(rect[3]-y1)/(y2-y1);
+ y = rect[3];
+ }
+
+ if (c == c1) {
+ x1=x; y1=y; c1=rectcode(x1,y1,rect);
+ } else {
+ x2=x; y2=y; c2=rectcode(x2,y2,rect);
+ }
+ }
+ return true; // we have clipped the line, and it is good
+}
+
+
+/* NAME
+ * rectcode INTERNAL routine to return position of
+ * point relative to a rectangle
+ *
+ * SYNOPSIS
+ * c = rectcode (x, y, rect)
+ * int c Position of point relative to the window
+ * double x, y Point to test against window
+ * double rect[4] Coordinates of rectangle extent
+ * Ordered [xmin, ymin, xmax, ymax]
+ */
+
+static int
+rectcode (double x, double y, const double rect[4])
+{
+ int c = 0;
+
+ if (x < rect[0])
+ c = 1;
+ else if (x > rect[2])
+ c = 2;
+ if (y < rect[1])
+ c += 4;
+ else if (y > rect[3])
+ c += 8;
+ return (c);
+}
--- /dev/null
+/*****************************************************************************
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: consoleio.cpp,v 1.1 2000/06/19 02:58:08 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+#include <stdio.h>
+#include <string.h>
+#include "cio.h"
+#include "kstddef.h"
+#include "cio.h"
+
+
+
+/* NAME
+ * cio_put_c Put a character on screen
+ *
+ * SYNOPSIS
+ * cio_put_c (c)
+ * char c Character to write
+ *
+ * NOTES
+ * Color of character is determined by the global variable, crtv.text_attr.
+ *
+ * SIDE EFFECTS
+ * Cursor is advanced by one. If necessary, the cursor will wrap around
+ * and maybe the screen will scroll
+ */
+
+void
+cio_put_c (int c)
+{
+ fputc(c, stdout);
+}
+
+
+
+/* NAME
+ * cio_put_cc Put a char on screen count times
+ *
+ * SYNOPSIS
+ * cio_put_cc (c, count)
+ * char c Character to write
+ * int count Number of characters to write
+ */
+
+void
+cio_put_cc (int c, int count)
+{
+ for (int i = 0; i < count; i++)
+ cio_put_c (c);
+}
+
+
+void
+cio_put_str (const char *str)
+{
+ fputs (str, stdout);
+}
+
+
+
+/* NAME
+ * kb_getc Get a character from the keyboard
+ *
+ * SYNOPSIS
+ * key = kb_getc()
+ *
+ * DESCRIPTION
+ * 1. This routine returns an EXTENTED ASCII code,
+ * the extended codes have a low byte of 0 and a distinctive
+ * high byte, such as 0x2D00 and 0x3200
+ * 2. This routine waits until a key has been typed
+ * 2. The keystroke will not be echoed.
+ */
+
+unsigned int cio_kb_getc(void)
+{
+ return fgetc(stdin);
+}
+
+void
+cio_kb_ungetc (unsigned int c)
+{
+ ungetc(c, stdin);
+}
+
+/* NAME
+ * kb_gets Get a string from the keyboard
+ *
+ * SYNOPSIS
+ * str = kb_gets (str, maxlen)
+ * char *str Space to store input string
+ * int maxlen Maximum number of characters to read
+ * (Not including EOS)
+ * NOTES
+ * Backspace - erases character to the right
+ * Escape - erases to beginning of line
+ * Return - ends string (no not cause a linefeed)
+ */
+
+char *
+cio_kb_gets (char *str, int maxlen)
+{
+ return fgets(str, maxlen, stdin);
+}
+
+/* NAME
+ * kb_waitc Wait for a character from the keyboard
+ *
+ * SYNOPSIS
+ * key = kb_waitc (astr, estr, beep)
+ * int key Keystroke entered
+ * char *astr String of valid ascii characters
+ * bool beep If TRUE, beep when user hits invalid key
+ *
+ */
+
+
+unsigned int
+cio_kb_waitc (const char *astr, int beep_on_error)
+{
+ unsigned int c;
+ do {
+ c = cio_kb_getc ();
+ if (strchr (astr, c) != NULL)
+ break;
+ if (beep_on_error)
+ cio_beep();
+ } while (1);
+
+ return (c);
+}
+
+
--- /dev/null
+/*****************************************************************************
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: filefuncs.cpp,v 1.1 2000/06/19 02:58:08 kevin Exp $
+**
+** Revision 1.1.1.1 2000/04/28 13:02:44 kevin
+** Initial CVS import for first public release
+**
+** 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 <string.h>
+#include <stdio.h>
+#include <ctype.h>
+#include <iostream>
+#include "kstddef.h"
+#include "cio.h"
+
+
+
+const char*
+fileBasename (const char* const filename)
+{
+ const char* p = strrchr (filename, '/');
+ return (p ? p + 1 : filename);
+}
+
+
+/* NAME
+ * file_exists Checks if a specified disk fie exists
+ *
+ * SYNOPSIS
+ * exist = file_exists (fname)
+ * bool exist TRUE if specified file exists
+ */
+
+bool
+file_exists (const char *fname)
+{
+ FILE *fp;
+ bool exist;
+
+ if (strlen(fname) == 0)
+ exist = false;
+ else if ((fp = fopen(fname, "r")) == NULL)
+ exist = false;
+ else {
+ fclose (fp);
+ exist = true;
+ }
+
+ return (exist);
+}
+
+/*-----------------------------------------------------------------------------
+ *
+ * FUNCTION IDENTIFICATION
+ *
+ * Name: sys_fopen Open a file for user
+ * Date: 12-17-84
+ * Programmer: Kevin Rosenberg
+ *
+ * SYNOPSIS
+ * fp = sys_fopen (filename, mode, progname)
+ * FILE *fp Standard pointer to 'C' file
+ * char *filename Name of file to open
+ * If user enters a new name, it goes here
+ * char *mode Mode to open file (std. 'C')
+ * char *progname Name of program calling this routine
+ *
+ * DESCRIPTION
+ * This routine opens a file using the standard C fopen() routine. If
+ * the file is not found, the user is given to option to:
+ *
+ * 1 - Retry opening file with same name
+ * 2 - Enter new file name
+ * 3 - Abort and return to DOS
+ *
+ * CAUTIONS
+ * If the the requested file is not found, the name of the file given
+ * entered at keyboard will be returned in filename. So, make sure there
+ * is room for a maximum length filename (MAXFULLNAME)
+ *
+ *---------------------------------------------------------------------------*/
+
+FILE *
+sys_fopen (const char *filename, const char *mode, const char *progname)
+{
+ FILE *fp;
+ char fname[256]; /* name used for call to fopen() */
+ char c; /* keyboard response */
+
+ strncpy (fname, filename, sizeof(fname));
+
+ do {
+ if ((fp = fopen (fname, mode)) == NULL) {
+ cerr << endl;
+ cerr << "Can't open file " << fname << " [" << progname << "]" << endl;
+ cerr << "Enter: <R> - Retry | <N> - New name | <A> - Abort program --> ";
+ c = cio_kb_waitc ("RrNnAa", TRUE);
+ c = tolower(c);
+ cerr << c << endl;
+
+ if (c == 'r') // Retry -- Nothing to do here
+ ;
+ else if (c == 'a') // Abort -- Exit to OS
+ exit (1);
+ else if (c == 'n') { // New name - get from console
+ cerr << "Enter new file name -- ";
+ fgets (fname, sizeof(fname), stdin);
+ str_wrm_tail (fname);
+ cerr << endl;
+ }
+ }
+ } while (fp == NULL);
+
+ return (fp);
+}
--- /dev/null
+/*****************************************************************************
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: normangle.cpp,v 1.1 2000/06/19 02:58:08 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+#include "kmath.h"
+
+
+/* NAME
+ * norm_ang Normalize angle to 0 to 2 * PI range
+ *
+ * SYNOPSIS
+ * t = norm_ang (theta)
+ * double t Normalized angle
+ * double theta Input angle
+ */
+
+double
+norm_ang (double theta)
+{
+ while (theta < 0.)
+ theta += TWOPI;
+ while (theta >= TWOPI)
+ theta -= TWOPI;
+
+ return (theta);
+}
--- /dev/null
+/*****************************************************************************
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: simpson.cpp,v 1.1 2000/06/19 02:58:08 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+#include "kmath.h"
+
+
+/* NAME
+ * integrateSimpson Integrate array of data by Simpson's rule
+ *
+ * SYNOPSIS
+ * double integrateSimpson (xmin, xmax, y, np)
+ * double xmin, xmax Extent of integration
+ * double y[] Function values to be integrated
+ * int np number of data points
+ * (must be an odd number and at least 3)
+ *
+ * RETURNS
+ * integrand of function
+ */
+
+double
+integrateSimpson (const double xmin, const double xmax, const double *y, const int np)
+{
+ if (np < 2)
+ return (0.);
+ else if (np == 2)
+ return ((xmax - xmin) * (y[0] + y[1]) / 2);
+
+ double area = 0;
+ int nDiv = (np - 1) / 2; // number of divisions
+ double width = (xmax - xmin) / (double) (np - 1); // width of cells
+
+ for (int i = 1; i <= nDiv; i++) {
+ int xr = 2 * i;
+ int xl = xr - 2; // 2 * (i - 1) == 2 * i - 2 == xr - 2
+ int xm = xr - 1; // (xl+xr)/2 == (xr+xr-2)/2 == (2*xr-2)/2 = xr-1
+
+ area += (width / 3.0) * (y[xl] + 4.0 * y[xm] + y[xr]);
+ }
+
+ if ((np & 1) == 0) /* do last trapazoid */
+ area += width * (y[np-2] + y[np-1]) / 2;
+
+ return (area);
+}
--- /dev/null
+/*****************************************************************************
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: strfuncs.cpp,v 1.1 2000/06/19 02:58:08 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+#include <ctype.h>
+#include "kstddef.h"
+
+
+/* NAME
+ * str_skip_head Skip leading characters of string
+ *
+ * SYNOPSIS
+ * shortened = str_skip_head (str, charlist)
+ * OUT shortened Start of shortened string
+ * IN char *str String to have beginning skipped
+ * IN char *charlist List of characters to skip over
+ *
+ * NOTES
+ * This routine returns the position in a string (str) of the
+ * first character that is not in an specified string of characters
+ * (charlist).
+ */
+
+
+const char*
+str_skip_head (const char* str, const char* const charlist)
+{
+ const char* p = str;
+
+ while (*p && (strchr (charlist, *p) != NULL))
+ p++;
+
+ return (p);
+}
+
+char*
+str_skip_head (char* str, char* charlist)
+{
+ char* p = str;
+
+ while (*p && (strchr (charlist, *p) != NULL))
+ p++;
+
+ return (p);
+}
+
+
+/* NAME
+ * str_lower Convert a string to lower case
+ *
+ * SYNOPSIS
+ * str = str_lower (str)
+ * char *str String to be converted
+ */
+
+char *
+str_lower (char *s)
+{
+ char *p = s;
+
+ while (*p) { /* while (*p != EOS) */
+ *p = tolower(*p);
+ ++p;
+ }
+
+ return (s);
+}
+
+/* NAME
+ * str_rm_tail Remove characters from end of string
+ *
+ * SYNOPSIS
+ * str = str_rm_tail (str, charlist)
+ * char *str String to have end removed
+k * char *charlist List of characters to remove from string
+ *
+ */
+
+
+char *
+str_rm_tail (char *str, const char* const charlist)
+{
+ int i;
+
+ for (i = strlen(str) - 1; i >= 0; i--)
+ if (strchr (charlist, str[i]) != NULL)
+ str[i] = EOS;
+ else
+ break; /* found non-specified char, all done */
+
+ return (str);
+}
+
+/* NAME
+ * str_wrm_tail Remove white space from end of string
+ *
+ * SYNOPSIS
+ * str = str_wrm_tail (str)
+ * char *str String to have white space removed
+ *
+ */
+
+char *
+str_wrm_tail (char *str)
+{
+ return (str_rm_tail(str, "\b\t\n\r"));
+}
+
+/* NAME
+ * str_upper Convert a string to upper case
+ *
+ * SYNOPSIS
+ * str = str_upper (str)
+ * char *str String to be converted
+ */
+
+char *
+str_upper (char *s)
+{
+ char *p = s;
+
+ while (*p) { /* while (*s != EOS) */
+ *p = toupper(*p);
+ p++;
+ }
+
+ return (s);
+}
+
+
+#ifdef TEST
+int
+main (void)
+{
+ string str, clist;
+ char *skip;
+
+ printf ("Test program for str_skip_head\n");
+ printf ("\n");
+ printf ("Enter string that will have its head skipped -- ");
+ gets (str);
+ printf ("Enter list of characters to be skipped -- ");
+ gets (clist);
+ printf ("\n");
+
+ skip = str_skip_head (str, clist);
+
+ printf ("Shortened string = '%s'\n", skip);
+}
+#endif
+
+
--- /dev/null
+/*****************************************************************************
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: syserror.cpp,v 1.1 2000/06/19 02:58:08 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+#include <iostream>
+#include <stdarg.h>
+#include <ctype.h>
+#include "kstddef.h"
+#include "cio.h"
+
+
+/* NAME
+ * sys_error System error handler
+ *
+ * SYNOPSIS
+ * sys_error (severity, msg, args . . .)
+ * int severity Severity of error
+ * char *msg Error message
+ * args Argument list, direct transfer to printf stack
+ * Can take 24 byte transfer
+ */
+
+static int errorlevel = ERR_WARNING; /* Set error reporting level */
+
+
+void sys_error (int severity, const char *msg, ...)
+{
+ va_list arg;
+
+ va_start(arg, msg);
+
+ sys_verror (severity, msg, arg);
+
+ va_end(arg);
+}
+
+static int nErrorCount = 0;
+const static int MAX_ERROR_COUNT = 20;
+
+
+void sys_verror (int severity, const char *msg, va_list arg)
+{
+ if (severity < errorlevel)
+ return; /* ignore error if less than max level */
+
+ nErrorCount++;
+ if (severity != ERR_FATAL) {
+ if (nErrorCount > MAX_ERROR_COUNT)
+ return;
+ else if (nErrorCount == MAX_ERROR_COUNT) {
+ cout << "*****************************************************************" << endl;
+ cout << "*** M A X I M U M E R R O R C O U N T R E A C H E D ***" << endl;
+ cout << "*** ***" << endl;
+ cout << "*** No further errors will be reported ***" << endl;
+ cout << "*****************************************************************" << endl;
+ return;
+ }
+ }
+
+ switch (severity) {
+ case ERR_FATAL:
+ cout << "FATAL ERROR: ";
+ break;
+ case ERR_SEVERE:
+ cout << "SEVERE ERROR: ";
+ break;
+ case ERR_WARNING:
+ cout << "WARNING ERROR: ";
+ break;
+ default:
+ sys_error (ERR_FATAL, "illegal error code #%d [sys_error]", severity);
+ }
+
+ vfprintf (stdout, msg, arg);
+
+ cout << "\n";
+
+ if (severity == ERR_FATAL)
+ exit(1);
+
+#if INTERACTIVE_ERROR_DISPLAY
+ cout << "A - Abort C - Continue W - Turn off warnings? ";
+ // fflush(stderr);
+ do
+ {
+ int c = cio_kb_waitc("AaBbCcWw", TRUE); /* get code from keyboard */
+ c = tolower (c);
+ fputc (c, stderr);
+ fputc (NEWLINE, stderr);
+
+ if (c == 'a')
+ exit (1);
+ else if (c == 'c')
+ return;
+ else if (c == 'w')
+ {
+ sys_error_level (ERR_SEVERE); /* report severe & fatal errors */
+ break;
+ }
+ } while (TRUE);
+#endif
+}
+
+
+/* NAME
+ * sys_error_level Set error reporting level
+ *
+ * SYNOPSIS
+ * sys_error_level (severity)
+ * int severity Report all error as serious as severity and beyond
+ *
+ * DESCRIPTION
+ * Causes the system to ignore all error below the level of severity
+ * For example, if severity == ERR_SEVERE, then report severe & fatal
+ * error and ignore warnings
+ */
+
+void
+sys_error_level (int severity)
+{
+ if (severity == ERR_FATAL ||
+ severity == ERR_SEVERE ||
+ severity == ERR_WARNING)
+ errorlevel = severity;
+ else
+ errorlevel = ERR_WARNING;
+}
+
--- /dev/null
+/*****************************************************************************
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: timedate.cpp,v 1.1 2000/06/19 02:58:08 kevin Exp $
+** $Log: timedate.cpp,v $
+** Revision 1.1 2000/06/19 02:58:08 kevin
+** *** empty log message ***
+**
+** Revision 1.1 2000/06/13 16:20:31 kevin
+** finished c++ conversions
+**
+** Revision 1.5 2000/06/05 01:33:11 kevin
+** *** empty log message ***
+**
+** Revision 1.4 2000/05/11 01:05:51 kevin
+** Changed sprintf to snprintf, changed index to strchr
+**
+** Revision 1.3 2000/04/28 18:00:55 kevin
+** remove unused files
+**
+** Revision 1.2 2000/04/28 17:38:26 kevin
+** Removed unused files
+**
+** Revision 1.1.1.1 2000/04/28 13:02:44 kevin
+** Initial CVS import for first public release
+**
+**
+**
+** 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
+******************************************************************************/
+
+/*----------------------------------------------------------------------*/
+/* ROUTINES THAT MANAGE TIME */
+/*----------------------------------------------------------------------*/
+
+#include <stdio.h>
+#include "kstddef.h"
+#include <time.h>
+
+/* NAME
+ * td_get_time Return time of day
+ *
+ * SYNOPSIS
+ * t = td_get_time (t)
+ * TIME *t; Returned time
+ */
+
+TIME *
+td_get_time (TIME *t)
+{
+ time_t currtime;
+ struct tm *tm;
+
+ currtime = time(NULL);
+ tm = localtime(&currtime);
+
+ t->hour = tm->tm_hour;
+ t->minute = tm->tm_min;
+ t->second = tm->tm_sec;
+ t->ms = 0;
+
+ return (t);
+}
+
+
+DATE *
+td_get_date (DATE *d)
+{
+ time_t t = time(NULL);
+ struct tm *lt = localtime(&t);
+
+ d->year = lt->tm_year + 1900;
+ d->month = lt->tm_mon;
+ d->month++;
+ d->date = lt->tm_mday;
+ d->dow = lt->tm_wday;
+ d->dow++;
+
+ return (d);
+}
+
+double td_current_sec (void)
+{
+ TIME t;
+
+ td_get_time (&t);
+ return (td_time_to_sec (&t));
+}
+
+
+double
+td_time_to_sec (TIME *t)
+{
+ double ts;
+
+ ts = t->hour * 3600.;
+ ts += t->minute * 60.;
+ ts += t->second;
+ ts += t->ms / 1000.;
+
+ return (ts);
+}
+
+
+TIME *
+td_time_sub (const TIME *t1, const TIME *t2, TIME *tdiff)
+{
+ tdiff->hour = t1->hour - t2->hour;
+ tdiff->minute = t1->minute - t2->minute;
+ tdiff->second = t1->second - t2->second;
+ tdiff->ms = t1->ms - t2->ms;
+
+ return td_time_norm (tdiff);
+}
+
+
+TIME *
+td_time_add (const TIME *t1, const TIME *t2, TIME *tsum)
+{
+ tsum->ms = t1->ms + t2->ms;
+ tsum->second = t1->second + t2->second;
+ tsum->minute = t1->minute + t2->minute;
+ tsum->hour = t1->hour + t2->hour;
+
+ return td_time_norm (tsum);
+}
+
+
+TIME *
+td_time_copy (TIME *to, const TIME *from)
+{
+ to->ms = from->ms;
+ to->second = from->second;
+ to->minute = from->minute;
+ to->hour = from->hour;
+
+ return (to);
+}
+
+
+/* NAME
+ * td_time_norm Normalize time in structure
+ *
+ * SYNOPSIS
+ * t = td_time_norm (t)
+ * TIME *t Time to be normalized
+ */
+
+TIME *
+td_time_norm (TIME *t)
+{
+ while (t->ms < 0) {
+ t->ms += 1000;
+ t->second--;
+ }
+ while (t->second < 0) {
+ t->second += 60;
+ t->minute--;
+ }
+ while (t->minute < 0) {
+ t->minute += 60;
+ t->hour--;
+ }
+
+ while (t->ms > 1000) {
+ t->ms -= 1000;
+ t->second++;
+ }
+ while (t->second > 60) {
+ t->second -= 60;
+ t->minute++;
+ }
+ while (t->minute > 60) {
+ t->minute -= 60;
+ t->hour++;
+ }
+
+ return (t);
+}
+
+
+/* NAME
+ * td_get_tmdt Return current time and date
+ *
+ * SYNOPSIS
+ * get_tmdt (td)
+ * TIMEDATE *td Pointer to structure that holds time & date
+ */
+
+void
+td_get_tmdt (TIMEDATE *td)
+{
+ td_get_date (&td->d);
+ td_get_time (&td->t);
+}
+
+
+/* NAME
+ * td_str_tmdt Put time & date into string
+ *
+ * SYNOPSIS
+ * str = td_str_tmdt (td)
+ * TIMEDATE *td Pointer ot time & date structure
+ * char *str Pointer to output string (Width = 21 chars)
+ *
+ * NOTES
+ * str points to a area of memory devoted to this routine
+ * DO NOT make any changes to str
+ */
+
+const char *
+td_str_tmdt (const TIMEDATE *td)
+{
+ static char str[80];
+
+ strncpy (str, td_str_date(&td->d), sizeof(str));
+ strncat (str, " ", sizeof(str));
+ strncat (str, td_str_stime(&td->t), sizeof(str));
+
+ return (str);
+}
+
+
+/* NAME
+ * td_str_time Convert time into long string form
+ *
+ * SYNOPSIS
+ * str = td_str_time (t)
+ * char *str Output string (Field width = 12 chars)
+ * TIME *t Time to be converted
+ */
+
+const char *
+td_str_time (const TIME *t)
+{
+ static char str[80];
+ char am_pm;
+ int hour;
+
+ hour = t->hour;
+ if (hour < 12) {
+ am_pm = 'a';
+ if (hour == 0)
+ hour = 12;
+ } else if (hour >= 12) {
+ am_pm = 'p';
+ if (hour > 12)
+ hour -= 12;
+ }
+
+ snprintf (str, sizeof(str), "%2d:%02d:%02d.%03d%c",
+ hour, t->minute, t->second, t->ms, am_pm);
+
+ return (str);
+}
+
+
+
+/* NAME
+ * td_str_stime Convert time into short string form
+ *
+ * SYNOPSIS
+ * str = td_str_stime (t)
+ * char *str Output string (Field width = 6 chars)
+ * TIME *t Time to be converted
+ */
+
+const char *
+td_str_stime (const TIME *t)
+{
+ static char str[80];
+ char am_pm;
+ int hour;
+
+ hour = t->hour;
+ if (hour < 12) {
+ am_pm = 'a';
+ if (hour == 0)
+ hour = 12;
+ } else if (hour >= 12) {
+ am_pm = 'p';
+ if (hour > 12)
+ hour -= 12;
+ }
+
+ snprintf (str, sizeof(str), "%2d:%02d%c",
+ hour, t->minute, am_pm);
+
+ return (str);
+}
+
+
+/* NAME
+ * td_str_date Convert date to string form
+ *
+ * SYNOPSIS
+ * str = td_str_date (d)
+ * DATE *d Date to be converted
+ * char *str Pointer to output string (Width = 8 chars)
+ *
+ * NOTES
+ * DO NOT make any changes to str. It belongs to this routine
+ */
+
+const char *
+td_str_date (const DATE *d)
+{
+ static char str[80];
+
+ snprintf (str, sizeof(str), "%2d-%02d-%02d", d->month, d->date, d->year - 1900);
+
+ return (str);
+}
+
+/*-----------------------------------------------------------------------------
+ * NAME
+ * td_str_cdate Convert date to a character string
+ *
+ * SYNOPSIS
+ * str = td_str_cdate (d)
+ * DATE *d Date to convert
+ * char *str Pointer to date in character format
+ *
+ * DESCRIPTION
+ * The date is put in the form:
+ * <day name>, <month name> <date>, <year>
+ *
+ * Field width ranges from 17 to 29 characters
+ *
+ * NOTES
+ * str belongs to this routine, do NOT alter str
+ *---------------------------------------------------------------------------*/
+
+
+char *
+td_str_cdate (DATE *d)
+{
+ static char str[50];
+ char temp[50];
+
+ strcpy (str, "");
+
+ if (d->dow != 0) { /* only print day name if dow != 0 */
+ strcat (str, td_day_name(d->dow));
+ strcat (str, ", ");
+ }
+
+ snprintf (temp, sizeof(temp), "%s %d, %d", td_month_name(d->month), d->date, d->year);
+ strcat (str, temp);
+
+ return (str);
+}
+
+
+/*--------------------------------------------------------------*/
+/* td_month_name(int n) */
+/* return pointer to name of month given month number, n */
+/*--------------------------------------------------------------*/
+char *
+td_month_name ( /* return name of n-th month */
+ int n
+)
+{
+ static char *name[] = {
+ "",
+ "January",
+ "February",
+ "March",
+ "April",
+ "May",
+ "June",
+ "July",
+ "August",
+ "September",
+ "October",
+ "November",
+ "December"
+ };
+
+ return ((n < 1 || n > 12) ? name[0] : name[n]);
+}
+
+
+/*--------------------------------------------------------------*/
+/* td_day_name(int n) */
+/* return pointer to name of day given day number, n (1..7)*/
+/*--------------------------------------------------------------*/
+char *
+td_day_name (int n)
+{
+ static char *name[] = {
+ "",
+ "Monday",
+ "Tuesday",
+ "Wednesday",
+ "Thursday",
+ "Friday",
+ "Saturday",
+ "Sunday"
+ };
+
+ return ((n < 1 || n > 7) ? name[0] : name[n]);
+}
--- /dev/null
+/*****************************************************************************
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: xform.cpp,v 1.1 2000/06/19 02:58:08 kevin Exp $
+**
+** This program is free software; you can redistribute it and/or modify
+** it under the terms of the GNU General Public License (version 2) as
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+#include "kmath.h"
+
+/* NAME
+ * rotate2d Rotates an array of 2 dimensional vectors
+ *
+ * SYNOPSIS
+ * rotate2d (x, y, n, angle)
+ * double x[], y[] Array of points
+ * int n Number of points in array
+ * double angle Rotation angle (counter-clockwise)
+ */
+
+void
+rotate2d (double x[], double y[], int n, double angle)
+{
+ double cos_theta = cos (angle);
+ double sin_theta = sin (angle);
+
+ for (int i = 0; i < n; i++) {
+ double xrot = x[i] * cos_theta - y[i] * sin_theta;
+ double yrot = x[i] * sin_theta + y[i] * cos_theta;
+ x[i] = xrot;
+ y[i] = yrot;
+ }
+}
+
+
+/* NAME
+ * xlat2d Translates an array of 2 dimensional vectors
+ *
+ * SYNOPSIS
+ * xlat2d (x, y, n, xoffset, yoffset)
+ * double x[], y[] Array of points
+ * int n Number of points in array
+ * double xoffset, yoffset Offset to translate points by
+ */
+
+void
+xlat2d (double x[], double y[], int n, double xoffset, double yoffset)
+{
+ for (int i = 0; i < n; i++) {
+ x[i] += xoffset;
+ y[i] += yoffset;
+ }
+}
+
+
+/* NAME
+ * scale2d Scale an array of 2 dimensional vectors
+ *
+ * SYNOPSIS
+ * scale2d (x, y, n, xoffset, yoffset)
+ * double x[], y[] Array of points
+ * int n Number of points in array
+ * double xfact, yfact x & y scaling factors
+ */
+
+void
+scale2d (double x[], double y[], int n, double xfact, double yfact)
+{
+ for (int i = 0; i < n; i++) {
+ x[i] *= xfact;
+ y[i] *= yfact;
+ }
+}
+
+
+void
+indent_mtx2 (GRFMTX_2D m)
+{
+ m[0][0] = 1.0; m[0][1] = 0.0; m[0][2] = 0.0;
+ m[1][0] = 0.0; m[1][1] = 1.0; m[1][2] = 0.0;
+ m[2][0] = 0.0; m[2][1] = 0.0; m[2][2] = 1.0;
+}
+
+void
+xlat_mtx2 (GRFMTX_2D m, const double x, const double y)
+{
+ indent_mtx2 (m);
+ m[2][0] = x;
+ m[2][1] = y;
+}
+
+void
+scale_mtx2 (GRFMTX_2D m, const double sx, const double sy)
+{
+ indent_mtx2 (m);
+ m[0][0] = sx;
+ m[1][1] = sy;
+}
+
+void
+rot_mtx2 (GRFMTX_2D m, const double theta)
+{
+ double c = cos(theta);
+ double s = sin(theta);
+
+ indent_mtx2 (m);
+ m[0][0] = c; m[0][1] = s;
+ m[1][0] = -s; m[1][1] = c;
+}
+
+void
+mult_mtx2 (GRFMTX_2D m1, GRFMTX_2D m2, GRFMTX_2D result)
+{
+ GRFMTX_2D temp;
+
+ for (int row = 0; row < 3; row++)
+ for (int col = 0; col < 3; col++) {
+ temp[row][col] = 0;
+ for (int calc = 0; calc < 3; calc++)
+ temp[row][col] += m1[row][calc] * m2[calc][col];
+ }
+
+ for (int row = 0; row < 3; row++)
+ for (int col = 0; col < 3; col++)
+ result[row][col] = temp[row][col];
+}
+
+void
+xform_mtx2 (GRFMTX_2D m, double& x, double& y)
+{
+ double xt = x * m[0][0] + y * m[1][0] + m[2][0];
+ double yt = x * m[0][1] + y * m[1][1] + m[2][1];
+
+ x = xt;
+ y = yt;
+}
projects / ctsim.git / commitdiff