PinholeSPECTUB_Tools.h

/*
    Copyright (C) 2022, Matthew Strugari
    Copyright (C) 2014, Biomedical Image Group (GIB), Universitat de Barcelona, Barcelona, Spain. All rights reserved.
    Copyright (C) 2014, 2021, University College London
    This file is part of STIR.

    SPDX-License-Identifier: Apache-2.0

    See STIR/LICENSE.txt for details

    \author Carles Falcon
    \author Matthew Strugari
*/

#ifndef _WM_SPECT_mph_H
#define _WM_SPECT_mph_H

#include <string>
#include <vector>

namespace SPECTUB_mph
{

//::: structures ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

//.....structure for distribution function information

typedef struct
{
  int max_dimx; // maximum size in bins (for allocation purpose)
  int max_dimz; // maximum size in bins (for allocation purpose)

  int dimx; // actual number of bins forming the PSF (dimx)
  int dimz; // actual number of bins forming the PSF (dimz)

  int ib0; // first i index for the bins of the PSF
  int jb0; // first j index for the bins of the PSF

  float lngxcmd2; // half of the lenght in cm, dimx
  float lngzcmd2; // half of the lenght in cm, dimz

  float xc;
  float zc;

  float sum; // sum of values (for intensity normalization)

  float** val; // array of values
               // float res1;     // resolution dimx
  // float res2;     // resolution dimz

} psf2d_type;

//... collimator's holes parametre structure................................

typedef struct
{
  // int idet;        // index of the detel the hole projects to

  //... coordinates ......................

  float acy;  // angular coordinate of the hole in volume framework: angle (deg) (cyl)
  float acyR; // angular coordinate of the hole in rotated framework: ang (deg) (cyl)

  float x1; // x coordinate of the center of the hole (cm): rotated framework
  float y1; // y coordinate of the center of the hole (cm): rotated framework
  float z1; // z coordinate of the center of the hole (cm): rotated framework

  //... axial, edge and acceptance angles ...............................

  float ahx; // x angle of the hole axis (deg)
  float ahz; // z angle of the hole axis (deg)

  float aa_x; // acceptance angle x (aperture with respect to hole axis: half of the total acceptance angle)
  float aa_z; // acceptance angle z (aperture with respect to hole axis: half of the total acceptance angle)

  float egx; // absolute angle of the hole edge (x, minor) in rotated framework
  float Egx; // absolute angle of the hole edge (x, major) in rotated framework
  float egz; // absolute angle of the hole edge (z, minor) in rotated framework
  float Egz; // absolute angle of the hole edge (z, major) in rotated framework

  float ax_M; // maximum angle of acceptance x
  float ax_m; // minimum angle of acceptance x
  float az_M; // maximum angle of acceptance z
  float az_m; // minimum angle of acceptance z

  //... others ....................

  std::string shape; // hole shape { rect, round }
  bool do_round;     // true: round shape || false: rectangular shape
  float dxcm;        // horizontal size of the hole (cm): horizontal axis, diameter
  float dzcm;        // vertical size of the hole (cm): vertical axis, diameter

} hole_type;

//... structure for collimator information

typedef struct
{
  std::string model; // cylindrical (cyl) or polygonal prism (pol)

  float rad; // radius of cylinder containig holes (cyl) or apothem (pol)
  float L;   // collimator thickness
  float Ld2; // half of the collimator thickness

  int Nht;                      // total number of holes
  std::vector<hole_type> holes; // array hole_type structure

} mphcoll_type;

//.....structure for ring elements information

typedef struct
{
  int nh;               // number of holes that project to this detel
  std::vector<int> who; // array of indices of holes projecting to this detel

  float x0; // cartesian coordinates of the center (unrotated system of reference): x
  float y0; // cartesian coordinates of the center (unrotated system of reference): y
  float z0; // cartesian coordinates of the center (unrotated system of reference): z

  float xbin0; // x coordinate for the first bin of the detection row corresponding to this angle
  float ybin0; // y coordinate for the first bin of the detection row corresponding to this angle
  float zbin0; // z coordinate for the first bin of the detection row corresponding to this angle

  float incx; // increment in x to the following bin in detector row
  float incy; // increment in y to the following bin in detector row
  float incz; // increment in z to the following detector row

  float theta; // theta: in-plane angle radius vs x-axis. longitude (deg)
  float costh; // cosinus of theta
  float sinth; // sinus of theta

} detel_type;

//.....structure for distribution function information

typedef struct
{
  int dim;    // length (in discretization intervals) (odd number)
  int i_max;  // last i-index = dim -1
  float res;  // spatial resolution of distfunc (discretization interval)
  float* val; // double array of values

} discrf1d_type;

//.....structure for distribution function information

typedef struct
{
  int dim;     // length (in discretization intervals) (odd number)
  int i_max;   // last i-index = dim -1
  int j_max;   // last j-index = dim -1
  float res;   // spatial resolution of distfunc (discretization interval)
  float** val; // double array of values

} discrf2d_type;

//.....structure for projection information

typedef struct
{
  int Nbin; // number of bins per row
  int Nsli; // number of slices

  int Ndt; // number of detels (detector elements)
  int Nbd; // number of bins per detel
  int Nbt; // total number of bins

  float szcm;   // bin size in cm
  float szcmd2; // bin size in cm divided by 2
  float thcm;   // slice thickness (cm)
  float thcmd2; // slice thickness in cm divided by 2

  float crth;      // crystal thickness (cm) to correction for depth
  float crth_2;    // power 2 of the above value
  float crattcoef; // attenuation coefficient of crystal
  float max_dcr;   // maximum distance of a ray inside the crystal

  float FOVxcmd2; // FOVcmx divided by 2
  float FOVzcmd2; // FOVcmz divided by 2

  float rad;  // ring radius
  float radc; // extended ring radius = radius + crystal thickness

  int NOS;  // number of subsets
  int NdOS; // Number of detels per subset = Ndt/NOS
  int NbOS; // total number of bins per subset = Nbt/NOS

  int* order;  // order of the angles of projection (array formed by indexs of angles belonging to consecutive subsets)
  float sgm_i; // sigma of intrinsic PSF (cm)

  float* val;

} prj_mph_type;

//... structure for bin information..................................

typedef struct
{

  int Dimx; // number of columns
  int Dimy; // number of rows
  int Dimz; // number of slices

  int Npix; // number of pixels (voxels) per axial slice
  int Nvox; // number of voxels (the whole volume)

  int first_sl; // first slice to reconstruct (0->Nslic-1)
  int last_sl;  // last slice to reconstruct + 1 (end of the 'for' loop) (1->Nslic)

  float FOVxcmd2; // half of the size of the volume, dimension x (cm);
  float FOVcmyd2; // half of the size of the volume, dimension y (cm);
  float FOVzcmd2; // half of the size of the volume, dimension z (cm);

  float szcm; // voxel size (side length in cm)
  float thcm; // voxel thickness (cm)

  float x0; // x coordinate (cm, ref center of volume) of the first voxel
  float y0; // y coordinate (cm, ref center of volume) of the first voxel
  float z0; // z coordinate (cm, ref center of volume) of the first voxel

  float* val; // array of values

} volume_type;

//... matrix header information ............................

typedef struct
{
  int subsamp;   // bin subsampling factor for accurate PSF and convolution calculations (typically 2 to 5)
  float mn_w;    // minimum weight to be taken into account
  float highres; // high spatial resolution of continous distributions in PSF calculation

  float Nsigm;  // number of sigmas in PSF calculation
  float mndvh2; // squared minimum distance voxel-hole (cm). Reference for efficiency
  float ro;     // radius of the object

  float max_hsxcm; // maximum hole size dimx (for allocation purposes)
  float max_hszcm; // maximum hole size dimz (for allocation purposes)
  float max_amp;   // maximum amplification
  float tmax_aix;  // tangent of the maximum incidence angle x (for allocation purposes)
  float tmax_aiz;  // tangent of the maximum incidence angle z (for allocation purposes)

  bool do_psfi;     // true: correct for intrinsic PSF
  bool do_depth;    // true: correct for impact depth
  bool do_att;      // true: correct for attenuation
  bool do_full_att; // true: coef att for each bin of PSF || false: same att factor for all bins of PSF (central line)
  bool do_msk_att;  // true: to use att map as a mask
  bool do_msk_file; // true: explicit mask

  // internal booleans variables

  bool do_subsamp;
  bool do_round_cumsum;
  bool do_square_cumsum;

  std::string att_fn;      // attenuation map filename
  std::string msk_fn;      // explicit mask filename
  std::string detector_fn; // ring parameter filename
  std::string collim_fn;   // collimator parameter filename

  volume_type vol;               // structure with information of volume
  std::vector<detel_type> detel; // structure with detection elements information
  prj_mph_type prj;              // structure with detection rings information
  mphcoll_type collim;           // structure with the collimator information

} wmh_mph_type;

//.......weight_mat structure definition. Structure for reading weight matrix

typedef struct
{
  //... weight matrix dimensions .....................

  int ne;   // nonzero elements
  int Nbt;  // dimension 2 (rows) of the weight matrix (NbOS or Nbt)
  int Nvox; // dimension 1 (columns) of the weight matrix (Nvox)

  //... weight matrix values .........................

  float* ar; // array of nonzero elements of weight matrix (by rows)
  int* ja;   // array of the column index of the above elements
  int* ia;   // array containing the indexes of the previous vector where a row change happens

} wm_type;

//.......weight_mat_da structure definition. Structure for generating weight matrix

typedef struct
{
  int Nbt;     // dimension 2 (rows) of the weight matrix (NbOS or Nbt)
  int Nvox;    // dimension 1 (columns) of the weight matrix (Nvox)
  float** val; // double array to store weights (index of the projection element, number of weight for that element)
  int** col;   // double array to store column indexs of the above element  (index of the projection element, number of weight for
               // that element)
  int* ne;     // array indicating how many elements has been stored for each element of projection

  //... filename .............................................

  std::string fn;     // matrix name
  std::string fn_hdr; // matrix header file name

  //... format ...............................................

  bool do_save_STIR; // to save weight matrix with STIR format

  //... indexs for STIR format ...............................

  int *na, *nb, *ns;       // indexs for projection elements (angle, bin, slice respec.)
  short int *nx, *ny, *nz; // indexs for image elements (x,y,z)

} wm_da_type;

//.....structure for distribution function information

typedef struct
{
  discrf2d_type square; // distribution for square shape hole
  discrf2d_type round;  // distribution for round shape hole
  discrf1d_type cr_att; // exponential to correct for crystal attenuation when do_depth

} pcf_type;

//.....structure for voxel information

typedef struct
{
  int ix; // column index
  int iy; // row index
  int iz; // slice index
  int ip; // inplane index (slice as an array)
  int iv; // volume index (considering the volume as an array) of the voxel

  float x; // x coordinate (cm, ref center of volume)
  float y; // y coordinate (cm, ref center of volume)
  float z; // z coordinate (cm, ref center of volume)

  float x1; // x coordinate in rotated framework (cm)
  float y1; // y coordinate in rotated framework (cm)

} voxel_type;

//.....structure for bin information

typedef struct
{
  float x; // x coordinate (cm, ref center of volume)
  float y; // y coordinate (cm, ref center of volume)
  float z; // z coordinate (cm, ref center of volume)

} bin_type;

//...... structure for LOR information

typedef struct
{
  //... all the following distances in cm .....................

  float x1d_l;  // x coordiante of intersection lor-detection plane in rotated reference system
  float z1d_l;  // z coordiante of intersection lor-detection plane in rotated reference system
  float x1dc_l; // x coordiante of intersection lor-detection plane + crystal in rotated reference system
  float z1dc_l; // z coordiante of intersection lor-detection plane + crystal in rotated reference system

  float hsxcm_d;    // size (cm) of the shadow of the hole at detector plane (x-azis)
  float hszcm_d;    // size (cm) of the shadow of the hole at detector plane (z-axis)
  float hsxcm_d_d2; // half of the size (cm) of the shadow of the hole at detector plane (x-azis)
  float hszcm_d_d2; // half of the size (cm) of the shadow of the hole at detector plane (z-axis)

  float hsxcm_dc;    // size (cm) of the shadow of the hole at detector + crystal plane (x-axis)
  float hszcm_dc;    // size (cm) of the shadow of the hole at detector + crystal plane (z-axis)
  float hsxcm_dc_d2; // half of the size (cm) of the shadow of the hole at detector plane (x-azis)
  float hszcm_dc_d2; // half of the size (cm) of the shadow of the hole at detector plane (z-axis)

  //... others .......................................................

  float eff; // effectiveness

} lor_type; // voxel-hole link

//::: functions :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

//... functions from wmtools_SPECT.cpp .........................................

void wm_alloc(const int* Nitems, wm_da_type& wm, const wmh_mph_type& wmh); // allocate wm

void read_prj_params_mph(wmh_mph_type& wmh); // read ring parameters from a file

void read_coll_params_mph(wmh_mph_type& wmh); // read collimator parameters from a file

// void which_hole();

void fill_pcf(const wmh_mph_type& wmh, pcf_type& pcf); // fill precalculated functions

void calc_cumsum(discrf2d_type* f);

void generate_msk_mph(bool* msk_3d,
                      const float* att,
                      const wmh_mph_type& wmh); // create a boolean mask for wm (no weights outside the msk)

void error_wmtools_SPECT_mph(int nerr, int ip, std::string txt); // error messages in wm_SPECT

} // namespace SPECTUB_mph

#endif //_WM_SPECT_H