From: Kevin M. Rosenberg Date: Mon, 19 Jun 2000 02:59:34 +0000 (+0000) Subject: r99: *** empty log message *** X-Git-Tag: debian-4.5.3-3~918 X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=commitdiff_plain;h=99dd1d6ed10db1f669a5fe6af71225a50fc0ddfb r99: *** empty log message *** --- diff --git a/libctsim/Makefile.am b/libctsim/Makefile.am new file mode 100644 index 0000000..327a13c --- /dev/null +++ b/libctsim/Makefile.am @@ -0,0 +1,6 @@ +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 + diff --git a/libctsim/Makefile.nt b/libctsim/Makefile.nt new file mode 100644 index 0000000..0760618 --- /dev/null +++ b/libctsim/Makefile.nt @@ -0,0 +1,35 @@ +# Makefile for libir + +!include + +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 diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp new file mode 100644 index 0000000..d05bcca --- /dev/null +++ b/libctsim/backprojectors.cpp @@ -0,0 +1,383 @@ +/***************************************************************************** +** 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(new BackprojectTrig (proj, im, interpType)); + else if (bjType == O_BPROJ_TABLE) + bj = static_cast(new BackprojectTable (proj, im, interpType)); + else if (bjType == O_BPROJ_DIFF) + bj = static_cast(new BackprojectDiff (proj, im, interpType)); + else if (bjType == O_BPROJ_DIFF2) + bj = static_cast(new BackprojectDiff2 (proj, im, interpType)); + else if (bjType == O_BPROJ_IDIFF2) + bj = static_cast(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 (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(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(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(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(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(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(L))]); +#endif + if (interpType == I_NEAREST) { + int iDetPos = iDetCenter + nearest (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(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 (xInc * sin (theta) / detInc * scale); + long int det_dy = nearest (yInc * cos (theta) / detInc * scale); + + // calculate L for first point in image (0, 0) + long int detPosColStart = nearest (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 +} diff --git a/libctsim/convolve.cpp b/libctsim/convolve.cpp new file mode 100644 index 0000000..52abd94 --- /dev/null +++ b/libctsim/convolve.cpp @@ -0,0 +1,84 @@ +/***************************************************************************** +** FILE IDENTIFICATION +** +** Name: Convolve.c Convolution functions +** Date Started: 13 Mar 86 +** +** This is part of the CTSim program +** Copyright (C) 1983-2000 Kevin Rosenberg +** +** This program is free software; you can redistribute it and/or modify +** it under the terms of the GNU General Public License (version 2) as +** published by the Free Software Foundation. +** +** This program is distributed in the hope that it will be useful, +** but WITHOUT ANY WARRANTY; without even the implied warranty of +** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +** GNU General Public License for more details. +** +** You should have received a copy of the GNU General Public License +** along with this program; if not, write to the Free Software +** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +******************************************************************************/ + +#include "ct.h" + + +/* NAME + * convolve Discrete convolution of two functions + * + * SYNOPSIS + * r = convolve (f1, f2, dx, n, np, func_type) + * double r Convolved result + * double f1[], f2[] Functions to be convolved + * double dx Difference between successive x values + * int n Array index to center convolution about + * int np Number of points in f1 array + * int func_type EVEN or ODD or EVEN_AND_ODD function f2 + * + * NOTES + * f1 is the projection data, its indices range from 0 to np - 1. + * The index for f2, the filter, ranges from -(np-1) to (np-1). + * There are 3 ways to handle the negative vertices of f2: + * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n]. + * All filters used in reconstruction are even. + * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n] + * 3. If f2 is both ODD AND EVEN, then we must store the value of f2 + * for negative indices. Since f2 must range from -(np-1) to (np-1), + * if we add (np - 1) to f2's array index, then f2's index will + * range from 0 to 2 * (np - 1), and the origin, x = 0, will be + * stored at f2[np-1]. + */ + +double +convolve (const double f1[], const double f2[], const double dx, const int n, const int np, const FunctionSymmetry func_type) +{ + double sum = 0.0; + + if (func_type == FUNC_BOTH) { +#if UNOPTIMIZED_CONVOLUTION + for (int i = 0; i < np; i++) + sum += f1[i] * f2[n - i + (np - 1)]; +#else + f2 += n + (np - 1); + for (int i = 0; i < np; i++) + sum += *f1++ * *f2--; +#endif + } else if (func_type == FUNC_EVEN) { + for (int i = 0; i < np; i++) { + int k = abs (n - i); + sum += f1[i] * f2[k]; + } + } else if (func_type == FUNC_ODD) { + for (int i = 0; i < np; i++) { + int k = n - i; + if (k < 0) + sum -= f1[i] * f2[k]; + else + sum += f1[i] * f2[k]; + } + } else + sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type); + + return (sum * dx); +} diff --git a/libctsim/dialogs.cpp b/libctsim/dialogs.cpp new file mode 100644 index 0000000..de12d98 --- /dev/null +++ b/libctsim/dialogs.cpp @@ -0,0 +1,212 @@ +/***************************************************************************** +** 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); +} + + + + diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp new file mode 100644 index 0000000..30d4718 --- /dev/null +++ b/libctsim/filter.cpp @@ -0,0 +1,319 @@ +/***************************************************************************** +** 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); +} + + diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp new file mode 100644 index 0000000..8756bc8 --- /dev/null +++ b/libctsim/imagefile.cpp @@ -0,0 +1,171 @@ +/***************************************************************************** +** 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 +} + + + + diff --git a/libctsim/options.cpp b/libctsim/options.cpp new file mode 100644 index 0000000..9bbb585 --- /dev/null +++ b/libctsim/options.cpp @@ -0,0 +1,271 @@ +/***************************************************************************** +** 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(-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(-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(-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(-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"); +} + + diff --git a/libctsim/phantom.cpp b/libctsim/phantom.cpp new file mode 100644 index 0000000..1508bea --- /dev/null +++ b/libctsim/phantom.cpp @@ -0,0 +1,724 @@ +/***************************************************************************** +** 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 (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); +} + + diff --git a/libctsim/phm2image.cpp b/libctsim/phm2image.cpp new file mode 100644 index 0000000..3e65e58 --- /dev/null +++ b/libctsim/phm2image.cpp @@ -0,0 +1,132 @@ +/***************************************************************************** +** 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 + +} + diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp new file mode 100644 index 0000000..8d36895 --- /dev/null +++ b/libctsim/projections.cpp @@ -0,0 +1,459 @@ +/***************************************************************************** +** 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); +} + + diff --git a/libctsim/reconstr.cpp b/libctsim/reconstr.cpp new file mode 100644 index 0000000..11b9fc5 --- /dev/null +++ b/libctsim/reconstr.cpp @@ -0,0 +1,181 @@ +/***************************************************************************** +** 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); +} + diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp new file mode 100644 index 0000000..1d7eaf7 --- /dev/null +++ b/libctsim/scanner.cpp @@ -0,0 +1,389 @@ +/***************************************************************************** +** 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(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()); +} + diff --git a/libctsupport/Makefile.am b/libctsupport/Makefile.am new file mode 100644 index 0000000..5fb10ba --- /dev/null +++ b/libctsupport/Makefile.am @@ -0,0 +1,7 @@ +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 + + + diff --git a/libctsupport/Makefile.nt b/libctsupport/Makefile.nt new file mode 100644 index 0000000..f85864e --- /dev/null +++ b/libctsupport/Makefile.nt @@ -0,0 +1,32 @@ +# Makefile for libk + +!include + +CC=cl +LD=link +CFLAGS=-O -nologo -I..\include -DENDIAN_HIGH=0 -DENDIAN_LOW=1 +LDFLAGS= +O=.obj + +# variables +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) + +all: libk.lib + +.obj: .c + $(CC) -c $(cvarsdll) $(CFLAGS) $*.c + + +libk.lib: $(OBJ1) + echo something to del > libk.lib + del libk.lib + lib /out:libk.lib $(OBJ1) + + +clean: + echo dummy > a.obj + echo dummy > a.exe + echo dummy > a.lib + del *.obj + del *.exe + del *.lib diff --git a/libctsupport/audio.cpp b/libctsupport/audio.cpp new file mode 100644 index 0000000..ad46a6d --- /dev/null +++ b/libctsupport/audio.cpp @@ -0,0 +1,58 @@ +/***************************************************************************** +** 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 +#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 +} diff --git a/libctsupport/byteorder.cpp b/libctsupport/byteorder.cpp new file mode 100644 index 0000000..f7734fa --- /dev/null +++ b/libctsupport/byteorder.cpp @@ -0,0 +1,342 @@ +#if HAVE_CONFIG_H +#include "config.h" +#endif + +#if HAVE_UNISTD_H +#include +#endif + +#include "kstddef.h" +#include "byteorder.h" + +inline void +SwapBytes2 (void* buffer) +{ + unsigned char* p = static_cast(buffer); + unsigned char temp = p[0]; + p[0] = p[1]; + p[1] = temp; +} + +// 0<->3 1<->2 = 0123 -> 3210 +inline void +SwapBytes4 (void* buffer) +{ + unsigned char* p = static_cast(buffer); + unsigned char temp = p[0]; + p[0] = p[3]; + p[3] = temp; + temp = p[1]; + p[1] = p[2]; + p[2] = temp; +} + +// 0<->7 1<->6 2<->5 3<->4 = 01234567 -> 76543210 +inline void +SwapBytes8 (void* buffer) +{ + unsigned char* p = static_cast(buffer); + unsigned char temp = p[0]; + p[0] = p[7]; + p[7] = temp; + temp = p[1]; + p[1] = p[6]; + p[6] = temp; + temp = p[2]; + p[2] = p[5]; + p[5] = temp; + temp = p[3]; + p[3] = p[4]; + p[4] = temp; +} + +void +ConvertNetworkOrder (void* buffer, size_t bytes) +{ +#if ! defined (WORDS_BIGENDIAN) + if (bytes < 2) + return; + + char* start = static_cast(buffer); + char* end = start + bytes - 1; // last byte + size_t nSwap = bytes / 2; + + while (nSwap-- > 0) { + unsigned char c = *start; + *start++ = *end; + *end-- = c; + } +#endif +} + +void +ConvertReverseNetworkOrder (void* buffer, size_t bytes) +{ +#if defined (WORDS_BIGENDIAN) + if (bytes < 2) + return; + + char* start = static_cast(buffer); + char* end = start + bytes - 1; // last byte + size_t nSwap = bytes / 2; + + while (nSwap-- > 0) { + unsigned char c = *start; + *start++ = *end; + *end-- = c; + } +#endif +} + +bool +write_nint16 (kuint16 const *n_in, int fd) +{ + kuint16 n = *n_in; + +#ifndef WORDS_BIGENDIAN + SwapBytes2 (&n); +#endif + + if (write (fd, &n, 2) != 2) + return false; + else + return true; +} + +bool +read_nint16 (kuint16 *n_out, int fd) +{ + if (read (fd, n_out, 2) != 2) + return false; + +#ifndef WORDS_BIGENDIAN + SwapBytes2 (n_out); +#endif + + return true; +} + +bool +write_nint32 (kuint32 const *n_in, int fd) +{ + kuint32 n = *n_in; + +#ifndef WORDS_BIGENDIAN + SwapBytes4(&n); +#endif + + if (write (fd, &n, 4) != 4) + return false; + else + return true; +} + +bool +read_nint32 (kuint32 *n_out, int fd) +{ + if (read (fd, n_out, 4) != 4) + return false; + +#ifndef WORDS_BIGENDIAN + SwapBytes4 (n_out); +#endif + + return true; +} + +bool +write_nfloat32 (kfloat32 const *f_in, int fd) +{ + kfloat32 f = *f_in; + +#ifndef WORDS_BIGENDIAN + SwapBytes4 (&f); +#endif + + if (write (fd, &f, 4) != 4) + return false; + else + return true; +} + +bool +read_nfloat32 (kfloat32 *f_out, int fd) +{ + if (read (fd, f_out, 4) != 4) + return false; + +#ifndef WORDS_BIGENDIAN + SwapBytes4(f_out); +#endif + + return true; +} + +bool +write_nfloat64 (kfloat64 const *f_in, int fd) +{ + kfloat64 f = *f_in; + +#ifndef WORDS_BIGENDIAN + SwapBytes8 (&f); +#endif + + if (write (fd, &f, 8) != 8) + return false; + else + return true; +} + +bool +read_nfloat64 (kfloat64 *f_out, int fd) +{ + if (read (fd, f_out, 8) != 8) + return false; + +#ifndef WORDS_BIGENDIAN + SwapBytes8 (f_out); +#endif + + return true; +} + + + +onetorderstream& onetorderstream::writeInt16 (kuint16 n) { +#ifndef WORDS_BIGENDIAN + SwapBytes2 (&n); +#endif + write (&n, 2); + return (*this); +} + +inetorderstream& inetorderstream::readInt16 (kuint16& n) { + read (&n, 2); +#ifndef WORDS_BIGENDIAN + SwapBytes2 (&n); +#endif + return (*this); +} + +onetorderstream& onetorderstream::writeInt32 (kuint32 n) { +#ifndef WORDS_BIGENDIAN + SwapBytes4(&n); +#endif + write (&n, 4); + return (*this); +} + +inetorderstream& inetorderstream::readInt32 (kuint32& n) { + read (&n, 4); +#ifndef WORDS_BIGENDIAN + SwapBytes4 (&n); +#endif + return (*this); +} + +onetorderstream& onetorderstream::writeFloat32 (kfloat32 n) { +#ifndef WORDS_BIGENDIAN + SwapBytes4 (&n); +#endif + write (&n, 4); + return (*this); +} + +inetorderstream& inetorderstream::readFloat32 (kfloat32& n) { + read (&n, 4); +#ifndef WORDS_BIGENDIAN + SwapBytes4 (&n); +#endif + return (*this); +} + +onetorderstream& onetorderstream::writeFloat64 (kfloat64 n) { +#ifndef WORDS_BIGENDIAN + SwapBytes8 (&n); +#endif + write (&n, 8); + return (*this); +} + +inetorderstream& inetorderstream::readFloat64 (kfloat64& n) { + read (&n, 8); +#ifndef WORDS_BIGENDIAN + SwapBytes8 (&n); +#endif + return (*this); +} + + + +void +write_rnint16 (kuint16 n, ostream& ostr) +{ +#ifdef WORDS_BIGENDIAN + SwapBytes2 (&n); +#endif + ostr.write (&n, 2); +} + +void +read_rnint16 (kuint16& n, istream& istr) +{ + istr.read (&n, 2); +#ifdef WORDS_BIGENDIAN + SwapBytes2 (&n); +#endif +} + +void +write_rnint32 (kuint32 n, ostream& ostr) +{ +#ifdef WORDS_BIGENDIAN + SwapBytes4(&n); +#endif + ostr.write (&n, 4); +} + +void +read_rnint32 (kuint32& n, istream istr) +{ + istr.read (&n, 4); +#ifdef WORDS_BIGENDIAN + SwapBytes4 (&n); +#endif +} + +void +write_rnfloat32 (kfloat32 n, ostream ostr) +{ +#ifdef WORDS_BIGENDIAN + SwapBytes4 (&n); +#endif + ostr.write (&n, 4); +} + +void +read_rnfloat32 (kfloat32& n, istream istr) +{ + istr.read (&n, 4); +#ifdef WORDS_BIGENDIAN + SwapBytes4 (&n); +#endif +} + +void +write_rnfloat64 (kfloat64 n, ostream ostr) +{ +#ifdef WORDS_BIGENDIAN + SwapBytes8 (&n); +#endif + ostr.write (&n, 8); +} + +void +read_rnfloat64 (kfloat64& n, istream istr) +{ + istr.read (&n, 8); +#ifdef WORDS_BIGENDIAN + SwapBytes8 (&n); +#endif +} + diff --git a/libctsupport/clip.cpp b/libctsupport/clip.cpp new file mode 100644 index 0000000..4aa6f01 --- /dev/null +++ b/libctsupport/clip.cpp @@ -0,0 +1,401 @@ +/***************************************************************************** +** 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); +} diff --git a/libctsupport/consoleio.cpp b/libctsupport/consoleio.cpp new file mode 100644 index 0000000..93d3160 --- /dev/null +++ b/libctsupport/consoleio.cpp @@ -0,0 +1,149 @@ +/***************************************************************************** +** 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 +#include +#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); +} + + diff --git a/libctsupport/filefuncs.cpp b/libctsupport/filefuncs.cpp new file mode 100644 index 0000000..100960a --- /dev/null +++ b/libctsupport/filefuncs.cpp @@ -0,0 +1,131 @@ +/***************************************************************************** +** 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 +#include +#include +#include +#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: - Retry | - New name | - 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); +} diff --git a/libctsupport/normangle.cpp b/libctsupport/normangle.cpp new file mode 100644 index 0000000..b04ee93 --- /dev/null +++ b/libctsupport/normangle.cpp @@ -0,0 +1,42 @@ +/***************************************************************************** +** 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); +} diff --git a/libctsupport/simpson.cpp b/libctsupport/simpson.cpp new file mode 100644 index 0000000..234e062 --- /dev/null +++ b/libctsupport/simpson.cpp @@ -0,0 +1,62 @@ +/***************************************************************************** +** 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); +} diff --git a/libctsupport/strfuncs.cpp b/libctsupport/strfuncs.cpp new file mode 100644 index 0000000..8033df4 --- /dev/null +++ b/libctsupport/strfuncs.cpp @@ -0,0 +1,168 @@ +/***************************************************************************** +** 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 +#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 + + diff --git a/libctsupport/syserror.cpp b/libctsupport/syserror.cpp new file mode 100644 index 0000000..dc3851c --- /dev/null +++ b/libctsupport/syserror.cpp @@ -0,0 +1,144 @@ +/***************************************************************************** +** 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 +#include +#include +#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; +} + diff --git a/libctsupport/timedate.cpp b/libctsupport/timedate.cpp new file mode 100644 index 0000000..0b920e0 --- /dev/null +++ b/libctsupport/timedate.cpp @@ -0,0 +1,415 @@ +/***************************************************************************** +** 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 +#include "kstddef.h" +#include + +/* 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: + * , , + * + * 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]); +} diff --git a/libctsupport/xform.cpp b/libctsupport/xform.cpp new file mode 100644 index 0000000..8d58f76 --- /dev/null +++ b/libctsupport/xform.cpp @@ -0,0 +1,148 @@ +/***************************************************************************** +** 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; +}