Merge branch 'RTSMarker' into 'ScoutRT_Qt6'

RTSMarker bug

See merge request cpt_bioeng/drt!43

reg23Topograms/itkDTRrecon/CMakeLists.txt

@@ -1,4 +1,4 @@
-cmake_minimum_required(VERSION 2.8.12)
+cmake_minimum_required(VERSION 3.16)
 
 SET(LIB_NAME itkDTRrecon)
 

reg23Topograms/itkDTRrecon/autoreg/itkDRTHelpers.h

@@ -6,12 +6,6 @@
 #include <cstring>
 #include <iterator>
 
-#ifdef __GNUC__
-#define VARIABLE_IS_NOT_USED __attribute__ ((unused))
-#else
-#define VARIABLE_IS_NOT_USED
-#endif
-
 //Function to get method and class name for logging purposes.
 template <size_t FL, size_t PFL>
 const char* computeMethodName(const char (&function)[FL], const char (&prettyFunction)[PFL])
@@ -30,7 +24,7 @@ const char* computeMethodName(const char (&function)[FL], const char (&prettyFun
 
 namespace itk {
 /* reference string required for comparison with tag values */
-static const char VARIABLE_IS_NOT_USED *ImageOrientationStrings[] = {
+static const char* ImageOrientationStrings[] = {
     "NotDefined",
     "HFS",
     "FFS",
@@ -46,49 +40,49 @@ static const bool verbose = false;
 /* this is in the end a IEC to HFS...
  * but we keep this for ourselves...
  */
-static double VARIABLE_IS_NOT_USED Standard_DRT2LPS[9] = {
+static double Standard_DRT2LPS[9] = {
     1, 0, 0,
     0, 0, -1,
     0, 1, 0
 };
 
-static double VARIABLE_IS_NOT_USED PAElementsIEC[9] = {
+static double PAElementsIEC[9] = {
     1, 0, 0,
     0, -1, 0,
     0, 0, -1
 };
 
-static double VARIABLE_IS_NOT_USED LATElementsIEC[9] = {
+static double LATElementsIEC[9] = {
     0, 0, -1,
     0, -1, 0,
     -1, 0, 0
 };
 
-static double VARIABLE_IS_NOT_USED HFS2IEC[9] = {
+static double HFS2IEC[9] = {
     1, 0, 0,
     0, 0, 1,
     0, -1, 0
 };
 
-static double VARIABLE_IS_NOT_USED FFS2IEC[9] = {
+static double FFS2IEC[9] = {
     -1, 0, 0,
     0, 0, -1,
     0, -1, 0
 };
 
-static double VARIABLE_IS_NOT_USED HFP2IEC[9] = {
+static double HFP2IEC[9] = {
     -1, 0, 0,
     0, 0, 1,
     0, 1, 0
 };
 
-static double VARIABLE_IS_NOT_USED FFP2IEC[9] = {
+static double FFP2IEC[9] = {
     1, 0, 0,
     0, 0, -1,
     0, 1, 0
 };
 
-static double VARIABLE_IS_NOT_USED PAT2IEC[9] = {
+static double PAT2IEC[9] = {
     1, 0, 0,
     0, 0, 1,
     0, -1, 0

reg23Topograms/itkDTRrecon/autoreg/itkMutualInformationTwoImageToOneImageMetric.hxx

@@ -143,7 +143,7 @@ MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue()
     //std::cout<< "-----> Mutual :: fixedImageRegion1: "<< metric1->GetFixedImageRegion() << std::endl;
 
     auto movingImageRegion = this->m_Filter1->GetOutput()->GetBufferedRegion();
-    //unsigned int numberOfPixels = movingImageRegion.GetNumberOfPixels();
+    unsigned int numberOfPixels = movingImageRegion.GetNumberOfPixels();
 
     //auto numberOfSamples = static_cast<unsigned int>(numberOfPixels * 0.50); //100%
     // since we have a ROI, then we should not set allPixels to TRUE.
@@ -180,7 +180,7 @@ MutualInformationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetValue()
     //std::cout<< "-----> Mutual :: fixedImageRegion2: "<< metric2->GetFixedImageRegion() << std::endl;
 
     movingImageRegion = this->m_Filter2->GetOutput()->GetBufferedRegion();
-    //numberOfPixels = movingImageRegion.GetNumberOfPixels();
+    numberOfPixels = movingImageRegion.GetNumberOfPixels();
 
     //numberOfSamples = static_cast<unsigned int>(numberOfPixels * 0.50); //100%
     //metric2->SetNumberOfSpatialSamples(numberOfSamples);

reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.h

@@ -20,7 +20,9 @@
 
 #include "itkCovariantVector.h"
 #include "itkMacro.h"
-#include "itkNormalizedCorrelationImageToImageMetric.hxx"
+//#include "itkNormalizedCorrelationImageToImageMetric.hxx"
+#include "itkNormalizedCorrelationImageToImageMetric.h"
+
 #include "itkPoint.h"
 #include "itkgTwoImageToOneImageMetric.h"
 

reg23Topograms/itkDTRrecon/autoreg/itkNormalizedCorrelationTwoImageToOneImageMetric.hxx

@@ -262,7 +262,7 @@ NormalizedCorrelationTwoImageToOneImageMetric<TFixedImage, TMovingImage>::GetVal
     }
 
     // Calculate the measure value between fixed image 1 and the moving image
-    //auto oldprecision = std::cout.precision();
+    auto oldprecision = std::cout.precision();
 
 
 //    std::cout<<"Region "  <<this->GetFixedImageRegion1() <<std::endl;

reg23Topograms/itkDTRrecon/itkImageProcessor.cpp

@@ -2127,13 +2127,14 @@ void itkImageProcessor::Write2DImages(){
 
 vtkImageData* itkImageProcessor::GetLocalizer1VTK()
 {
-
+    // Rescale the intensity of the projection images to 0-32768 for output.
     using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
 
     RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
     rescaler1->SetOutputMinimum(0);
-    rescaler1->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
-    rescaler1->SetInput( m_PASourceDupli->GetOutput() );
+    rescaler1->SetOutputMaximum(32768);
+    rescaler1->SetInput(m_PASourceDupli->GetOutput());
+
     rescaler1->Update();
 
     toVTKLocalizer1->SetInput(rescaler1->GetOutput());
@@ -2186,18 +2187,22 @@ vtkImageData* itkImageProcessor::GetLocalizer1VTK()
 vtkImageData* itkImageProcessor::GetLocalizer2VTK()
 {
 
+
+    // Rescale the intensity of the projection images to 0-32768 for output.
     using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
 
     RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
     rescaler2->SetOutputMinimum(0);
-    rescaler2->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
-    rescaler2->SetInput( m_LATSourceDupli->GetOutput() );
+    rescaler2->SetOutputMaximum(32768);
+    rescaler2->SetInput(m_LATSourceDupli ->GetOutput());
+
     rescaler2->Update();
 
     toVTKLocalizer2->SetInput(rescaler2->GetOutput());
     toVTKLocalizer2->Update();
 
 
+
 //    if(true) {
 //        using ImageRegionType3D = ImageType3D::RegionType;
 //        using SizeType3D = ImageRegionType3D::SizeType;
@@ -2241,7 +2246,6 @@ vtkImageData* itkImageProcessor::GetLocalizer2VTK()
             toVTKLocalizer2->GetOutput();
 }
 
-
 vtkImageData* itkImageProcessor::GetProjection1VTK()
 {
 
@@ -2252,60 +2256,33 @@ vtkImageData* itkImageProcessor::GetProjection1VTK()
         //TODO
     }
 
+    // Rescale the intensity of the projection images to 0-32768 for output.
     using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
 
     RescaleFilterType::Pointer rescaler1 = RescaleFilterType::New();
     rescaler1->SetOutputMinimum(0);
-    rescaler1->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
+    rescaler1->SetOutputMaximum(32768);
     rescaler1->SetInput( imageDRT1In );
     rescaler1->Update();
 
+    //    using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator<InternalImageType>;
+    //    auto imageCalculatorFilter = ImageCalculatorFilterType::New();
+    //    imageCalculatorFilter->SetImage(imageDRT1In);
+    //    imageCalculatorFilter->Compute();
+
+    //    std::cout<< "itkImageProcessor::imageDRT1In() " <<
+    //                imageCalculatorFilter <<std::endl;
+
+    //    using ImageCalculatorFilterType2 = itk::MinimumMaximumImageCalculator<OutputImageType>;
+    //    auto imageCalculatorFilter2 = ImageCalculatorFilterType2::New();
+    //    imageCalculatorFilter2->SetImage(rescaler1->GetOutput());
+    //    imageCalculatorFilter2->Compute();
+    //    std::cout<< "itkImageProcessor::imageDRT2In() " <<
+    //                imageCalculatorFilter2 <<std::endl;
+
     toVTK2D1->SetInput(rescaler1->GetOutput());
     toVTK2D1->Update();
 
-    return
-            toVTK2D1->GetOutput();
-}
-
-vtkImageData* itkImageProcessor::GetProjection1VTKToWrite()
-{
-
-    if(m_DRTImage1MetaInfo == NULL ||
-            m_DRTGeometryMetaInfo == NULL ||
-            m_TransformMetaInfo == NULL     ){
-        return NULL;
-        //TODO
-    }
-
-    using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator<InternalImageType>;
-    auto imageCalculatorFilter = ImageCalculatorFilterType::New();
-    imageCalculatorFilter->SetImage(imageDRT1In);
-    imageCalculatorFilter->Compute();
-
-    using IntWindowType = itk::IntensityWindowingImageFilter<InternalImageType, OutputImageType>;
-    auto intWindowFilter = IntWindowType::New();
-    intWindowFilter->SetWindowMinimum(imageCalculatorFilter->GetMinimum());
-    intWindowFilter->SetWindowMaximum(imageCalculatorFilter->GetMaximum());
-
-    intWindowFilter->SetOutputMinimum(0);
-    if(imageCalculatorFilter->GetMaximum() > SHRT_MAX){
-        intWindowFilter->SetOutputMaximum(SHRT_MAX);
-    }else
-        intWindowFilter->SetOutputMaximum(imageCalculatorFilter->GetMaximum());
-
-
-    intWindowFilter->SetInput(imageDRT1In);
-    intWindowFilter->Update();
-
-    toVTK2D1->SetInput(intWindowFilter->GetOutput());
-    toVTK2D1->Update();
-
-//    using ImageCalculatorFilterType2 = itk::MinimumMaximumImageCalculator<OutputImageType>;
-//    auto imageCalculatorFilter2 = ImageCalculatorFilterType2::New();
-//    imageCalculatorFilter2->SetImage(intWindowFilter->GetOutput());
-//    imageCalculatorFilter2->Compute();
-//    std::cout<< "itkImageProcessor::imageCalculatorFilter2() " <<
-//                imageCalculatorFilter2->GetMinimum() << " " << imageCalculatorFilter2->GetMaximum() <<std::endl;
 
     //    using ImageRegionType3D = ImageType3D::RegionType;
     //    using SizeType3D = ImageRegionType3D::SizeType;
@@ -2410,63 +2387,19 @@ vtkImageData* itkImageProcessor::GetProjection2VTK()
         //TODO
     }
 
-
+    // Rescale the intensity of the projection images to 0-32768 for output.
     using RescaleFilterType = itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;
 
     RescaleFilterType::Pointer rescaler2 = RescaleFilterType::New();
     rescaler2->SetOutputMinimum(0);
-    rescaler2->SetOutputMaximum(IMG_VIS_MAXIMUM_RANGE);
+    rescaler2->SetOutputMaximum(32768);
     rescaler2->SetInput( imageDRT2In );
-    rescaler2->Update();
 
+    rescaler2->Update();
     toVTK2D2->SetInput(rescaler2->GetOutput());
     toVTK2D2->Update();
 
-    return
-            toVTK2D2->GetOutput();
-}
-
-
-
-vtkImageData* itkImageProcessor::GetProjection2VTKToWrite()
-{
-    if(m_DRTImage2MetaInfo == NULL ||
-            m_DRTGeometryMetaInfo == NULL ||
-            m_TransformMetaInfo == NULL     ){
-        return NULL;
-        //TODO
-    }
-
-    using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator<InternalImageType>;
-    auto imageCalculatorFilter = ImageCalculatorFilterType::New();
-    imageCalculatorFilter->SetImage(imageDRT2In);
-    imageCalculatorFilter->Compute();
-
-    using IntWindowType = itk::IntensityWindowingImageFilter<InternalImageType, OutputImageType>;
-    auto intWindowFilter = IntWindowType::New();
-    intWindowFilter->SetWindowMinimum(imageCalculatorFilter->GetMinimum());
-    intWindowFilter->SetWindowMaximum(imageCalculatorFilter->GetMaximum());
-
-    intWindowFilter->SetOutputMinimum(0);
-    if(imageCalculatorFilter->GetMaximum() > SHRT_MAX)
-        intWindowFilter->SetOutputMaximum(SHRT_MAX);
-    else
-        intWindowFilter->SetOutputMaximum(imageCalculatorFilter->GetMaximum());
-
-
-    intWindowFilter->SetInput(imageDRT2In);
-    intWindowFilter->Update();
-
-
-    toVTK2D2->SetInput(intWindowFilter->GetOutput());
-    toVTK2D2->Update();
-
-//    using ImageCalculatorFilterType2 = itk::MinimumMaximumImageCalculator<OutputImageType>;
-//    auto imageCalculatorFilter2 = ImageCalculatorFilterType2::New();
-//    imageCalculatorFilter2->SetImage(intWindowFilter->GetOutput());
-//    imageCalculatorFilter2->Compute();
-//    std::cout<< "itkImageProcessor::imageCalculatorFilter2() " <<
-//                imageCalculatorFilter2->GetMinimum() << " " << imageCalculatorFilter2->GetMaximum() <<std::endl;
+    //    rescaler2->Print(std::cout);
 
     //    using ImageRegionType3D = ImageType3D::RegionType;
     //    using SizeType3D = ImageRegionType3D::SizeType;

reg23Topograms/itkDTRrecon/itkImageProcessor.h

@@ -22,7 +22,6 @@ gfattori 08.11.2021
 #include "itkResampleImageFilter.h"
 #include "itkCastImageFilter.h"
 #include "itkRescaleIntensityImageFilter.h"
-#include "itkIntensityWindowingImageFilter.h"
 #include "itkChangeInformationImageFilter.h"
 #include "itkGDCMImageIO.h"
 #include "itkMetaDataObject.h"
@@ -50,8 +49,6 @@ gfattori 08.11.2021
 #include "itkReg23.h"
 #include "itkReg23MetaInformation.h"
 
-#define IMG_VIS_MAXIMUM_RANGE 2048
-
 namespace itk
 {
 
@@ -245,10 +242,6 @@ public:
     vtkImageData* GetProjection1VTK();
     vtkImageData* GetProjection2VTK();
 
-    /** Conveniency methods to get VTK images for rendering */
-    vtkImageData* GetProjection1VTKToWrite();
-    vtkImageData* GetProjection2VTKToWrite();
-
     vtkImageData* GetLocalizer1VTK();
     vtkImageData* GetLocalizer2VTK();
 

reg23Topograms/itkDTRrecon/itkgSiddonJacobsRayCastInterpolateImageFunction.hxx

@@ -48,7 +48,6 @@ CIT 6, 89-94 (1998).
 #define itkgSiddonJacobsRayCastInterpolateImageFunction_hxx
 
 #include "itkgSiddonJacobsRayCastInterpolateImageFunction.h"
-#include "itkContinuousIndex.h"
 
 #include "itkMath.h"
 #include <cstdlib>
@@ -129,7 +128,7 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
     float alphaX, alphaY, alphaZ, alphaCmin, alphaCminPrev;
     float alphaUx, alphaUy, alphaUz;
     float alphaIntersectionUp[3], alphaIntersectionDown[3];
-    float d12, value,valuetril;
+    float d12, value;
     float firstIntersectionIndex[3];
     int   firstIntersectionIndexUp[3], firstIntersectionIndexDown[3];
     int   iU, jU, kU;
@@ -156,12 +155,9 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
     }
 
     PointType PointReq = point;
-    //std::cout<<"PointReq: "<<point[0] <<" "<<point[1] <<" "<<point[2] <<" "<<std::endl;
     PointReq[0] += m_PanelOffset;
 
     drrPixelWorld = m_InverseTransform->TransformPoint(PointReq);
-    //std::cout<<"drrPixelWorld: "<<drrPixelWorld[0] <<" "<<drrPixelWorld[1] <<" "<<drrPixelWorld[2] <<" "<<std::endl;
-
 
     // Get ths input pointers
     InputImageConstPointer inputPtr = this->GetInputImage();
@@ -185,7 +181,6 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
     SlidingSourcePoint[2] = 0.;
 
     PointType SourceWorld = m_InverseTransform->TransformPoint(SlidingSourcePoint);
-    //std::cout<<"SourceWorld: "<<SourceWorld[0] <<" "<<SourceWorld[1] <<" "<<SourceWorld[2] <<" "<<std::endl;
 
     PointType O(3);
     O[0] = -ctPixelSpacing[0]/2.;
@@ -435,7 +430,6 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
     cIndex[1] = firstIntersectionIndexDown[1];
     cIndex[2] = firstIntersectionIndexDown[2];
 
-
     while (alphaCmin < alphaMax) /* Check if the ray is still in the CT volume */
     {
         /* Store the current ray position */
@@ -463,104 +457,17 @@ gSiddonJacobsRayCastInterpolateImageFunction<TInputImage, TCoordRep>::Evaluate(c
             alphaZ = alphaZ + alphaUz;
         }
 
-
         if ((cIndex[0] >= 0) && (cIndex[0] < static_cast<IndexValueType>(sizeCT[0])) && (cIndex[1] >= 0) &&
-            (cIndex[1] < static_cast<IndexValueType>(sizeCT[1])) && (cIndex[2] >= 0) &&
-            (cIndex[2] < static_cast<IndexValueType>(sizeCT[2])))
+                (cIndex[1] < static_cast<IndexValueType>(sizeCT[1])) && (cIndex[2] >= 0) &&
+                (cIndex[2] < static_cast<IndexValueType>(sizeCT[2])))
         {
-            // Calculate entry and exit points using alphaCmin and alphaCminPrev
+            /* If it is a valid index, get the voxel intensity. */
             value = static_cast<float>(inputPtr->GetPixel(cIndex));
-            // Accumulate the interpolated intensity along the ray path
             if (value > m_Threshold) /* Ignore voxels whose intensities are below the threshold. */
             {
-
-                // Move along the ray by alphaCminPrev to find the entry point of this voxel
-                PointType entryPoint;
-                entryPoint[0] = SourceWorld[0] + alphaCminPrev * rayVector[0] ;
-                entryPoint[1] = SourceWorld[1] + alphaCminPrev * rayVector[1] ;
-                entryPoint[2] = SourceWorld[2] + alphaCminPrev * rayVector[2] ;
-
-                // Move along the ray by alphaCmin to find the exit point of this voxel
-                PointType exitPoint;
-                exitPoint[0] = SourceWorld[0] + alphaCmin * rayVector[0] ;
-                exitPoint[1] = SourceWorld[1] + alphaCmin * rayVector[1] ;
-                exitPoint[2] = SourceWorld[2] + alphaCmin * rayVector[2] ;
-
-                // Get the mid-point of the voxel / ray interception
-                PointType midpoint;
-                midpoint[0]= (entryPoint[0] + exitPoint[0]) * 0.5;
-                midpoint[1]= (entryPoint[1] + exitPoint[1]) * 0.5;
-                midpoint[2]= (entryPoint[2] + exitPoint[2]) * 0.5;
-
-                // Ray is computed in 'Siddon' geometry with zero origin,
-                // whereas the continuous index will be computed from the input
-                // image. The origin and shifts to the voxel edge are to be account for
-                midpoint[0] += ctOrigin[0] ;//+ O[0];
-                midpoint[1] += ctOrigin[1] ;//+ O[1];
-                midpoint[2] += ctOrigin[2] ;//+ O[2];
-
-                // Convert mid-point phyisical point into continuous index
-                // We need to use this position to find the neighbouring voxels
-                // for trilinear interpolation - not the cIndex!
-                itk::ContinuousIndex <double,3> continuousIndex;
-                inputPtr->TransformPhysicalPointToContinuousIndex(midpoint,continuousIndex);
-
-                // Get the baseIndex by flooring the continuous index
-                IndexType baseIndex;
-                baseIndex[0]=static_cast<IndexValueType>(std::floor(continuousIndex[0]));
-                baseIndex[1]=static_cast<IndexValueType>(std::floor(continuousIndex[1]));
-                baseIndex[2]=static_cast<IndexValueType>(std::floor(continuousIndex[2]));
-
-                // Calculate fractional parts for interpolation at the midpoint
-                double x_frac = continuousIndex[0] - baseIndex[0];
-                double y_frac = continuousIndex[1] - baseIndex[1];
-                double z_frac = continuousIndex[2] - baseIndex[2];
-
-                // Perform boundary checks for trilinear interpolation
-                bool within_bounds = (baseIndex[0] >= 0 && baseIndex[0] < sizeCT[0] - 1) &&
-                                     (baseIndex[1] >= 0 && baseIndex[1] < sizeCT[1] - 1) &&
-                                     (baseIndex[2] >= 0 && baseIndex[2] < sizeCT[2] - 1);
-                if (within_bounds)
-                {
-                    // Fetch intensities from neighboring voxels for trilinear interpolation
-                    float c000 = static_cast<float>(inputPtr->GetPixel(baseIndex));
-                    baseIndex[0] += 1;
-                    float c100 = static_cast<float>(inputPtr->GetPixel(baseIndex));
-                    baseIndex[1] += 1;
-                    float c110 = static_cast<float>(inputPtr->GetPixel(baseIndex));
-                    baseIndex[0] -= 1;
-                    float c010 = static_cast<float>(inputPtr->GetPixel(baseIndex));
-                    baseIndex[2] += 1;
-                    float c011 = static_cast<float>(inputPtr->GetPixel(baseIndex));
-                    baseIndex[0] += 1;
-                    float c111 = static_cast<float>(inputPtr->GetPixel(baseIndex));
-                    baseIndex[1] -= 1;
-                    float c101 = static_cast<float>(inputPtr->GetPixel(baseIndex));
-                    baseIndex[0] -= 1;
-                    float c001 = static_cast<float>(inputPtr->GetPixel(baseIndex));
-
-                    // NOTE: do not use baseIndex anymore after this. It's messed up.
-
-                    // Perform trilinear interpolation at the midpoint
-                    float c00 = c000 * (1 - x_frac) + c100 * x_frac;
-                    float c01 = c001 * (1 - x_frac) + c101 * x_frac;
-                    float c10 = c010 * (1 - x_frac) + c110 * x_frac;
-                    float c11 = c011 * (1 - x_frac) + c111 * x_frac;
-
-                    float c0 = c00 * (1 - y_frac) + c10 * y_frac;
-                    float c1 = c01 * (1 - y_frac) + c11 * y_frac;
-
-                    valuetril = c0 * (1 - z_frac) + c1 * z_frac;
-
-                    d12 += (alphaCmin - alphaCminPrev) * (valuetril - m_Threshold);
-
-                } else {
-                    d12 += (alphaCmin - alphaCminPrev) * (value - m_Threshold);
-                }
-
+                d12 += (alphaCmin - alphaCminPrev) * (value - m_Threshold);
             }
         }
-
     }
 
     if (d12 < minOutputValue)