Multipexing-hkem: allows the use of multiple "anatomical" images

examples/samples/KOSMAPOSL_hkem.par

@@ -57,18 +57,13 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
 
 input file := my_noisy_prompts_1ax_pos.hs
 
-; Daniel: here we have the possibility to choose the parameters which define the kernel matrix and the name of the MR image.
+; here we have the possibility to choose the parameters which define the kernel matrix and the name of the MR image.
 
 projector pair type := Matrix
 Projector Pair Using Matrix Parameters :=
 Matrix type := Ray Tracing
 Ray tracing matrix parameters :=
 number of rays in tangential direction to trace for each bin := 10
-do symmetry 90degrees min phi := 1
-do_symmetry_180degrees_min_phi:= 1
-do_symmetry_swap_s:= 1
-do_symmetry_swap_segment:= 1
-do_symmetry_shift_z:= 1
 End Ray tracing matrix parameters :=
 End Projector Pair Using Matrix Parameters :=
 

examples/samples/KOSMAPOSL_mhkem.par

@@ -0,0 +1,102 @@
+KOSMAPOSLParameters :=
+; Example file for using [Multiplexing] Hybrid Kernelized Expectation Maximisation (HKEM or KEM).
+; The difference with the normal KOSMAPOSL is that we use multiple anatomical images and sigma_m
+; for more info.
+
+;the following disable the alpha coefficient output:
+disable output :=1
+
+
+; here we have the possibility to choose the parameters which define the kernel
+; matrix and the name of the anatomical image. the following are the defaults values:
+
+; 1 (default): use hybrid kernel (prior from MR and PET estimate)
+; OR
+; 0 kernel is MR-only
+hybrid:=1
+
+; Gaussian scaling parameter for the anatomical prior (units of image intensity)
+; It controls the edge preservation from the anatomical image, the bigger the stronger
+; default: 1
+sigma_m:= {0.1,0.1}
+
+
+; Gaussian scaling parameter for the PET estimate (units of PET image intensity)
+; It controls the edge preservation from the functional image, the bigger the stronger
+; default: 1
+sigma_p:=1
+
+
+; NB: sigma_dm and sigma_dp should be the same 
+; Spatial Gaussian scaling parameter for the anatomical prior (mm)
+; default: 1 (usual range 1-5)
+sigma_dm:=3
+
+; Spatial Gaussian scaling parameter for the PET prior (mm)
+; default: 1 (usual range 1-5)
+sigma_dp:=3	
+
+; Number of neigbouring voxels to compare: (num_neighbours X num_neighbours X num_neighbours)
+; default: 3
+number of neighbours:= 3
+
+; Number of non-zero elements in the feature vectors
+; default: 1
+number of non-zero feature elements:=1
+
+; this is the file name of the anatomical image
+anatomical image filenames:= {MRbrain.hv,CTbrain.hv}
+
+; the following should be 1 if you want to reconstruct 2D data
+only_2D:=1
+
+; the following is the output prefix of the PET reconstructed image
+kernelised output filename prefix :=MKOSMAPOSL
+
+objective function type:= PoissonLogLikelihoodWithLinearModelForMeanAndProjData
+PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
+
+input file := my_noisy_prompts_1ax_pos.hs
+
+; here we have the possibility to choose the parameters which define the kernel matrix and the name of the MR image.
+
+projector pair type := Matrix
+Projector Pair Using Matrix Parameters :=
+Matrix type := Ray Tracing
+Ray tracing matrix parameters :=
+number of rays in tangential direction to trace for each bin := 10
+End Ray tracing matrix parameters :=
+End Projector Pair Using Matrix Parameters :=
+
+; additive projection data to handle randoms etc
+additive sinogram := my_additive_sinogram_1ax_pos.hs 
+
+; norm and acf
+Bin Normalisation type := From ProjData
+Bin Normalisation From ProjData :=
+normalisation projdata filename:= my_multfactors_1ax_pos.hs
+End Bin Normalisation From ProjData:=
+
+; if the next parameters are enabled, 
+; the sensitivity will be computed and saved
+; use %d where you want the subset-number (a la printf)
+; we do this here for illustration, but also for re-use later on (to save some time)
+; CAREFUL: use correct number of subsets in name to avoid confusion
+
+subset sensitivity filenames:= sens_%d.hv
+recompute sensitivity := 1
+
+xy output image size (in pixels) := 256
+end PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
+
+
+; if you want to continue from a particular image
+; initial estimate:= reconML.hv
+
+; Number of subsets should usually be a divisor of num_views/8
+; the following is an example for the Siemens mMR
+number of subsets:= 21
+number of subiterations:= 630
+save estimates at subiteration intervals:= 21
+
+END KOSMAPOSLParameters:=

src/include/stir/KOSMAPOSL/KOSMAPOSLReconstruction.h

@@ -52,7 +52,10 @@ class PoissonLogLikelihoodWithLinearModelForMean;
   This class implements the iterative algorithm obtained using the Kernel method (KEM) and Hybrid kernel method (HKEM).
   This implementation corresponds to the one presented by Deidda D et al, ``Hybrid PET-MR list-mode kernelized expectation maximization reconstruction",
   Inverse Problems, 2019, DOI: https://doi.org/10.1088/1361-6420/ab013f. However, this allows
