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

examples/samples/KOSMAPOSL_mhkem.par

@@ -0,0 +1,107 @@
+KOSMAPOSLParameters :=
+; Example file for using [Hybrid] Kernelized Expectation Maximisation (HKEM or KEM).
+; See documentation of KOSMAPOSLreconstruction
+; 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
+
+; Daniel: 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 :=
+
+; 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

@@ -150,30 +150,30 @@ class KOSMAPOSLReconstruction:
     {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>);
+  void set_anatomical_prior_sptr(shared_ptr<TargetT>&, int index);
+  void set_anatomical_image_filenames(std::string&, int index);
 
-  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);
+  void set_sigma_m(double&, int &index);
   void set_sigma_p(const double);
   void set_sigma_dp(const double);
   void set_sigma_dm(const double);
@@ -196,11 +196,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 +224,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 +235,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 +271,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,