libctsim/scanner.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **   Name:          scanner.cpp
   5 **   Purpose:       Classes for CT scanner
   6 **   Programmer:    Kevin Rosenberg
   7 **   Date Started:  1984
   8 **
   9 **  This is part of the CTSim program
  10 **  Copyright (C) 1983-2000 Kevin Rosenberg
  11 **
  12 **  $Id: scanner.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $
  13 **
  14 **  This program is free software; you can redistribute it and/or modify
  15 **  it under the terms of the GNU General Public License (version 2) as
  16 **  published by the Free Software Foundation.
  17 **
  18 **  This program is distributed in the hope that it will be useful,
  19 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  20 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21 **  GNU General Public License for more details.
  22 **
  23 **  You should have received a copy of the GNU General Public License
  24 **  along with this program; if not, write to the Free Software
  25 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  26 ******************************************************************************/
  27
  28 #include "ct.h"
  29
  30
  31 // NAME
  32 //   DetectorArray       Construct a DetectorArray
  33
  34 DetectorArray::DetectorArray (const int nDet)
  35 {
  36   m_nDet = nDet;
  37   m_detValues = new DetectorValue [m_nDet];
  38 }
  39
  40
  41 // NAME
  42 //   ~DetectorArray             Free memory allocated to a detector array
  43
  44 DetectorArray::~DetectorArray (void)
  45 {
  46   delete [] m_detValues;
  47 }
  48
  49
  50
  51 /* NAME
  52  *   Scanner::Scanner           Construct a user specified detector structure
  53  *
  54  * SYNOPSIS
  55  *   Scanner (phm, nDet, nView, nSample)
  56  *   Phantom& phm               PHANTOM that we are making detector for
  57  *   int geomety                Geometry of detector
  58  *   int nDet                   Number of detector along detector array
  59  *   int nView                  Number of rotated views
  60  *   int nSample                Number of rays per detector
  61  */
  62
  63 Scanner::Scanner (const Phantom& phm, const ScannerGeometry geometry, int nDet, int nView, int nSample, const double rot_anglen)
  64 {
  65   m_phmLen = phm.maxAxisLength();      // maximal length along an axis
  66
  67   if (nView < 1)
  68     nView = 1;
  69   if (nSample < 1)
  70     m_nSample = 1;
  71   if (nDet < 1)
  72     nDet = 1;
  73   if ((nDet % 2) == 0)
  74     ++nDet;             // ensure odd number of detectors
  75
  76   m_geometry = geometry;
  77   m_nDet     = nDet;
  78   m_nView    = nView;
  79   m_nSample  = nSample;
  80   m_detLen   = SQRT2 * m_phmLen * ((m_nDet + N_EXTRA_DETECTORS) / static_cast<double>(m_nDet));
  81
  82   m_rotLen   =  rot_anglen;
  83
  84   m_radius  = m_detLen / 2;
  85   m_detInc  = m_detLen / m_nDet;
  86   m_rotInc  = m_rotLen / m_nView;
  87
  88   m_initPos.xd1 = m_detLen / 2;
  89   m_initPos.yd1 = -m_detLen / 2;
  90   m_initPos.xd2 = m_detLen / 2;
  91   m_initPos.yd2 = m_detLen / 2;
  92   m_initPos.xs1 = -m_detLen / 2;
  93   m_initPos.ys1 = -m_detLen / 2;
  94   m_initPos.xs2 = -m_detLen / 2;
  95   m_initPos.ys2 = m_detLen / 2;
  96   m_initPos.angle = 0.0;
  97 }
  98
  99 Scanner::~Scanner (void)
 100 {
 101 }
 102
 103
 104
 105 /* NAME
 106  *   raysum_collect                             Calculate ray sums for a Phantom
 107  *
 108  * SYNOPSIS
 109  *   rs = raysum_collect (det, phm, start_view, trace, unit_pulse)
 110  *   Scanner& det                       Scanner specifications**
 111  *   RAYSUM *rs                         Calculated ray sum matrix
 112  *   Phantom& phm                       Phantom we are taking ray sums of
 113  *   int trace                          Boolean flag to signal ray sum tracing
 114 */
 115
 116 void
 117 Scanner::collectProjections (Projections& proj, const Phantom& phm, const int start_view, const int trace)
 118 {
 119   GRFMTX_2D rotmtx_initial, temp;
 120   GRFMTX_2D rotmtx_incr;
 121
 122   double start_angle = start_view * proj.rotInc();
 123   double xcent = phm.xmin() + (phm.xmax() - phm.xmin()) / 2;
 124   double ycent = phm.ymin() + (phm.ymax() - phm.ymin()) / 2;
 125
 126   double xd1 = xcent + m_initPos.xd1;
 127   double yd1 = ycent + m_initPos.yd1;
 128   double xd2 = xcent + m_initPos.xd2;
 129   double yd2 = ycent + m_initPos.yd2;
 130   double xs1 = xcent + m_initPos.xs1;
 131   double ys1 = ycent + m_initPos.ys1;
 132   double xs2 = xcent + m_initPos.xs2;
 133   double ys2 = ycent + m_initPos.ys2;
 134
 135   m_trace = trace;
 136
 137 #ifdef HAVE_SGP
 138   SGP_ID gid;
 139   if (m_trace >= TRACE_PHM) {
 140     double wsize = 1.42 * m_phmLen / 2; /* sqrt(2) * radius */
 141
 142     gid = sgp2_init (512, 512, "RayCollect");
 143     sgp2_color (C_LTBLUE);
 144     sgp2_window (xcent - wsize, ycent - wsize, xcent + wsize, ycent + wsize);
 145     sgp2_color (C_BROWN);
 146 #if RADIUS
 147     sgp2_draw_circle (m_phmLen / 2);
 148 #else
 149     sgp2_draw_rect (xcent - m_phmLen / 2, ycent - m_phmLen / 2,
 150                     xcent + m_phmLen / 2, ycent + m_phmLen / 2);
 151 #endif
 152     sgp2_color (C_BROWN);
 153     sgp2_move_abs (0., 0.);
 154     sgp2_draw_circle (wsize);
 155     //      raysum_trace_menu_column = (crt->xsize * crt->asp) / 8 + 3;
 156     traceShowParam ("X-Ray Simulator", "%s", RAYSUM_TRACE_ROW_TITLE, 8+C_LTWHITE, " ");
 157     traceShowParam ("---------------", "%s", RAYSUM_TRACE_ROW_TITLE2, 8+C_LTWHITE, " ");
 158     traceShowParam ("Phantom:",       "%s", RAYSUM_TRACE_ROW_PHANT_ID, C_YELLOW, " Herman");
 159     traceShowParam ("Chomaticity  :", "%s", RAYSUM_TRACE_ROW_CHROMATIC, C_LTGREEN, "Mono");
 160     traceShowParam ("Scatter      :", "%5.1f", RAYSUM_TRACE_ROW_SCATTER, C_LTGREEN, 0.);
 161     traceShowParam ("Photon Uncert:", "%5.1f", RAYSUM_TRACE_ROW_PHOT_STAT, C_LTGREEN, 0.);
 162     traceShowParam ("Num Scanners:", "%5d", RAYSUM_TRACE_ROW_NDET, C_LTRED, proj.nDet());
 163     traceShowParam ("Num Views    :", "%5d", RAYSUM_TRACE_ROW_NVIEW, C_LTRED, proj.nView());
 164     traceShowParam ("Samples / Ray:", "%5d", RAYSUM_TRACE_ROW_SAMPLES, C_LTRED, m_nSample);
 165
 166     sgp2_color (C_LTGREEN);
 167     phm.draw();
 168
 169     initmarker (BDIAMOND, 129);
 170   }
 171 #endif
 172
 173 /* Calculate initial rotation matrix */
 174   xlat_mtx2 (rotmtx_initial, -xcent, -ycent);
 175   rot_mtx2 (temp, start_angle);
 176   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
 177   xlat_mtx2 (temp, xcent, ycent);
 178   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
 179
 180   xform_mtx2 (rotmtx_initial, xd1, yd1);        /* rotate detector endpoints */
 181   xform_mtx2 (rotmtx_initial, xd2, yd2);      /* to initial view_angle */
 182   xform_mtx2 (rotmtx_initial, xs1, ys1);
 183   xform_mtx2 (rotmtx_initial, xs2, ys2);
 184
 185 /* Calculate incrementatal rotation matrix */
 186   xlat_mtx2 (rotmtx_incr, -xcent, -ycent);
 187   rot_mtx2 (temp, proj.rotInc());
 188   mult_mtx2 (rotmtx_incr, temp, rotmtx_incr);
 189   xlat_mtx2 (temp, xcent, ycent);
 190   mult_mtx2 (rotmtx_incr, temp, rotmtx_incr);
 191
 192   int iview;
 193   double viewAngle;
 194   for (iview = 0, viewAngle = start_angle;  iview < proj.nView(); iview++, viewAngle += proj.rotInc()) {
 195       DetectorArray& detArray = proj.getDetectorArray( iview );
 196
 197 #ifdef HAVE_SGP
 198     if (m_trace >= TRACE_PHM) {
 199       sgp2_move_abs (xd1, yd1);
 200       sgp2_line_abs (xd2, yd2);
 201       sgp2_move_abs (xs1, ys1);
 202       sgp2_line_abs (xs2, ys2);
 203     }
 204 #endif
 205     if (m_trace)
 206       traceShowParam ("Current View :", "%5d", RAYSUM_TRACE_ROW_CURR_VIEW, C_LTMAGENTA, iview);
 207
 208     projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2);
 209     detArray.setViewAngle (viewAngle);
 210
 211 #ifdef HAVE_SGP
 212     if (m_trace >= TRACE_PHM) {
 213       //        rs_plot (detArray, xd1, yd1, xcent, ycent, theta);
 214       sgp2_move_abs (xd1, yd1);
 215       sgp2_line_abs (xd2, yd2);
 216       sgp2_move_abs (xs1, ys1);
 217       sgp2_line_abs (xs2, ys2);
 218     }
 219 #endif
 220     xform_mtx2 (rotmtx_incr, xd1, yd1);  // rotate detector endpoints
 221     xform_mtx2 (rotmtx_incr, xd2, yd2);
 222     xform_mtx2 (rotmtx_incr, xs1, ys1);
 223     xform_mtx2 (rotmtx_incr, xs2, ys2);
 224   } /* for each iview */
 225 }
 226
 227
 228 /* NAME
 229  *    rayview                   Calculate raysums for a view at any angle
 230  *
 231  * SYNOPSIS
 232  *    rayview (phm, detArray, xd1, nSample, yd1, xd2, yd2, xs1, ys1, xs2, ys2)
 233  *    Phantom& phm              Phantom to scan
 234  *    DETARRAY *detArray                Storage of values for detector array
 235  *    Scanner& det              Scanner parameters
 236  *    double xd1, yd1, xd2, yd2 Beginning & ending detector positions
 237  *    double xs1, ys1, xs2, ys2 Beginning & ending source positions
 238  *
 239  * RAY POSITIONING
 240  *         For each detector, have there are a variable number of rays traced.
 241  *     The source of each ray is the center of the source x-ray cell. The
 242  *     detector positions are equally spaced within the cell
 243  *
 244  *         The increments between rays are calculated so that the cells start
 245  *     at the beginning of a detector cell and they end on the endpoint
 246  *     of the cell.  Thus, the last cell starts at (xd2-ddx),(yd2-ddy).
 247  *         The exception to this is if there is only one ray per detector.
 248  *     In that case, the detector position is the center of the detector cell.
 249  */
 250
 251 void
 252 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)
 253 {
 254   double ddx = (xd2 - xd1) / detArray.nDet();  // change in coords between detectors
 255   double ddy = (yd2 - yd1) / detArray.nDet();
 256   double sdx = (xs2 - xs1) / detArray.nDet();  // change in coords between source
 257   double sdy = (ys2 - ys1) / detArray.nDet();
 258
 259   double ddx2 = ddx / m_nSample;        // Incr. between rays with detector cell
 260   double ddy2 = ddy / m_nSample;        // Doesn't include detector endpoints
 261   double ddx2_ofs = ddx2 / 2;           // offset of 1st ray from start of detector cell
 262   double ddy2_ofs = ddy2 / 2;
 263
 264   double xd_maj = xd1 + ddx2_ofs;       // Incr. between detector cells
 265   double yd_maj = yd1 + ddy2_ofs;
 266   double xs_maj = xs1 + (sdx / 2);      // put ray source in center of cell
 267   double ys_maj = ys1 + (sdy / 2);
 268
 269   DetectorValue* detval = detArray.detValues();
 270
 271   if (phm.getComposition() == P_UNIT_PULSE) {  // put unit pulse in center of view
 272     for (int d = 0; d < detArray.nDet(); d++)
 273       if (detArray.nDet() / 2 == d && (d % 2) == 1)
 274         detval[d] = 1;
 275       else
 276         detval[d] = 0;
 277   } else {
 278     for (int d = 0; d < detArray.nDet(); d++) {
 279       double xd = xd_maj;
 280       double yd = yd_maj;
 281       double xs = xs_maj;
 282       double ys = ys_maj;
 283       double sum = 0.0;
 284       for (int i = 0; i < m_nSample; i++) {
 285 #ifdef HAVE_SGP
 286         if (m_trace >= TRACE_RAYS) {
 287           sgp2_move_abs (xs, ys);
 288           sgp2_line_abs (xd, yd);
 289         }
 290 #endif
 291         sum += projectSingleLine (phm, xd, yd, xs, ys);
 292
 293         if (m_trace >= TRACE_RAYS)
 294           traceShowParam ("Attenuation  :", "%5.2f", RAYSUM_TRACE_ROW_ATTEN, C_LTMAGENTA, "sum");
 295
 296 #ifdef HAVE_SGP
 297         if (m_trace >= TRACE_RAYS) {
 298           sgp2_move_abs (xs, ys);
 299           sgp2_line_abs (xd, yd);
 300         }
 301 #endif
 302         xd += ddx2;
 303         yd += ddy2;
 304       }
 305
 306       detval[d] = sum / m_nSample;
 307       xd_maj += ddx;
 308       yd_maj += ddy;
 309       xs_maj += sdx;
 310       ys_maj += sdy;
 311     } /* for each detector */
 312   } /* if not unit pulse */
 313 }
 314
 315
 316 void
 317 Scanner::traceShowParam (const char *label, const char *fmt, int row, int color, ...)
 318 {
 319   char s[80];
 320   va_list arg;
 321
 322   va_start(arg, color);
 323   //  cio_set_cpos (raysum_trace_menu_column, row);
 324   snprintf (s, sizeof(s), label, "%s");
 325   //  cio_set_text_clr (color - 8, 0);
 326   cio_put_str (s);
 327   vsnprintf (s, sizeof(s), fmt, arg);
 328   //  cio_set_text_clr (color, 0);
 329   cio_put_str (s);
 330
 331   va_end(arg);
 332 }
 333
 334
 335 /* NAME
 336  *    projectSingleLine                 INTERNAL: Calculates raysum along a line for a Phantom
 337  *
 338  * SYNOPSIS
 339  *    rsum = phm_ray_attenuation (phm, x1, y1, x2, y2)
 340  *    double rsum               Ray sum of Phantom along given line
 341  *    Phantom& phm;             Phantom from which to calculate raysum
 342  *    double *x1, *y1, *x2, y2  Endpoints of ray path (in Phantom coords)
 343  */
 344
 345 double
 346 Scanner::projectSingleLine (const Phantom& phm, const double x1, const double y1, const double x2, const double y2)
 347 {
 348   // check ray against each pelem in Phantom
 349   double rsum = 0.0;
 350   for (PElemConstIterator i = phm.listPElem().begin(); i != phm.listPElem().end(); i++)
 351     rsum += projectLineAgainstPElem (**i, x1, y1, x2, y2);
 352
 353   return (rsum);
 354 }
 355
 356
 357 /* NAME
 358  *   pelem_ray_attenuation              Calculate raysum of an pelem along one line
 359  *
 360  * SYNOPSIS
 361  *   rsum = pelem_ray_attenuation (pelem, x1, y1, x2, y2)
 362  *   double rsum                Computed raysum
 363  *   PhantomElement& pelem              Pelem to scan
 364  *   double x1, y1, x2, y2      Endpoints of raysum line
 365  */
 366
 367 double
 368 Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double y1, double x2, double y2)
 369 {
 370   if (! pelem.clipLineWorldCoords (x1, y1, x2, y2)) {
 371     if (m_trace == TRACE_CLIPPING)
 372       cio_tone (1000., 0.05);
 373     return (0.0);
 374   }
 375
 376 #ifdef HAVE_SGP
 377   if (m_trace == TRACE_CLIPPING) {
 378     sgp2_move_abs (x1, y1);
 379     sgp2_line_abs (x2, y2);
 380     cio_tone (8000., 0.05);
 381     sgp2_move_abs (x1, y1);
 382     sgp2_line_abs (x2, y2);
 383   }
 384 #endif
 385
 386   double len = lineLength (x1, y1, x2, y2);
 387   return (len * pelem.atten());
 388 }
 389