#ifndef ALGOSTIROSMAPOSL_H #define ALGOSTIROSMAPOSL_H /** @file AlgoStirOSMAPOSL.h @brief Class which is interfacing the OSMAPOST algorithm of STIR with the TIRIUS user interface. @author Jean-Daniel Leroux, @version 1.1 @date June 2008 */ /* -------------------------------------------------------------------------- Copyright (C) 2006, 2008, Jean-Daniel Leroux Copyright (C) 2006, 2008, Universite de Sherbrooke This file is part of the TIRIUS project. TIRIUS is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or any later version. TIRIUS is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with TIRIUS. If not, see . ----------------------------------------------------------------------------- */ #include "Algo.h" #include "AlgoGlobal.h" #include "ErrorHandler.h" #include "InterfileConverter.h" #include "stir/OSMAPOSL/OSMAPOSLReconstruction.h" #include "stir/Succeeded.h" #include "stir/LogLikBased/common.h" #include "stir/recon_buildblock/distributable.h" #include "stir/recon_buildblock/ProjectorByBinPairUsingProjMatrixByBin.h" #include "stir/ProjData.h" #include "stir/VoxelsOnCartesianGrid.h" #include "stir/utilities.h" #define CURRENT_PROJECTION_DATA_FILENAME "stir_data/currentSinogram" #define CURRENT_NORM_FACTOR_DATA_FILENAME "stir_data/currentSinoNormFactor" #define CURRENT_ATT_FACTOR_DATA_FILENAME "stir_data/currentSinoAttFactor" #define CURRENT_SENSIBILITY_DATA_FILENAME "stir_data/currentImageSensitivity" #define CURRENT_IMAGE_OUTPUT_DATA_FILENAME "stir_data/currentImageOutput" #define DEFAULT_ALGO_OSMAPOSL_MATRIX_ID MATRIX_MLEM_ID #define NB_OSMAPOSL_VERSION 5 ///< Number of active version of the OSMAPOSL algorithm #define STR_OSMAPOSL_VERSION1 "Matrix" ///< Name of the first OSMAPOSL algorithm version #define STR_OSMAPOSL_VERSION2 "Ray tracing" ///< Name of the second OSMAPOSL algorithm version #define STR_OSMAPOSL_VERSION3 "Matrix + post-filter" ///< Name of the thirth OSMAPOSL algorithm version #define STR_OSMAPOSL_VERSION4 "Matrix + Q-prior" ///< Name of the fourth OSMAPOSL algorithm version #define STR_OSMAPOSL_VERSION5 "Matrix + MR-prior" ///< Name of the fifth OSMAPOSL algorithm version #define OSMAPOSL_PARAM_FILENAME_MAX_SIZE 200 #define DEFAULT_OSMAPOSL_PARAM_FILENAME_V1 "stir_parameters/Default_OSMAPOSL_Parameters_v1.par" #define DEFAULT_OSMAPOSL_PARAM_FILENAME_V2 "stir_parameters/Default_OSMAPOSL_Parameters_v2.par" #define DEFAULT_OSMAPOSL_PARAM_FILENAME_V3 "stir_parameters/Default_OSMAPOSL_Parameters_v3.par" #define DEFAULT_OSMAPOSL_PARAM_FILENAME_V4 "stir_parameters/Default_OSMAPOSL_Parameters_v4.par" #define DEFAULT_OSMAPOSL_PARAM_FILENAME_V5 "stir_parameters/Default_OSMAPOSL_Parameters_v5.par" // enMLEMVersion typedef enum {OSMAPOSL_version1, OSMAPOSL_version2, OSMAPOSL_version3, OSMAPOSL_version4, OSMAPOSL_version5} enOSMAPOSLVersion; ///< Enumeration to select the algorithm version we want to use. USING_NAMESPACE_STIR //! Interface to set or access information about the OSMAPOSL reconstruction process. /** The Ordered Subset MAP One Step Late (OSMAPOSL) algorithm comes from the Software for Tomographic * Image Reconstruction (STIR) program. Accordingly, the user should read the STIR documentation for * knowing exacly the available options and how to set the algorithm parameters. * * This class act like a wrapper in order to connect the TIRIUS user-interface with the OSMAPOSL algorithm. * All the functionalities and parameters of the OSMAPOSL algorithm are controlled by deriving this class * from the OSMAPOSLReconstruction. We find it easier to do it this way but we know that it could have been * cleaner to simply created a OSMAPOSLReconstruction object... but we were facing some problems. * * As you will see, some parameters of the OSMAPOSL algorithm can be modified from the TIRIUS user-interface. * Nevertheless, since there are many options that can be set in OSMAPOSL, not all of them can be modified from * TIRIUS. For more complex parameter setting like selecting filter or prior one can load a parameter file (.par) * in interfile format which contain all the desired settings. In fact, there is two way to do so. One is to modified * one of the parameter file associated with one of the OSMAPOSL version (see the DEFAULT_OSMAPOSL_PARAM_FILENAME * definitions). The other way to proceed is to load directly a parameter file using SetAlgoOptions() function * associated to the first algorithm parameter. * * Another thing to know about the OSMAPOSL algorithm is that it is using a system matrix (or if you prefer projectors) * which are implemented direclty in the STIR software. Accordingly, you do not need to load a system matrix * before performing an image reconstruction with this algorithm. Moreover, since this algorithm may be based on * system matrix (or projector) which are computed on the flag, they may be more suitable to solve bigger tomographic * problems where the available RAM memory may be restrictive. * */ class cAlgoStirOSMAPOSL : public OSMAPOSLReconstruction, public cAlgo { public: cAlgoStirOSMAPOSL(); cAlgoStirOSMAPOSL(mINT matrixID); ~cAlgoStirOSMAPOSL(); virtual mINT DeleteAll(); virtual mINT SetAlgoOptions(float Value, int noParam); virtual mINT SetAlgoOptions(int Sel, int noParam); virtual mINT SetAlgoOptions(char *FileName, int noParam); virtual mINT GetCopyOfImage(float *CopyImage, float *minVal, float *maxVal, int selPlane, bool bUse3DContrast); mINT ReconstructNewImage(char *AlgoParamFileName); virtual mINT ReconstructNewImage(stSinogram *NewSinogram, stMatrixPtr MatPtr, stMatrixParams MatParams); virtual mINT ReconstructImage(); private: enOSMAPOSLVersion mAlgoVersion; ///< Variable used to select the version of the algorithm we want to use to reconstruct the image // Parameters that can be modified by the user throught the interface int mNbSubset; bool mbRandomisedSubsetOrder; ///< Flag to determine if the subset order should be randomized (see STIR documentation). bool mbEnforceInitialPositivity; ///< Flag to determine if positivity should be enforce on the image estimate (see STIR documentation). bool mbZeroSegEndPlane; ///< Flag to determine if the end segment plane are zeroed (see STIR documentation). bool mbDoRimTruncation; ///< Flag to determine if rim truncation is performed (see STIR documentation). float mFOVDiameter; ///< Useful field of view (FOV) diameter in millimeter. int mMaxRingDifference; ///< Maximum ring difference for which the coincidence are accept in the axial direction int mNbPixelXY; ///< Number of pixels in the X-axis and Y-axis direction (tangential) int mNbPixelZ; ///< Number of pixels in the Z-axis (axial) direction. float mPixelSizeXY; ///< Size of the pixel (in millimeter) in the X-axis and Y-axis direction (tangential). float mPixelSizeZ; ///< Size of the pixels in the Z-axis (axial) direction. float mImageZoom; ///< Zoom factor which is used to rescale the pixel size (X-axis and Y-axis) according to the default bin size (mBinSize) int mNbSubIterations; ///< Number of sub-iteration to be used for the image reconstruction. When set to zero we use the number of iteration set throught the TIRIUS user interface. // Variables that are required to modified some user selectable parameters of the class int mNbRing; ///< Number ring position in the apparatus. float mRingDiameter; ///< Scanner ring diameter (in millimeter). float mFOVHeight; ///< Useful field of view (FOV) height in millimeter. float mBinSize; ///< Width (in millimeter) of the parallel bin of the apparatus. mINT mCurrentScannerTypeID; char mCurrentOSMAPOSLParamFileName[OSMAPOSL_PARAM_FILENAME_MAX_SIZE]; // Private Functions //------------------------------------------ mINT InitDefaultParams(); mINT SetDefaultAlgoOptions(); mINT UpdateAlgoParams(); mINT WriteUserAlgoParamsSelection(); mINT ComputeSensitivityImage(stSinogram *NewSinogram); void compute_sensitivity_image(DiscretisedDensity<3,float>& result, shared_ptr const& proj_data_ptr, shared_ptr > const& attenuation_image_ptr, const bool do_attenuation, shared_ptr const& atten_proj_data_ptr, shared_ptr const& norm_proj_data_ptr); // Version 1 mINT ReconstructNewImageOSMAPOSL(stSinogram *NewSinogram, mINT nbIter, mINT nbLoop); mINT ReconstructImageOSMAPOSL(mINT nbIter); mINT SaveOneImageBin(char *FileName, float *Image, int nbPixelXY, int nbPixelZ); }; #endif