iOptimizerAPPGML.cc

Go to the documentation of this file.

#include "iOptimizerAPPGML.hh"
#include "sOutputManager.hh"

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

iOptimizerAPPGML::iOptimizerAPPGML() : vOptimizer()
{
  // ---------------------------
  // Mandatory member parameters
  // ---------------------------

  // Initial value at 1
  m_initialValue = 1.;
  // Only one backward image for APPGML
  m_nbBackwardImages = 1;
  // APPGML accepts penalties
  m_requiredPenaltyDerivativesOrder = 2;
  // APPGML is only compatible with histogram data
  m_listmodeCompatibility = false;
  m_histogramCompatibility = true;
  // APPGML is only compatible with emission data
  m_emissionCompatibility = true;
  m_transmissionCompatibility = false;
  // APPGML needs pre process loop to convolve region of interest once
  m_needPreIteration = true;

  // --------------------------
  // Specific member parameters
  // --------------------------

  m_dataSpaceDenominatorThreshold = -1.;
  m_minimumImageUpdateFactor = -1.;
  m_maximumImageUpdateFactor = -1.;
  m_bound = -1.;
  m_regionImagePath = "";
  m4p_imageNotShifted = NULL;
  m4p_firstDerivativePenaltyImage = NULL;
  m4p_secondDerivativePenaltyImage = NULL;
  mp_specificRegionOfInterest = NULL;
  mp_specificRegionOfInterestConvolved = NULL;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

iOptimizerAPPGML::~iOptimizerAPPGML()
{
  // Delete the image (not shifted)
  if (m4p_imageNotShifted)
  {
    // Loop over time basis functions
    for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
    {
      if (m4p_imageNotShifted[tbf])
      {
        // Loop over respiratory basis functions
        for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
        {
          if (m4p_imageNotShifted[tbf][rbf])
          {
            free(m4p_imageNotShifted[tbf][rbf]);
          }
        }
        free(m4p_imageNotShifted[tbf]);
      }
    }
    free(m4p_imageNotShifted);
  }
  
  // Delete the penalty image
  if (m4p_firstDerivativePenaltyImage)
  {
    // Loop over time basis functions
    for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
    {
      if (m4p_firstDerivativePenaltyImage[tbf])
      {
        // Loop over respiratory basis functions
        for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
        {
          if (m4p_firstDerivativePenaltyImage[tbf][rbf])
          {
            free(m4p_firstDerivativePenaltyImage[tbf][rbf]);
          }
        }
        free(m4p_firstDerivativePenaltyImage[tbf]);
      }
    }
    free(m4p_firstDerivativePenaltyImage);
  }
  
  // Delete the second order derivative penalty image
  if (m4p_secondDerivativePenaltyImage)
  {
    // Loop over time basis functions
    for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
    {
      if (m4p_secondDerivativePenaltyImage[tbf])
      {
        // Loop over respiratory basis functions
        for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
        {
          if (m4p_secondDerivativePenaltyImage[tbf][rbf])
          {
            free(m4p_secondDerivativePenaltyImage[tbf][rbf]);
          }
        }
        free(m4p_secondDerivativePenaltyImage[tbf]);
      }
    }
    free(m4p_secondDerivativePenaltyImage);
  }
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

void iOptimizerAPPGML::ShowHelpSpecific()
{
  cout << "This optimizer is the A-Penalized Preconditioned Gradient algorithm (A-PPGML) from M. Millardet et al, IEEE TRPMS," << endl;
  cout << "2021, vol. 6, pp. 629-640. It is based on the same idea as AML, which is allowing for negative values in the image" << endl;
  cout << "by shifting an algorithm enforcing positive values. APPGML can be seen as a shifted version of PPGML with a lower" << endl;
  cout << "bound A. This algorithm allows for negative image values in case the provided bound is also negative. The shift is" << endl;
  cout << "applied uniformly in the whole image unless a shift image is provided. This shift image allows to customize the" << endl;
  cout << "shift value applied in each voxel. The algorithm can also use penalized reconstruction because it is based on PPGML." << endl;
  cout << "The following options can be used (in this particular order when provided as a list):" << endl;
  cout << "  initial image value: to set the uniform voxel value for the initial image" << endl;
  cout << "  denominator threshold: to set the threshold of the data space denominator under which the ratio is set to 1" << endl;
  cout << "  minimum image update: to set the minimum of the image update factor under which it stays constant (0 or a negative value" << endl;
  cout << "                        means no minimum thus allowing a 0 update)" << endl;
  cout << "  maximum image update: to set the maximum of the image update factor over which it stays constant (0 or a negative value means" << endl;
  cout << "                        no maximum)" << endl;
  cout << "  bound: to set the lower bound value to apply the shift (must be negative to authorize negative values)" << endl;
  cout << "  shift image path: provide an image used to applied non-uniform shift within the reconstructed image. The values of the" << endl;
  cout << "                    shift image weight the shift applied to every voxel (0 means the voxel is not shifted, 1 means the" << endl;
  cout << "                    voxel is shifted by the bound, and t means the voxel is shifted by t times the bound). If no shift" << endl;
  cout << "                    image is provided, the shift is equal to the bound and is applied uniformly on the whole image" << endl;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerAPPGML::ReadConfigurationFile(const string& a_configurationFile)
{
  string key_word = "";
  string buffer = "";
  
  // Read the initial image value option
  key_word = "initial image value";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_initialValue, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerAPPGML::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Read the denominator threshold option
  key_word = "denominator threshold";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_dataSpaceDenominatorThreshold, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerAPPGML::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Read the minimum image update option
  key_word = "minimum image update";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_minimumImageUpdateFactor, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerAPPGML::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Read the maximum image update option
  key_word = "maximum image update";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_maximumImageUpdateFactor, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerAPPGML::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Read the bound value
  key_word = "bound";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_bound, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerAPPGML::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Read the path to image defining region to apply the shift
  key_word = "region image path";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_regionImagePath, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerAPPGML::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Normal end
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerAPPGML::ReadOptionsList(const string& a_optionsList)
{
  // Too complicated to do it that way
  Cerr("***** iOptimizerAPPGML::ReadOptionsList() -> Options can be specified only using a configuration file !" << endl);
  return 1;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerAPPGML::CheckSpecificParameters()
{
  // Check that initial image value is strictly positive
  if (m_initialValue<=0.)
  {
    Cerr("***** iOptimizerAPPGML->CheckSpecificParameters() -> Provided initial image value (" << m_initialValue << ") must be strictly positive !" << endl);
    return 1;
  }
  // Check that denominator threshold value is strictly positive
  if (m_dataSpaceDenominatorThreshold<=0.)
  {
    Cerr("***** iOptimizerAPPGML->CheckSpecificParameters() -> Provided data space denominator threshold (" << m_dataSpaceDenominatorThreshold << ") must be strictly positive !" << endl);
    return 1;
  }
  // Check that maximum image update factor is higher than the minimum
  if (m_minimumImageUpdateFactor>0. && m_maximumImageUpdateFactor>0. && m_maximumImageUpdateFactor<m_minimumImageUpdateFactor)
  {
    Cerr("***** iOptimizerAPPGML->CheckSpecificParameters() -> Provided minimum/maximum (" << m_minimumImageUpdateFactor << "/" << m_maximumImageUpdateFactor << " are inconsistent !" << endl);
    return 1;
  }
  // Check that the initial image value is above the provided bound
  if (m_initialValue<m_bound)
  {
    Cerr("***** iOptimizerAPPGML::CheckSpecificParameters() -> The initial image value (" << m_initialValue << ") must be higher than or equal to the provided bound value (" << m_bound << ") !" << endl);
    return 1;
  }
  // Normal end
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerAPPGML::InitializeSpecific()
{
  // Allocate and create the image (not shifted)
  m4p_imageNotShifted = (FLTNB****)malloc(mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions()*sizeof(FLTNB***));
  // Loop over time basis functions
  for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
  { 
    m4p_imageNotShifted[tbf] = (FLTNB***)malloc(mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions()*sizeof(FLTNB**));
    // Loop over respiratory basis functions
    for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
    {
      m4p_imageNotShifted[tbf][rbf] = (FLTNB**)malloc(mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions()*sizeof(FLTNB*));
      // Loop over cardiac basis functions
      for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++)
      {
        m4p_imageNotShifted[tbf][rbf][cbf] = mp_ImageSpace -> AllocateMiscellaneousImage(); //Get a pointer to a newly allocated image
        //for (int v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++) m4p_imageNotShifted[tbf][rbf][cbf][v] = -1;
      }
    }
  }
  
  // Allocate and create the penalty image
  m4p_firstDerivativePenaltyImage = (FLTNB****)malloc(mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions()*sizeof(FLTNB***));
  // Loop over time basis functions
  for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
  { 
    m4p_firstDerivativePenaltyImage[tbf] = (FLTNB***)malloc(mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions()*sizeof(FLTNB**));
    // Loop over respiratory basis functions
    for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
    {
      m4p_firstDerivativePenaltyImage[tbf][rbf] = (FLTNB**)malloc(mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions()*sizeof(FLTNB*));
      // Loop over cardiac basis functions
      for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++)
      {
        m4p_firstDerivativePenaltyImage[tbf][rbf][cbf] = mp_ImageSpace -> AllocateMiscellaneousImage(); //Get a pointer to a newly allocated image
      }
    }
  }
  
  // Allocate and create the second order derivative penalty image
  m4p_secondDerivativePenaltyImage = (FLTNB****)malloc(mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions()*sizeof(FLTNB***));
  // Loop over time basis functions
  for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
  { 
    m4p_secondDerivativePenaltyImage[tbf] = (FLTNB***)malloc(mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions()*sizeof(FLTNB**));
    // Loop over respiratory basis functions
    for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
    {
      m4p_secondDerivativePenaltyImage[tbf][rbf] = (FLTNB**)malloc(mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions()*sizeof(FLTNB*));
      // Loop over cardiac basis functions
      for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++)
      {
        m4p_secondDerivativePenaltyImage[tbf][rbf][cbf] = mp_ImageSpace -> AllocateMiscellaneousImage(); //Get a pointer to a newly allocated image
      }
    }
  }  

  // Allocate the specific region of interest
  mp_specificRegionOfInterest = mp_ImageSpace -> AllocateMiscellaneousImage(); // Get a pointer to a newly allocated image
  mp_specificRegionOfInterestConvolved = mp_ImageSpace -> AllocateMiscellaneousImage(); // Get a pointer to a newly allocated image
  // Initialize/Read the specific region of interest
  if (m_regionImagePath == "") // whole phantom to be considered, so uniform image
  {
    for (int v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++)
    {
        mp_specificRegionOfInterest[v] = 1.;
        mp_specificRegionOfInterestConvolved[v] = 1.;
    }
  }
  else
  {
    // Interfile image reading (INTF_LERP_ENABLED = interpolation allowed)
    if (IntfReadImage(m_regionImagePath, mp_specificRegionOfInterest, mp_ImageDimensionsAndQuantification, m_verbose, INTF_LERP_ENABLED))
    {
      Cerr("***** iOptimizerAPPGML::InitializeSpecific() -> Error reading interfile image '" << m_regionImagePath << "' !" << endl);
      return 1;
    }
    if (IntfReadImage(m_regionImagePath, mp_specificRegionOfInterestConvolved, mp_ImageDimensionsAndQuantification, m_verbose, INTF_LERP_ENABLED))
    {
      Cerr("***** iOptimizerAPPGML::InitializeSpecific() -> Error reading interfile image '" << m_regionImagePath << "' !" << endl);
      return 1;
    }
  }  
  // Verbose
  if (m_verbose>=2)
  {
    Cout("iOptimizerAPPGML::InitializeSpecific() -> Use the APPGML optimizer" << endl);
    if (m_verbose>=3)
    {
      Cout("  --> Initial image value: " << m_initialValue << endl);
      Cout("  --> Data space denominator threshold: " << m_dataSpaceDenominatorThreshold << endl);
      if (m_minimumImageUpdateFactor>0.) Cout("  --> Minimum image update factor: " << m_minimumImageUpdateFactor << endl);
      else Cerr("!!!!! The minimum update value is not set, if using subsets, voxels could be trapped in 0 value causing some negative bias !" << endl);
      if (m_maximumImageUpdateFactor>0.) Cout("  --> Maximum image update factor: " << m_maximumImageUpdateFactor << endl);
      Cout("  --> Bound value: " << m_bound << endl);
      if (m_regionImagePath != "") Cout("  --> Shift weighted for each voxel by provided image: " << m_regionImagePath << endl);
      else Cout("  --> Uniform shift across the whole image" << endl);
    }
  }
  // Normal end
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================


int iOptimizerAPPGML::PreDataUpdateSpecificStep()
{
  // Apply convolver in pre iteration to region of interest because image-based PSF not inside ForwardProject function
  if (m_isInPreIteration)
  {
    if (mp_ImageConvolverManager->ApplyConvolution(mp_specificRegionOfInterestConvolved))
    {
      Cerr("***** iOptimizerAPPGML::PreDataUpdateSpecificStep() -> A problem occurred while applying image convolver to region to apply shift !" << endl);
      return 1;
    }
  }
  if (!m_isInPreIteration)
  {
    // Set the number of threads
    #ifdef CASTOR_OMP
    omp_set_num_threads(mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation());
    #endif

    // Loop over time basis functions
    for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
    {
      // Loop over respiratory basis functions
      for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
      {
        // Loop over cardiac basis functions
        for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++)
        {
          int v;
          // multi-threading
          #pragma omp parallel for private(v) schedule(guided)
          for (v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++)
          {

            // Store the image without shift for penalty computation
            m4p_imageNotShifted[tbf][rbf][cbf][v] = mp_ImageSpace-> m4p_image[tbf][rbf][cbf][v];
            mp_ImageSpace-> m4p_image[tbf][rbf][cbf][v] -= mp_specificRegionOfInterest[v]*m_bound;
            mp_ImageSpace-> m4p_forwardImage[tbf][rbf][cbf][v] -= mp_specificRegionOfInterest[v]*m_bound;
          }
        }
      }
    }
  }
  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerAPPGML::SensitivitySpecificOperations( FLTNB a_data, FLTNB a_forwardModel, FLTNB* ap_weight,
                                                   FLTNB a_multiplicativeCorrections, FLTNB a_additiveCorrections, FLTNB a_blankValue,
                                                   FLTNB a_quantificationFactor, oProjectionLine* ap_Line )
{
  // Line weight here is simply 1
  *ap_weight = 1.;
  // That's all
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerAPPGML::DataSpaceSpecificOperations( FLTNB a_data, FLTNB a_forwardModel, FLTNB* ap_backwardValues,
                                                 FLTNB a_multiplicativeCorrections, FLTNB a_additiveCorrections, FLTNB a_blankValue,
                                                 FLTNB a_quantificationFactor, oProjectionLine* ap_Line )
{
  if (!m_isInPreIteration)
  {
    // Compute the bound projection
    FLTNB shift = ForwardProject(ap_Line,mp_specificRegionOfInterestConvolved) * m_bound;
    a_data -= shift;
    // Truncate data to 0 if negative
    if (a_data<0.) a_data = 0.;

    // Compute numerator
    FLTNB numerator = a_data - a_forwardModel;
    // Compute denominator that will be strictly positive
    FLTNB denominator = max(a_forwardModel,m_dataSpaceDenominatorThreshold);
    // Update backward values
    *ap_backwardValues = numerator / denominator;
  }
  // That's all
  return 0;

}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerAPPGML::PreImageUpdateSpecificStep()
{
  if (!m_isInPreIteration)
  {
    // ==========================================================================================
    // If no penalty, then exit (the penalty image term has been initialized to 0)
    if (mp_Penalty==NULL) return 0;
    // Set the number of threads
    #ifdef CASTOR_OMP
    omp_set_num_threads(mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation());
    #endif
    // Verbose
    if (m_verbose>=1) Cout("iOptimizerAPPGML::PreImageUpdateSpecificStep() -> Compute penalty term" << endl);
    // ==========================================================================================
    // Global precomputation step if needed by the penalty
    if (mp_Penalty->GlobalPreProcessingStep())
    {
      Cerr("***** iOptimizerAPPGML::PreImageUpdateSpecificStep() -> A problem occurred while computing the penalty pre-processing step !" << endl);
      return 1;
    }
    // ==========================================================================================
    // Loop over time basis functions
    for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
    {
      // Loop over respiratory basis functions
      for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
      {
        // Loop over cardiac basis functions
        for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++)
        {
          // In order to detect problems in the multi-threaded loop
          bool problem = false;
          // Voxel index
          INTNB v;
          // Multi-threading over voxels
          #pragma omp parallel for private(v) schedule(guided)
          for (v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++)
          { 
            // Get the thread index
            int th = 0;
            #ifdef CASTOR_OMP
            th = omp_get_thread_num();
            #endif
            // Local precomputation step if needed by the penalty
            if (mp_Penalty->LocalPreProcessingStep(tbf,rbf,cbf,v,th))
            {
              Cerr("***** iOptimizerAPPGML::PreImageUpdateSpecificStep() -> A problem occurred while computing the penalty local pre-processing step for voxel " << v << " !" << endl);
              problem = true;
            }
            // Compute first and second derivative order penalty terms on image (not shifted)
            m4p_firstDerivativePenaltyImage[tbf][rbf][cbf][v] = mp_Penalty->ComputeFirstDerivative(m4p_imageNotShifted[tbf][cbf][rbf],v,th);
            m4p_secondDerivativePenaltyImage[tbf][rbf][cbf][v] = mp_Penalty->ComputeSecondDerivative(m4p_imageNotShifted[tbf][cbf][rbf],v,th);
          }
          // Check for problems
          if (problem)
          {
            Cerr("***** iOptimizerAPPGML::PreImageUpdateSpecificStep() -> A problem occurred inside the multi-threaded loop, stop now !" << endl);
            return 1;
          }
        }
      }
    }
  }
  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerAPPGML::ImageSpaceSpecificOperations( FLTNB a_currentImageValue, FLTNB* ap_newImageValue,
                                                  FLTNB a_sensitivity, FLTNB* ap_correctionValues, INTNB a_voxel, int tbf, int rbf, int cbf )
{
  if (!m_isInPreIteration)
  {
    FLTNB image_update_factor = *ap_correctionValues - m4p_firstDerivativePenaltyImage[tbf][rbf][cbf][a_voxel] / (FLTNB)(mp_nbSubsets[m_currentIteration]); //divide the penalty by the number of subset, otherwise the regularization parameter has to be adapted when we change the number of subsets
    
    a_sensitivity += a_currentImageValue * m4p_secondDerivativePenaltyImage[tbf][rbf][cbf][a_voxel] / (FLTNB)(mp_nbSubsets[m_currentIteration]);
    
    image_update_factor = 1. + image_update_factor / a_sensitivity;
    // Apply minimum image update factor
    if ( m_minimumImageUpdateFactor > 0. && image_update_factor < m_minimumImageUpdateFactor ) image_update_factor = m_minimumImageUpdateFactor;
    // Apply maximum image update factor
    if ( m_maximumImageUpdateFactor > 0. && image_update_factor > m_maximumImageUpdateFactor ) image_update_factor = m_maximumImageUpdateFactor;
    // Update image
    *ap_newImageValue = a_currentImageValue * image_update_factor;
    // Check if it is a number, if no, keep the value unchanged
    if (!isfinite(*ap_newImageValue)) *ap_newImageValue = a_currentImageValue;
    // Apply the shift to region of interest to obtain final image (will be shifted in next iteration in PreDataUpdateSpecificStep)
    *ap_newImageValue += mp_specificRegionOfInterest[a_voxel] * m_bound;
  }
  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================