-  also sinogram-based reconstruction. Each voxel value of the image, \f$ \boldsymbol{\lambda}\f$, can be represented as a
+  also sinogram-based reconstruction. Furtermore, it is possible to use multiple side information images (they can be anatomical images or images from
+  different modalities). This extension is called multiplexing-HKEM and is described in Deidda et al. ``Multiplexing Kernelized Expectation Maximization
+  Reconstruction for PET-MR." IEEE NSS/MIC Proceedings (NSS/MIC) 2018, DOI: https://doi.org/10.1109/NSSMIC.2018.8824312.
+  Each voxel value of the image, \f$ \boldsymbol{\lambda}\f$, can be represented as a
   linear combination using the kernel method.  If we have an image with prior information, we can construct for each voxel
   \f$ j \f$ of the emission image a feature vector, $\f \boldsymbol{v}_j \f$, using the prior information. The voxel value,
   \f$\lambda_j\f$, can then be described using the kernel matrix
@@ -149,30 +152,38 @@ class KOSMAPOSLReconstruction:
   const KOSMAPOSLReconstruction& get_parameters() const
     {return *this;}
 
+  //@{
   //kernel
-  const std::string get_anatomical_image_filenames() const;
+  const std::vector<std::string> get_anatomical_image_filenames() const;
   const int get_num_neighbours() const;
   const int get_num_non_zero_feat() const;
-  const double get_sigma_m() const;
+  const std::vector<double> get_sigma_m() const;
   const double get_sigma_p() const;
   const double get_sigma_dp() const;
   const double get_sigma_dm()const;
   const bool get_only_2D() const;
   const bool get_hybrid()const;
 
-  shared_ptr<const TargetT> get_anatomical_prior_sptr() const;
+  std::vector<shared_ptr<TargetT> > get_anatomical_prior_sptr();
+//@}
 
     /*! \name Functions to set parameters
     This can be used as alternative to the parsing mechanism.
    \warning Be careful with setting shared pointers. If you modify the objects in 
    one place, all objects that use the shared pointer will be affected.
   */
+ //@{
+  void set_anatomical_prior_sptr(shared_ptr<TargetT>, int index);
+  //! sets all elements of vector anatomical_prior to the same value
+  void set_anatomical_prior_sptr(shared_ptr<TargetT> arg);
+  void set_anatomical_image_filenames(const std::string&, const int index);
+  void set_anatomical_image_filenames(const std::string&);
 
-  void set_anatomical_prior_sptr(shared_ptr<TargetT>);
 
-  void set_anatomical_image_filenames(const std::string&);
   void set_num_neighbours(const int);
   void set_num_non_zero_feat(const int);
+  void set_sigma_m(const double, const int index);
+  //! sets all elements of vector sigma_m to the same value
   void set_sigma_m(const double);
   void set_sigma_p(const double);
   void set_sigma_dp(const double);
@@ -182,7 +193,7 @@ class KOSMAPOSLReconstruction:
 
   //! boolean value to determine if the update images have to be written to disk
   void set_write_update_image(const int);
-
+  //@}
 
   //! prompts the user to enter parameter values manually
   virtual void ask_parameters();
@@ -196,11 +207,11 @@ class KOSMAPOSLReconstruction:
   //! Anatomical image filename
   std::vector<std::string> anatomical_image_filenames;
 
-  shared_ptr<TargetT> anatomical_prior_sptr;
-  shared_ptr<TargetT> kpnorm_sptr,kmnorm_sptr;
+  std::vector<shared_ptr<TargetT> > anatomical_prior_sptr,kmnorm_sptr;
+  shared_ptr<TargetT> kpnorm_sptr;
  //kernel parameters
   int num_neighbours,num_non_zero_feat,num_elem_neighbourhood,num_voxels,dimz,dimy,dimx;
-  double sigma_m;
+  std::vector<double> sigma_m;
   bool only_2D;
   bool hybrid;
   double sigma_p;
@@ -224,7 +235,7 @@ class KOSMAPOSLReconstruction:
   virtual void update_estimate (TargetT& current_image_estimate);
 
   int subiteration_counter;
-  double anatomical_sd;
+  std::vector<double> anatomical_sd;
   mutable Array<3,float> distance;
   /*! Create a matrix containing the norm of the difference between two feature vectors, \f$ \|  \boldsymbol{z}^{(n)}_j-\boldsymbol{z}^{(n)}_l \| \f$. */
   /*! This is done for the emission image which keeps changing*/
@@ -235,12 +246,12 @@ class KOSMAPOSLReconstruction:
 
   /*! Create a matrix similarly to calculate_norm_matrix() but this is done for the anatomical image, */
   /*! which does not  change over iteration.*/
-    void calculate_norm_const_matrix(TargetT &normm,
+    void calculate_norm_const_matrix(std::vector<shared_ptr<TargetT> > &normm,
                                 const int dimf_row,
                                 const int dimf_col);
 
   /*! Estimate the SD of the anatomical image to be used as normalisation for the feature vector */
-    double estimate_stand_dev_for_anatomical_image();
+    void estimate_stand_dev_for_anatomical_image(std::vector<double> &SD);
 
   /*! Compute for each voxel, jl, of the emission image the linear combination between the coefficient \f$ \alpha_{jl} \f$ and the kernel matrix \f$ k_{jl} \f$\f$ */
   /*! The information is stored in the image, kImage */
@@ -271,7 +282,8 @@ class KOSMAPOSLReconstruction:
                                 const double distance_dzdydx,
                                 const bool use_compact_implementation,
                                 const int l,
-                                const int m);
+                                const int m,
+                                const int index);
 
   double calc_kernel_from_precalculated(const double precalculated_norm_zxy,
                                         const double sq_sigma_int,