SPECTUB_Tools.h

/*
    Copyright (c) 2013, Biomedical Image Group (GIB), Universitat de Barcelona, Barcelona, Spain.
    Copyright (c) 2013, University College London
    This file is part of STIR.
 
    SPDX-License-Identifier: Apache-2.0
 
    See STIR/LICENSE.txt for details
 
  \author Carles Falcon
*/
 
#ifndef _WM_SPECTUB_H
#define _WM_SPECTUB_H
 
#include <string>
#include <vector>
 
namespace SPECTUB
{
 
//::: srtuctures ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

typedef struct
{
  int num;    // number of collimator (see weight_64.cpp for options)
  bool do_fb; // true: fanbeam collimator || false: parallel collimator
 
  //... parallel collimator parameters .....................
 
  float A; // linear factor for dependency of sigma on distance: sigma=A*dist+B (parallel, fanbeam-vertical)
  float B; // independent factor for dependency of sigma on distance: sigma=A*dist+B (parallel, fanbeam-vertical)
 
  //... fanbeam collimator parameters ......................
 
  float F;     // Focal length (fanbeam)
  float L;     // collimator to detector distance (?) (fanbeam)
  float A_h;   // linear factor for dependency of sigma on distance (fanbeam horizontal)
  float A_v;   // linear factor for dependency of sigma on distance (fanbeam horizontal)
  float D;     // collimator thicknes?? (fanbeam)
  float w;     // collimator thickness (?) (fanbeam)
  float insgm; // intrinsic sigma (cristal resolution) (fanbeam)
 
} collim_type;

typedef struct
{
  int Nrow; // number of rows
  int Ncol; // number of columns
  int Nsli; // number of slices
  int Npix; // number of pixels (axial planes)
  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 Nrowd2; // half of Nrow
  float Ncold2; // half of Ncol
  float Nslid2; // half of Nsli
 
  float Xcmd2; // Half of the size of the volume, dimension x (cm);
  float Ycmd2; // Half of the size of the volume, dimension y (cm);
  float Zcmd2; // 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;

typedef struct
{
  int Nbin;    // length of the detection line in bins (number of bins per line)
  float lngcm; // length of the detection line in cm.
  float szcm;  // bin size in cm
 
  int Nsli;   // number of slices
  float thcm; // slice thickness in cm
 
  int Nang;                  // number of projection angles
  std::vector<float> angles; // projection angles. degrees from upper detection plane (parallel to table). Negative for CW
                             // rotations (see manual)
 
  int NOS;    // number of subsets in which to split the matrix
  int NangOS; // Number of angles in each subset = Nang/NOS
  int Nbp;    // number of bins in a 2D projection=lng*Nsli
  int Nbt;    // total number of bins= Nbp*Nang
  int NbOS;   // total number of bins per subset= Nbp*NangOS = Nbt/NOS
  int* order; // order of the angles of projection (array formed by indexs of angles belonging to consecutive subsets)
 
  float Nbind2;  // length of the detection line (in bins) divided by 2
  float lngcmd2; // length of the detection line in cm divided by 2
  float Nslid2;  // number of slices divided by 2
 
} proj_type;

typedef struct
{
  int subset_ind; // subset index of this matrix (-1: all subsets in a single file)
  int* index;     // included angles into this subset (index. Multiply by increm to get corresponding angles in degrees)
 
  float* Rrad;   // Rotation radius (one value for each projection angle)
  float min_w;   // minimum weight to be taken into account
  float psfres;  // spatial resolution of continous distributions in PSF calculation
  float maxsigm; // maximum number of sigmas in PSF calculation
 
  bool fixed_Rrad;  // true: fixed radius of projection || false: variable radius of projection
  bool do_psf;      // true: to correct for PSF         || false: do not correct for PSF
  bool do_psf_3d;   // true: 3d correction for PSF      || false: 2d correction for PSF
  bool predef_col;  // true: predefined collimator      || false: user defined PSF parametres
  bool do_att;      // true: to correct for attenuation || false: do not correct for attenuation
  bool do_full_att; // true: diff att for each PSF bin  || false: the whole PSF has the same att. factor (central line)
  bool do_msk;      // true: weights just inside msk    || false: weights for the whole FOV
  bool do_msk_slc;  // true: weights for several slices || false: weights for all slices
  bool do_msk_cyl;  // true: to use cylinder as a mask  || false: not to use cylinder as a mask
  bool do_msk_att;  // true: to use att map as a mask   || false: not to use att map as a mask
  bool do_msk_file; // true: explicit mask              || false: not to use explicit mask
 
  std::string att_fn;  // attenuation map filename
  std::string msk_fn;  // explicit mask filename
  std::string col_fn;  // collimator parameters filename
  std::string Rrad_fn; // rotation radius file name
 
  volume_type vol; //
  proj_type prj;   //
  collim_type COL; // collimator structure (see weight3d_64b.cpp for options)
 
} wmh_type;

typedef struct
{
  // weight matrix dimensions
 
  int ne;   // nonzero elements
  int NbOS; // dimension 1 (rows) of the weight matrix (NbOS or NBt)
  int Nvox; // dimension 2 (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
 
  bool do_save_wmh; // to save or not to save weight_mat header info into weight_mat file
 
} wm_type;

typedef struct
{
  int NbOS;    // dimension 1 (rows) of the weight matrix (NbOS or NBt)
  int Nvox;    // dimension 2 (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 base name (filename without extension index)
  std::string OSfn;   // matrix filename
  std::string fn_hdr; // matrix header file name
 
  //... 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)
 
  //... format ...............................................
 
  bool do_save_wmh;  // to save or not to save weight_mat header info into weight_mat file
  bool do_save_STIR; // to save weight matrix with STIR format
 
} wm_da_type;

typedef struct
{
  int lng;    // length (in discretization intervals) (odd number)
  int lngd2;  // half of the length (in discretization intervals) (lng-1)/2
  float res;  // spatial resolution of distfunc (discretization interval)
  float* val; // array of values
  float* acu; // distribution function values (cumulative sum)
 
} discrf_type;

typedef struct
{
  int maxszb; // maximum size in bins (for allocation purposes)
 
  int di;   // discretization interval (to reduce spatial resolution to bin resolution). (int: #points)
  int* ind; // projection indexs for the bins of the PSF (horizontal)
  int Nib;  // actual number of bins forming the PSF (length of PSF in bins)
 
  float sgmcm;   // sigma of the PSF in cm
  float lngcm;   // length of PSF (in cm)
  float lngcmd2; // half of the length of PSF (in cm)
  float* val;    // array of values
  float efres;   // effective resolution (psfres rescaled to real psf length)
 
} psf1d_type;

typedef struct
{
  int maxszb_h; // maximum size in bins horizontal (for allocation purposes)
  int maxszb_v; // maximum size in bins vertical (for allocation purposes)
  int maxszb_t; // maximum size in bins total (for allocation purposes)
 
  int* ib; // projection indexs for the bins of the PSF (horizontal)
  int* jb; // projection indexs for the bins of the PSF (vertical)
  int Nib; // actual number of bins forming the PSF (length of PSF in bins)
 
  float* val; // array of values
 
} psf2da_type;

typedef struct
{
  int ind;      // index of angle considering the whole set of projections (sequential order: 0->Nang-1)
  int indOS;    // index of angle considering the subjet
  int iOS_proj; // index of the first bin for this angle (in subset of projections)
 
  float cos; // coninus of the angle
  float sin; // sinus of the angle
 
  // parametres for describng the trapezoidal projection of a square voxel
 
  float p;           // plateau higness
  float m;           // slope of the trapezoid
  float n;           // independent term of the slope
  int N1;            // index of the first vertice (end of plateau) in DX units
  int N2;            // index of the second vertice (end of the slope) in DX units
  discrf_type vxprj; // projection of a square voxel in this direction (for no PSF)
 
  // variable rotation radius
 
  float Rrad; // rotation radius for this angle
 
  // first bin position and increments
 
  float xbin0; // x coordinate for the first bin of the detection line corresponding to this angle
  float ybin0; // y coordinate for the first bin of the detection line corresponding to this angle
  float incx;  // increment in x to the following bin in detection line
  float incy;  // increment in y to the following bin in detection line
 
} angle_type;

typedef struct
{
  float szcm; // voxel size (side length in cm)
  float thcm; // voxel thickness (cm)
 
  int irow; // row index
  int icol; // column index
  int islc; // slice index
  int ip;   // in plane index (considering the slice as an array) of the voxel
  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 coordinade in rotated framework
 
  float dv2dp;  // distance from voxel to detection plane
  float costhe; // cosinus of theta (angle between focal-voxel line and line perpendicular to detection plane) (fanbeam)
  float xdc;    // distance (cm over detection line) from projected voxel to the center of the detection line
  float xd0;    // distance (cm over detection line) from projected voxel to the begin of the detection line
  float zd0;    // distance (cm) to the lowest plane of the volume
 
} voxel_type;

typedef struct
{
  float szcm;   // bin size (cm)
  float szcmd2; // half of the above value
  float thcm;   // bin thickness (cm)
  float thcmd2; // bin thickness (cm)
 
  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 szdx; // bin size in resolution units
  float thdx; // bin thickness in resolution units
 
} bin_type;

typedef struct
{
  float* dl;  // distance of attenuation path on each crossed voxel of the attenuation map
  int* iv;    // in-plane index (considering slices of attmap as an array) of any crossed voxel of the attenuation map
  int lng;    // number of elements in the attenuation path
  int maxlng; // maximum number of elements in the attenuation path (for allocation)
 
} attpth_type;
 
//::: functions :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 
//... functions from wmtools_SPECT.cpp .........................................
 
void index_calc(int* indexs, const wmh_type& wmh); // to calculate projection index order in subsets
 
//
// void col_params ( collim_type *COL );              // to fill collimator structure
 
void fill_ang(angle_type* ang, const wmh_type& wmh, const float* Rrad); // to fill angle structure
 
void generate_msk(bool* msk_3d,
                  bool* msk_2d,
                  const float* att,
                  const volume_type* vol,
                  const wmh_type& wmh); // to create a boolean mask for wm (no weights outside the msk)
 
int max_psf_szb(const angle_type* ang, const wmh_type& wmh);
 
float calc_sigma_h(voxel_type vox, collim_type COL);
 
float calc_sigma_v(voxel_type vox, collim_type COL);
 
void free_wm(wm_type* f); // to free weight_mat
 
void free_wm_da(wm_da_type* f); // to free weight_mat_da
 
} // namespace SPECTUB
 
#endif //_WM_SPECT_H