libctsim/phm2image.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **   Name:         phm2image.cpp
   5 **   Purpose:      Convert phantom objects to rasterized images
   6 **   Programmer:   Kevin Rosenberg
   7 **   Date Started: 1984
   8 **
   9 **  This part of the CTSim program
  10 **  Copyright (C) 1983-2000 Kevin Rosenberg
  11 **
  12 **  $Id: phm2image.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  *    pic_to_imagefile          Make image array from Phantom
  33  *
  34  * SYNOPSIS
  35  *    pic_to_imagefile (pic, im, nsample)
  36  *    Phantom& pic              Phantom definitions
  37  *    ImageFile  *im            Computed pixel array
  38  *    int nsample               Number of samples along each axis for each pixel
  39  *                              (total samples per pixel = nsample * nsample)
  40  */
  41
  42 void
  43 phm_to_imagefile (const Phantom& phm, ImageFile& im, const int col_start, const int col_count, const int in_nsample, const int trace)
  44 {
  45   int nsample = in_nsample;
  46   double dx = phm.xmax() - phm.xmin();
  47   double dy = phm.ymax() - phm.ymin();
  48   double xcent = phm.xmin() + dx / 2;
  49   double ycent = phm.ymin() + dy / 2;
  50   double phmlen = (dx > dy ? dx : dy);
  51
  52   double phmradius = phmlen / 2;
  53
  54   double xmin = xcent - phmradius;
  55   double xmax = xcent + phmradius;
  56   double ymin = ycent - phmradius;
  57   double ymax = ycent + phmradius;
  58
  59   im.setAxisExtent (xmin, xmax, ymin, ymax);
  60   int nx = im.nx();
  61   int ny = im.ny();
  62
  63 /* Each pixel holds the average of the intensity of the cell with (ix,iy)
  64  * at the center of the pixel
  65  *
  66  * Set major increments so that the last cell v[nx-1][ny-1] will start at
  67  * (xmax - xinc, ymax - yinc).
  68  *
  69  * Set minor increments so that sample points are centered in cell
  70  */
  71
  72   if (nx < 2 || ny < 2)
  73       return;
  74
  75   double xinc = (xmax - xmin) / nx;
  76   double yinc = (ymax - ymin) / ny;
  77   im.setAxisIncrement (xinc, yinc);
  78
  79   if (nsample < 1)
  80          nsample = 1;
  81
  82   double kxinc = xinc / nsample;                /* interval between samples */
  83   double kyinc = yinc / nsample;
  84   double kxofs = kxinc / 2;             /* offset of 1st point */
  85   double kyofs = kyinc / 2;
  86
  87 #ifdef HAVE_SGP
  88   SGP_ID gid = NULL;
  89   if (trace > TRACE_TEXT) {
  90     gid = sgp2_init (0, 0, "Phantom to Image");
  91     sgp2_window (xmin, ymin, xmax, ymax);
  92     sgp2_frame_vpt();
  93   }
  94 #endif
  95
  96   ImageFileArray v = im.getArray();
  97
  98   double x_start = xmin + (col_start * xinc);
  99   for (PElemConstIterator i = phm.listPElem().begin(); i != phm.listPElem().end(); i++) {
 100 #ifdef HAVE_SGP
 101     if (trace > TRACE_TEXT)
 102       sgp2_polyline_abs ((*i)->xOutline(), (*i)->yOutline(), (*i)->nOutlinePoints());
 103 #endif
 104     double x, y, xi, yi;
 105     int ix, iy, kx, ky;
 106     for (ix = 0, x = x_start; ix < col_count; ix++, x += xinc) {
 107       for (iy = 0, y = ymin; iy < ny; iy++, y += yinc) {
 108         for (kx = 0, xi = x + kxofs; kx < nsample; kx++, xi += kxinc) {
 109           for (ky = 0, yi = y + kyofs; ky < nsample; ky++, yi += kyinc)
 110             if ((*i)->isPointInside (xi, yi, PHM_COORD) == TRUE)
 111               v[ix][iy] += (*i)->atten();
 112         } // for kx
 113       } /* for iy */
 114     }  /* for ix */
 115   }  /* for pelem */
 116
 117
 118   if (nsample > 0) {
 119     double factor = 1.0 / (nsample * nsample);
 120
 121     for (int ix = 0; ix < col_count; ix++)
 122       for (int iy = 0; iy < ny; iy++)
 123         v[ix][iy] *= factor;
 124   }
 125
 126 #if HAVE_SGP
 127   if (trace > TRACE_TEXT)
 128       sgp2_close(gid);
 129 #endif
 130
 131 }
 132