Question: I'm implementing MATLAB's imreconstruct in both Python and C++. However, for my test case, the Python implementation matches the MATLAB's output, the C++ one doesn't. This is the Python implementation: def imReconstruct(marker: np.array, mask: np.array) -> np.array: """ Naive implementation of MatLAB's imReconstruct function works when `mask` consists of mostly background (global minumum) will be slow otherwise """ kernel = cv.getStructuringElement(cv.MORPH_RECT, (3, 3)) # calculate the extreme values of the mask image min_val, max_val, _, _ = cv.minMaxLoc(mask) # clip the marker by global extrema of mask _, marker = cv.threshold(marker, min_val, max_val, cv.THRESH_TRUNC | cv.THRESH_BINARY_INV) while True: expanded = cv.dilate(marker, kernel) expanded = np.minimum(expanded, mask) # return `expanded` when the difference is small if np.max(np.abs(expanded - marker)) < 1e-5: return expanded # set expanded to marker and repeat marker = expanded D = np.array( [[4.2426405, 3.6055512, 3.1622777, 3. , 3. ], [3.6055512, 2.828427 , 2.236068 , 2. , 2. ], [3.1622777, 2.236068 , 1.4142135, 1. , 1. ], [3. , 2. , 1. , 0. , 0. ], [3. , 2. , 1. , 0. , 0. ]], dtype=float32) imReconstruct(D-.85,D) # output # array([[3.3926406, 3.3926406, 3.1622777, 3. , 3. ], # [3.3926406, 2.828427 , 2.236068 , 2. , 2. ], # [3.1622777, 2.236068 , 1.4142135, 1. , 1. ], # [3. , 2. , 1. , 0. , 0. ], # [3. , 2. , 1. , 0. , 0. ]], dtype=float32)  And C++: Mat imReconstruct(Mat marker, Mat mask){ /************** * naive implementation of MatLAB's imReconstruct * works when `mask` consist of mostly background (global minimum) * will be slow otherwise *************/ Mat kernel = getStructuringElement(MORPH_RECT, Size(3,3)); // calculate the min and max values from mask double minMask, maxMask; minMaxLoc(mask, &minMask, &maxMask); // clip the marker by global extrema of mask threshold(marker, marker, minMask, maxMask, THRESH_TRUNC|THRESH_BINARY_INV); Mat expanded; // keep filling the holes with `dilate` // until there are no more changes while (1){ dilate(marker, expanded, kernel); expanded = min(expanded, mask); // compute the max difference minMaxLoc(expanded-marker, &minMask, &maxMask); // return image when changes are small if (maxMask<1e-5) return expanded; // set expanded as marker and continue looping marker = expanded; } } // test case cv::Mat D = (cv::Mat_<float> (5,5) << 4.2426405, 3.6055512, 3.1622777, 3. , 3. , 3.6055512, 2.828427 , 2.236068 , 2. , 2. , 3.1622777, 2.236068 , 1.4142135, 1. , 1. , 3. , 2. , 1. , 0. , 0. , 3. , 2. , 1. , 0. , 0. ); std::cout << imReconstruct(D - .85, D) << std::endl; // output // [3.3926406, 3.3926406, 3.1622777, 2.7555513, 2.3122778; // 3.3926406, 2.8284271, 2.236068, 2, 2; // 3.1622777, 2.236068, 1.4142135, 1, 1; // 2.7555513, 2, 1, 0, 0; // 2.3122778, 2, 1, 0, 0]  What is the cause of this discrepancy? I may have overlooked something simple, but I already spent few hours in vain without any positive results.

Answer: There are minor, but significant differences between the Python and C++ implementation of the while loop body you presented. Those, coupled with some nuances of how cv::Mat and many OpenCV functions work came together to bite you. The variant of Python code that accurately matches the C++ implementation would look something like the following: expanded = None while True: expanded = cv.dilate(marker, kernel, expanded) expanded = cv.min(expanded, mask, expanded) _, max_val, _, _ = cv.minMaxLoc(expanded - marker) if max_val < 1e-5: return expanded marker = expanded  This implementation would then suffer from the same issue as the C++ variant, producing the same (incorrect result). But why? Let's step through the relevant part of C++ code (the while loop) and discuss what happens. Analysis We begin with cv::Mat marker containing a 5x5 float array, and cv::Mat expanded empty. First iteration We call cv::dilate(marker, expanded, kernel);. Internally, cv::Mat::create is called on the destination expanded. Since it's empty, a new 5x5 float array is allocated. Dilation is performed and result written to expanded's buffer. Next step is expanded = cv::min(expanded, mask);. This invokes an overload of cv::Mat::operator= which "can reuse already allocated matrix if it has the right size and type to fit the matrix expression result". Hence, this is equivalent to call to cv::min, specifically cv::min(expanded, mask, expanded);. Like before, cv::Mat::create is called on the destination expanded. This time it already contains a 5x5 float array, so no allocation happens, and the buffer is reused. The calculation is done in-place, and after the call, expanded still refers to the same chunk of memory. Next the cv::minMaxLoc runs on result of expanded - marker. Nothing surprising here, the two are different arrays. The maxMask is not less than 1e-5, so we continue. Problem Starts Trouble begins with the final statement in the loop body, marker = expanded; — a shallow copy. This invokes a different overload of cv::Mat::operator=. "Matrix assignment is an O(1) operation. This means that no data is copied but the data is shared and the reference counter, if any, is incremented." That means that after this point, both marker and expanded refer the the same 5x5 float array, with reference count of 2. With this in mind, we proceed to second iteration. Second Iteration We call cv::dilate(marker, expanded, kernel);. Destination is already 5x5 float array, no allocation happens, and result is written to expanded's buffer (which is still shared with marker). Next, expanded = cv::min(expanded, mask);. Destination is already 5x5 float array, no allocation happens, and result is written to expanded's buffer (which is still shared with marker). Failure Finally, the cv::minMaxLoc runs on result of expanded - marker... but both of those really refer to the same data, so it's rather expanded - expanded, and that's just all zeros. Therefore, the test (maxMask < 1e-5) succeeds and the loop terminates prematurely.  Solution 1 Ensure that new expanded is allocated on each iteration, as happens in the Python version of the code. Do this by simply moving declaration of expanded into the loop body. Code: // ... while (true) { cv::Mat expanded; cv::dilate(marker, expanded, kernel); expanded = cv::min(expanded, mask); cv::minMaxLoc(expanded - marker, nullptr, &maxMask); if (maxMask < 1e-5) { return expanded; } marker = expanded; }  Output: [3.3926406, 3.3926406, 3.1622777, 3, 3; 3.3926406, 2.8284271, 2.236068, 2, 2; 3.1622777, 2.236068, 1.4142135, 1, 1; 3, 2, 1, 0, 0; 3, 2, 1, 0, 0]  Solution 2 Instead of shallow copy, perform a deep copy. This is done by calling cv::Mat::copyTo. Code: Mat expanded; while (true) { cv::dilate(marker, expanded, kernel); expanded = cv::min(expanded, mask); cv::minMaxLoc(expanded - marker, nullptr, &maxMask); if (maxMask < 1e-5) { return expanded; } expanded.copyTo(marker); }  Output is identical to that of solution 1. NB: The Python equivalent of this would be np.copyto(expanded, marker).