OpenCV canny vs my own C++ version - high difference

Question

I have implemented my own canny edge detector in C++. If I compare the result with the same thresholds to OpenCV (cv::Canny(img, res, 50, 100)), the difference is relatively high (left - my version, right - OpenCV). Is there any trick in OpenCV? I have tried to understand its code, but I got lost.

I have used test image of chameleon, where I can see the algorithm steps as well. My result looks comparable with this one on Wikipedia.

My code (I am using my own Image2d class, but the methods should be clear by its names)

template <typename T>
void CannyEdgeDetector<T>::Run(T low, T high, bool l2Magnitude){    
    this->l2Magnitude = l2Magnitude;

    if constexpr (std::is_same<T, uint8_t>::value){
        this->low = low / 255.0f;
        this->high = high / 255.0f;
    }
    else{
        this->low = low;
        this->high = high;
    }

  
    this->ComputeGradients();
    auto nonMaxSup = this->NonMaximaSuppression();
    this->Hysteresis(nonMaxSup);
}

template <typename T>
void CannyEdgeDetector<T>::ComputeGradients(){
    //When the size of the kernel is 3, the Sobel kernel shown above may produce noticeable inaccuracies 
    //(after all, Sobel is only an approximation of the derivative). 
    //OpenCV addresses this inaccuracy for kernels of size 3 by using the Scharr() function
    //use kernel of size 3 by default

    ImageDerivator derivator(ImageDerivator::FilterType::SCHARR);
    derivator.Run(input, ImageDerivator::Derivatives::dX |
        ImageDerivator::Derivatives::dY);

    dx = std::move(derivator.dx);
    dy = std::move(derivator.dy);
}

template <typename T>
Image2d<float> CannyEdgeDetector<T>::NonMaximaSuppression(){
    
    Image2d<float> mag = (l2Magnitude) ? 
        ImageUtils::CalcL2Magnitude(dx, dy) : //sqrt(x*x + y*y)
        ImageUtils::CalcL1Magnitude(dx, dy);  //abs(x) + abs(y)         
    mag.MapToInterval(0, 1); //normalize magnitudes to [0, 1] interval

    Image2d<float> angle = ImageUtils::CalcAngle(dx, dy, true); //calc atan2(y, x)

    Image2d<float> res = mag.CreateDeepCopy();

    int neighbor1X, neighbor1Y, neighbor2X, neighbor2Y;

    for (int y = 1; y < input.GetHeight() - 1; y++) {
        for (int x = 1; x < input.GetWidth() - 1; x++) {

            float gradientAngle = std::abs(angle.GetPixelStart(x, y)[0]);            
      
            //calculate center of intervals
            //eg. 22.5 < gradientAngle <= 67.5
            //center is 45
            //offset is: [1, 1] and in negative half [-1, -1] 

            if (gradientAngle <= 22.5) { 
                //0
                neighbor1X = x + 1; neighbor1Y = y;
                neighbor2X = x - 1; neighbor2Y = y;
            }
            else if (gradientAngle <= 67.5) {
                //45
                neighbor1X = x + 1; neighbor1Y = y + 1;
                neighbor2X = x - 1; neighbor2Y = y - 1;
            }
            else if (gradientAngle <= 112.5) {
                //90
                neighbor1X = x; neighbor1Y = y + 1;
                neighbor2X = x; neighbor2Y = y - 1;
            }
            else if (gradientAngle <= 157.5) {
                //135
                neighbor1X = x - 1; neighbor1Y = y + 1;
                neighbor2X = x + 1; neighbor2Y = y - 1;
            }
            else if (gradientAngle <= 202.5) {
                //180 (same as <= 22.5)
                neighbor1X = x + 1; neighbor1Y = y;
                neighbor2X = x - 1; neighbor2Y = y;
            }
            else {
                //gradient angle is nan - 0 / 0
                //no edge
                neighbor1X = x; neighbor1Y = y;
                neighbor2X = x; neighbor2Y = y;                
            }
            
            const float* magVal = mag.GetPixelStart(x, y);
            float* resVal = res.GetPixelStart(x, y);

            if (magVal[0] < mag.GetPixelStart(neighbor1X, neighbor1Y)[0])
            {
                resVal[0] = 0;
                continue;                
            }
            
            if (magVal[0] < mag.GetPixelStart(neighbor2X, neighbor2Y)[0])
            {
                resVal[0] = 0;
                continue;                
            }
        }
    }
    return res;
}

template <typename T>
bool CannyEdgeDetector<T>::IsStrong(float val){
    return val > high;
}

template <typename T>
bool CannyEdgeDetector<T>::IsWeak(float val){
    return val < low;
}

template <typename T>
void CannyEdgeDetector<T>::Hysteresis(Image2d<float>& nonMaxSup){
    Image2d<T> result = nonMaxSup.CreateEmpty<T>();
  
    NeighborhoodKernel n3 = NeighborhoodKernel::CreateFromSize(3);

    for (size_t i = 0; i < nonMaxSup.GetPixelsCount(); i++){
        if (IsStrong(nonMaxSup[i][0]) == false){
            continue;
        }
       
        std::stack<size_t> st;
        nonMaxSup[i][0] = std::numeric_limits<float>::lowest();
        st.push(i);
       
        while (st.empty() == false)
        {
            size_t ii = st.top();
            st.pop();
            
            int x, y;
            nonMaxSup.GetPositionFromIndex(ii, x, y);

            nonMaxSup.ForEachPixelInNeighborhood(x, y, n3, [&](float* val, int xx, int yy){
                if (IsWeak(val[0])){
                    return;
                }

                size_t neighIndex = nonMaxSup.GetIndexFromPosition(xx, yy);
                st.push(neighIndex);
                val[0] = std::numeric_limits<float>::lowest();
            });

            result.GetPixelStart(x, y)[0] = ColorSpace::GetMaxGrayValue<T>();            
        }
    }
}

My magnitude result is this, which looks correct: