CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00002 /* 00003 Implementation of class iImageConvolverStationaryGaussian 00004 00005 - separators: done 00006 - doxygen: done 00007 - default initialization: done 00008 - CASTOR_DEBUG: 00009 - CASTOR_VERBOSE: 00010 */ 00011 00018 #include "vImageConvolver.hh" 00019 #include "oImageDimensionsAndQuantification.hh" 00020 #include "oImageSpace.hh" 00021 #include "iImageConvolverStationaryGaussian.hh" 00022 00023 // ===================================================================== 00024 // --------------------------------------------------------------------- 00025 // --------------------------------------------------------------------- 00026 // ===================================================================== 00027 00028 iImageConvolverStationaryGaussian::iImageConvolverStationaryGaussian() : vImageConvolver() 00029 { 00030 m_nbSigmas = -1.; 00031 m_transFWHM = -1.; 00032 m_axialFWHM = -1.; 00033 m_dimKernelXY = -1; 00034 m_dimKernelXYZ = -1; 00035 } 00036 00037 // ===================================================================== 00038 // --------------------------------------------------------------------- 00039 // --------------------------------------------------------------------- 00040 // ===================================================================== 00041 00042 iImageConvolverStationaryGaussian::~iImageConvolverStationaryGaussian() 00043 { 00044 } 00045 00046 // ===================================================================== 00047 // --------------------------------------------------------------------- 00048 // --------------------------------------------------------------------- 00049 // ===================================================================== 00050 00051 void iImageConvolverStationaryGaussian::ShowHelp() 00052 { 00053 cout << "This convolver is based on a classic stationary Gaussian kernel." << endl; 00054 cout << "It can be anisotropic so that the transaxial and axial FWHM can be different." << endl; 00055 cout << "One also have to choose the number of sigmas of the Gaussian distribution that will be included in the kernel." << endl; 00056 cout << "The following options can be used (in this particular order when provided as a list, separated by commas):" << endl; 00057 cout << " transaxial FWHM: to set the transaxial FWHM (in mm)" << endl; 00058 cout << " axial FWHM: to set the axial FWHM (in mm)" << endl; 00059 cout << " number of sigmas: to set the number of sigmas included in the kernel (recommendations: at least 3. and maximum 5.)" << endl; 00060 } 00061 00062 // ===================================================================== 00063 // --------------------------------------------------------------------- 00064 // --------------------------------------------------------------------- 00065 // ===================================================================== 00066 00067 int iImageConvolverStationaryGaussian::ReadConfigurationFile(const string& a_configurationFile) 00068 { 00069 string key_word = ""; 00070 // Read the transaxial FWHM option 00071 key_word = "transaxial FWHM"; 00072 if (ReadDataASCIIFile(a_configurationFile, key_word, &m_transFWHM, 1, KEYWORD_MANDATORY)) 00073 { 00074 Cerr("***** iImageConvolverStationaryGaussian::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl); 00075 return 1; 00076 } 00077 // Read the axial FWHM option 00078 key_word = "axial FWHM"; 00079 if (ReadDataASCIIFile(a_configurationFile, key_word, &m_axialFWHM, 1, KEYWORD_MANDATORY)) 00080 { 00081 Cerr("***** iImageConvolverStationaryGaussian::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl); 00082 return 1; 00083 } 00084 // Read the number of sigma option 00085 key_word = "number of sigmas"; 00086 if (ReadDataASCIIFile(a_configurationFile, key_word, &m_nbSigmas, 1, KEYWORD_MANDATORY)) 00087 { 00088 Cerr("***** iImageConvolverStationaryGaussian::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl); 00089 return 1; 00090 } 00091 // Normal end 00092 return 0; 00093 } 00094 00095 // ===================================================================== 00096 // --------------------------------------------------------------------- 00097 // --------------------------------------------------------------------- 00098 // ===================================================================== 00099 00100 int iImageConvolverStationaryGaussian::ReadOptionsList(const string& a_optionsList) 00101 { 00102 // There are 3 floating point variables as options 00103 FLTNB options[3]; 00104 // Read them 00105 if (ReadStringOption(a_optionsList, options, 3, ",", "Stationary Gaussian convolver configuration")) 00106 { 00107 Cerr("***** iImageConvolverStationaryGaussian::ReadConfigurationFile() -> Failed to correctly read the list of options !" << endl); 00108 return 1; 00109 } 00110 // Affect options 00111 m_transFWHM = options[0]; 00112 m_axialFWHM = options[1]; 00113 m_nbSigmas = options[2]; 00114 // Normal end 00115 return 0; 00116 } 00117 00118 // ===================================================================== 00119 // --------------------------------------------------------------------- 00120 // --------------------------------------------------------------------- 00121 // ===================================================================== 00122 00123 int iImageConvolverStationaryGaussian::CheckSpecificParameters() 00124 { 00125 // Check number of sigmas 00126 if (m_nbSigmas<=0.) 00127 { 00128 Cerr("***** iImageConvolverStationaryGaussian::CheckSpecificParameters() -> Number of sigmas included in the kernel must be strictly positive !" << endl); 00129 return 1; 00130 } 00131 // Check the transaxial FWHM 00132 if (m_transFWHM<0.) 00133 { 00134 Cerr("***** iImageConvolverStationaryGaussian::CheckSpecificParameters() -> Transaxial FWHM is negative !" << endl); 00135 return 1; 00136 } 00137 // Check the axial FWHM 00138 if (m_axialFWHM<0.) 00139 { 00140 Cerr("***** iImageConvolverStationaryGaussian::CheckSpecificParameters() -> Axial FWHM is negative !" << endl); 00141 return 1; 00142 } 00143 // Normal end 00144 return 0; 00145 } 00146 00147 // ===================================================================== 00148 // --------------------------------------------------------------------- 00149 // --------------------------------------------------------------------- 00150 // ===================================================================== 00151 00152 int iImageConvolverStationaryGaussian::BuildConvolutionKernel() 00153 { 00154 // Verbose 00155 if (m_verbose>=2) 00156 { 00157 Cout("iImageConvolverStationaryGaussian::BuildConvolutionKernel() -> Compute convolution kernel" << endl); 00158 Cout(" --> Transaxial FWHM: " << m_transFWHM << " mm" << endl); 00159 Cout(" --> Axial FWHM: " << m_axialFWHM << " mm" << endl); 00160 Cout(" --> Included sigmas: " << m_nbSigmas << endl); 00161 } 00162 00163 // -------------------------------- 00164 // Stationary kernel 00165 // -------------------------------- 00166 00167 // Set the number of kernels to 1 00168 m_nbKernels = 1; 00169 // Specify that it is stationary (we do not overload the Convolve() function here) 00170 m_stationary = true; 00171 // Compute kernel half size based on the number of sigmas included in the kernel 00172 FLTNB two_sqrt_two_ln_two = 2.*sqrt(2.*log(2.)); 00173 FLTNB half_kern_floatX = m_transFWHM * m_nbSigmas / (two_sqrt_two_ln_two*mp_ImageDimensionsAndQuantification->GetVoxSizeX()); 00174 FLTNB half_kern_floatY = m_transFWHM * m_nbSigmas / (two_sqrt_two_ln_two*mp_ImageDimensionsAndQuantification->GetVoxSizeY()); 00175 FLTNB half_kern_floatZ = m_axialFWHM * m_nbSigmas / (two_sqrt_two_ln_two*mp_ImageDimensionsAndQuantification->GetVoxSizeZ()); 00176 // For each half size, if the decimal part is below 0.5, then it will be taken into account into the central kernel's voxel, 00177 // otherwise, we had a voxel to the half size 00178 INTNB half_kern_dimX = ((INTNB)(half_kern_floatX)); 00179 if (half_kern_floatX-((FLTNB)half_kern_dimX)>0.5) half_kern_dimX++; 00180 INTNB half_kern_dimY = ((INTNB)(half_kern_floatY)); 00181 if (half_kern_floatY-((FLTNB)half_kern_dimY)>0.5) half_kern_dimY++; 00182 INTNB half_kern_dimZ = ((INTNB)(half_kern_floatZ)); 00183 if (half_kern_floatZ-((FLTNB)half_kern_dimZ)>0.5) half_kern_dimZ++; 00184 // Allocate kernel dimensions 00185 mp_dimKernelX = (INTNB*)malloc(1*sizeof(INTNB)); 00186 mp_dimKernelY = (INTNB*)malloc(1*sizeof(INTNB)); 00187 mp_dimKernelZ = (INTNB*)malloc(1*sizeof(INTNB)); 00188 // Compute kernel dimensions 00189 mp_dimKernelX[0] = half_kern_dimX * 2 + 1; 00190 mp_dimKernelY[0] = half_kern_dimY * 2 + 1; 00191 mp_dimKernelZ[0] = half_kern_dimZ * 2 + 1; 00192 m_dimKernelXY = mp_dimKernelX[0] * mp_dimKernelY[0]; 00193 m_dimKernelXYZ = m_dimKernelXY * mp_dimKernelZ[0]; 00194 // Verbose 00195 if (m_verbose>=2) Cout (" --> Kernel dimensions: [" << mp_dimKernelX[0] << ";" << mp_dimKernelY[0] << ";" << mp_dimKernelZ[0] << "]" << endl); 00196 // Allocate kernel 00197 m2p_kernel = (FLTNB**)malloc(1*sizeof(FLTNB*)); 00198 m2p_kernel[0] = (FLTNB*)malloc(m_dimKernelXYZ*sizeof(FLTNB)); 00199 00200 // -------------------------------- 00201 // Compute the kernel 00202 // -------------------------------- 00203 00204 // Variables 00205 FLTNB sigma_x = m_transFWHM / two_sqrt_two_ln_two; 00206 FLTNB sigma_y = m_transFWHM / two_sqrt_two_ln_two; 00207 FLTNB sigma_z = m_axialFWHM / two_sqrt_two_ln_two; 00208 FLTNB mu_x = (FLTNB)(half_kern_dimX); 00209 FLTNB mu_y = (FLTNB)(half_kern_dimY); 00210 FLTNB mu_z = (FLTNB)(half_kern_dimZ); 00211 // Compute kernel by integration of the Gaussian distribution (this method is not exactly exact... but is closer to it compared to single Gaussian points) 00212 FLTNB sum_kernel = 0.; 00213 for (int z=0, index=0; z<mp_dimKernelZ[0]; z++) 00214 { 00215 FLTNB kern_Z = 0.5 * ( erf( ( (((FLTNB)z)-mu_z+0.5)*mp_ImageDimensionsAndQuantification->GetVoxSizeZ()) / (sqrt(2.)*sigma_z) ) 00216 - erf( ( (((FLTNB)z)-mu_z-0.5)*mp_ImageDimensionsAndQuantification->GetVoxSizeZ()) / (sqrt(2.)*sigma_z) ) ); 00217 for (int y=0; y<mp_dimKernelY[0]; y++) 00218 { 00219 FLTNB kern_ZY = 0.5 * ( erf( ( (((FLTNB)y)-mu_y+0.5)*mp_ImageDimensionsAndQuantification->GetVoxSizeY()) / (sqrt(2.)*sigma_y) ) 00220 - erf( ( (((FLTNB)y)-mu_y-0.5)*mp_ImageDimensionsAndQuantification->GetVoxSizeY()) / (sqrt(2.)*sigma_y) ) ) 00221 * kern_Z; 00222 for (int x=0; x<mp_dimKernelX[0]; x++, index++) 00223 { 00224 m2p_kernel[0][index] = 0.5 * ( erf( ( (((FLTNB)x)-mu_x+0.5)*mp_ImageDimensionsAndQuantification->GetVoxSizeX()) / (sqrt(2.)*sigma_x) ) 00225 - erf( ( (((FLTNB)x)-mu_x-0.5)*mp_ImageDimensionsAndQuantification->GetVoxSizeX()) / (sqrt(2.)*sigma_x) ) ) 00226 * kern_ZY; 00227 sum_kernel += m2p_kernel[0][index]; 00228 } 00229 } 00230 } 00231 // Normalize kernel 00232 for (int v=0; v<m_dimKernelXYZ; v++) m2p_kernel[0][v] /= sum_kernel; 00233 // Normal end 00234 return 0; 00235 } 00236 00237 // ===================================================================== 00238 // --------------------------------------------------------------------- 00239 // --------------------------------------------------------------------- 00240 // =====================================================================