in modules/imgproc/src/canny.cpp [251:564]
void operator()() const
{
#if CV_SSE2
bool haveSSE2 = checkHardwareSupport(CV_CPU_SSE2);
#endif
const int type = src.type(), cn = CV_MAT_CN(type);
Mat dx, dy;
ptrdiff_t mapstep = src.cols + 2;
// In sobel transform we calculate ksize2 extra lines for the first and last rows of each slice
// because IPPDerivSobel expects only isolated ROIs, in contrast with the opencv version which
// uses the pixels outside of the ROI to form a border.
uchar ksize2 = aperture_size / 2;
if (boundaries.start == 0 && boundaries.end == src.rows)
{
Mat tempdx(boundaries.end - boundaries.start + 2, src.cols, CV_16SC(cn));
Mat tempdy(boundaries.end - boundaries.start + 2, src.cols, CV_16SC(cn));
memset(tempdx.ptr<short>(0), 0, cn * src.cols*sizeof(short));
memset(tempdy.ptr<short>(0), 0, cn * src.cols*sizeof(short));
memset(tempdx.ptr<short>(tempdx.rows - 1), 0, cn * src.cols*sizeof(short));
memset(tempdy.ptr<short>(tempdy.rows - 1), 0, cn * src.cols*sizeof(short));
Sobel(src, tempdx.rowRange(1, tempdx.rows - 1), CV_16S, 1, 0, aperture_size, 1, 0, BORDER_REPLICATE);
Sobel(src, tempdy.rowRange(1, tempdy.rows - 1), CV_16S, 0, 1, aperture_size, 1, 0, BORDER_REPLICATE);
dx = tempdx;
dy = tempdy;
}
else if (boundaries.start == 0)
{
Mat tempdx(boundaries.end - boundaries.start + 2 + ksize2, src.cols, CV_16SC(cn));
Mat tempdy(boundaries.end - boundaries.start + 2 + ksize2, src.cols, CV_16SC(cn));
memset(tempdx.ptr<short>(0), 0, cn * src.cols*sizeof(short));
memset(tempdy.ptr<short>(0), 0, cn * src.cols*sizeof(short));
Sobel(src.rowRange(boundaries.start, boundaries.end + 1 + ksize2), tempdx.rowRange(1, tempdx.rows),
CV_16S, 1, 0, aperture_size, 1, 0, BORDER_REPLICATE);
Sobel(src.rowRange(boundaries.start, boundaries.end + 1 + ksize2), tempdy.rowRange(1, tempdy.rows),
CV_16S, 0, 1, aperture_size, 1, 0, BORDER_REPLICATE);
dx = tempdx.rowRange(0, tempdx.rows - ksize2);
dy = tempdy.rowRange(0, tempdy.rows - ksize2);
}
else if (boundaries.end == src.rows)
{
Mat tempdx(boundaries.end - boundaries.start + 2 + ksize2, src.cols, CV_16SC(cn));
Mat tempdy(boundaries.end - boundaries.start + 2 + ksize2, src.cols, CV_16SC(cn));
memset(tempdx.ptr<short>(tempdx.rows - 1), 0, cn * src.cols*sizeof(short));
memset(tempdy.ptr<short>(tempdy.rows - 1), 0, cn * src.cols*sizeof(short));
Sobel(src.rowRange(boundaries.start - 1 - ksize2, boundaries.end), tempdx.rowRange(0, tempdx.rows - 1),
CV_16S, 1, 0, aperture_size, 1, 0, BORDER_REPLICATE);
Sobel(src.rowRange(boundaries.start - 1 - ksize2, boundaries.end), tempdy.rowRange(0, tempdy.rows - 1),
CV_16S, 0, 1, aperture_size, 1, 0, BORDER_REPLICATE);
dx = tempdx.rowRange(ksize2, tempdx.rows);
dy = tempdy.rowRange(ksize2, tempdy.rows);
}
else
{
Mat tempdx(boundaries.end - boundaries.start + 2 + 2*ksize2, src.cols, CV_16SC(cn));
Mat tempdy(boundaries.end - boundaries.start + 2 + 2*ksize2, src.cols, CV_16SC(cn));
Sobel(src.rowRange(boundaries.start - 1 - ksize2, boundaries.end + 1 + ksize2), tempdx,
CV_16S, 1, 0, aperture_size, 1, 0, BORDER_REPLICATE);
Sobel(src.rowRange(boundaries.start - 1 - ksize2, boundaries.end + 1 + ksize2), tempdy,
CV_16S, 0, 1, aperture_size, 1, 0, BORDER_REPLICATE);
dx = tempdx.rowRange(ksize2, tempdx.rows - ksize2);
dy = tempdy.rowRange(ksize2, tempdy.rows - ksize2);
}
int maxsize = std::max(1 << 10, src.cols * (boundaries.end - boundaries.start) / 10);
std::vector<uchar*> stack(maxsize);
uchar **stack_top = &stack[0];
uchar **stack_bottom = &stack[0];
AutoBuffer<uchar> buffer(cn * mapstep * 3 * sizeof(int));
int* mag_buf[3];
mag_buf[0] = (int*)(uchar*)buffer;
mag_buf[1] = mag_buf[0] + mapstep*cn;
mag_buf[2] = mag_buf[1] + mapstep*cn;
// calculate magnitude and angle of gradient, perform non-maxima suppression.
// fill the map with one of the following values:
// 0 - the pixel might belong to an edge
// 1 - the pixel can not belong to an edge
// 2 - the pixel does belong to an edge
for (int i = boundaries.start - 1; i <= boundaries.end; i++)
{
int* _norm = mag_buf[(i > boundaries.start) - (i == boundaries.start - 1) + 1] + 1;
short* _dx = dx.ptr<short>(i - boundaries.start + 1);
short* _dy = dy.ptr<short>(i - boundaries.start + 1);
if (!L2gradient)
{
int j = 0, width = src.cols * cn;
#if CV_SSE2
if (haveSSE2)
{
__m128i v_zero = _mm_setzero_si128();
for ( ; j <= width - 8; j += 8)
{
__m128i v_dx = _mm_loadu_si128((const __m128i *)(_dx + j));
__m128i v_dy = _mm_loadu_si128((const __m128i *)(_dy + j));
v_dx = _mm_max_epi16(v_dx, _mm_sub_epi16(v_zero, v_dx));
v_dy = _mm_max_epi16(v_dy, _mm_sub_epi16(v_zero, v_dy));
__m128i v_norm = _mm_add_epi32(_mm_unpacklo_epi16(v_dx, v_zero), _mm_unpacklo_epi16(v_dy, v_zero));
_mm_storeu_si128((__m128i *)(_norm + j), v_norm);
v_norm = _mm_add_epi32(_mm_unpackhi_epi16(v_dx, v_zero), _mm_unpackhi_epi16(v_dy, v_zero));
_mm_storeu_si128((__m128i *)(_norm + j + 4), v_norm);
}
}
#elif CV_NEON
for ( ; j <= width - 8; j += 8)
{
int16x8_t v_dx = vld1q_s16(_dx + j), v_dy = vld1q_s16(_dy + j);
vst1q_s32(_norm + j, vaddq_s32(vabsq_s32(vmovl_s16(vget_low_s16(v_dx))),
vabsq_s32(vmovl_s16(vget_low_s16(v_dy)))));
vst1q_s32(_norm + j + 4, vaddq_s32(vabsq_s32(vmovl_s16(vget_high_s16(v_dx))),
vabsq_s32(vmovl_s16(vget_high_s16(v_dy)))));
}
#endif
for ( ; j < width; ++j)
_norm[j] = std::abs(int(_dx[j])) + std::abs(int(_dy[j]));
}
else
{
int j = 0, width = src.cols * cn;
#if CV_SSE2
if (haveSSE2)
{
for ( ; j <= width - 8; j += 8)
{
__m128i v_dx = _mm_loadu_si128((const __m128i *)(_dx + j));
__m128i v_dy = _mm_loadu_si128((const __m128i *)(_dy + j));
__m128i v_dx_ml = _mm_mullo_epi16(v_dx, v_dx), v_dx_mh = _mm_mulhi_epi16(v_dx, v_dx);
__m128i v_dy_ml = _mm_mullo_epi16(v_dy, v_dy), v_dy_mh = _mm_mulhi_epi16(v_dy, v_dy);
__m128i v_norm = _mm_add_epi32(_mm_unpacklo_epi16(v_dx_ml, v_dx_mh), _mm_unpacklo_epi16(v_dy_ml, v_dy_mh));
_mm_storeu_si128((__m128i *)(_norm + j), v_norm);
v_norm = _mm_add_epi32(_mm_unpackhi_epi16(v_dx_ml, v_dx_mh), _mm_unpackhi_epi16(v_dy_ml, v_dy_mh));
_mm_storeu_si128((__m128i *)(_norm + j + 4), v_norm);
}
}
#elif CV_NEON
for ( ; j <= width - 8; j += 8)
{
int16x8_t v_dx = vld1q_s16(_dx + j), v_dy = vld1q_s16(_dy + j);
int16x4_t v_dxp = vget_low_s16(v_dx), v_dyp = vget_low_s16(v_dy);
int32x4_t v_dst = vmlal_s16(vmull_s16(v_dxp, v_dxp), v_dyp, v_dyp);
vst1q_s32(_norm + j, v_dst);
v_dxp = vget_high_s16(v_dx), v_dyp = vget_high_s16(v_dy);
v_dst = vmlal_s16(vmull_s16(v_dxp, v_dxp), v_dyp, v_dyp);
vst1q_s32(_norm + j + 4, v_dst);
}
#endif
for ( ; j < width; ++j)
_norm[j] = int(_dx[j])*_dx[j] + int(_dy[j])*_dy[j];
}
if (cn > 1)
{
for(int j = 0, jn = 0; j < src.cols; ++j, jn += cn)
{
int maxIdx = jn;
for(int k = 1; k < cn; ++k)
if(_norm[jn + k] > _norm[maxIdx]) maxIdx = jn + k;
_norm[j] = _norm[maxIdx];
_dx[j] = _dx[maxIdx];
_dy[j] = _dy[maxIdx];
}
}
_norm[-1] = _norm[src.cols] = 0;
// at the very beginning we do not have a complete ring
// buffer of 3 magnitude rows for non-maxima suppression
if (i <= boundaries.start)
continue;
uchar* _map = map + mapstep*i + 1;
_map[-1] = _map[src.cols] = 1;
int* _mag = mag_buf[1] + 1; // take the central row
ptrdiff_t magstep1 = mag_buf[2] - mag_buf[1];
ptrdiff_t magstep2 = mag_buf[0] - mag_buf[1];
const short* _x = dx.ptr<short>(i - boundaries.start);
const short* _y = dy.ptr<short>(i - boundaries.start);
if ((stack_top - stack_bottom) + src.cols > maxsize)
{
int sz = (int)(stack_top - stack_bottom);
maxsize = std::max(maxsize * 3/2, sz + src.cols);
stack.resize(maxsize);
stack_bottom = &stack[0];
stack_top = stack_bottom + sz;
}
#define CANNY_PUSH(d) *(d) = uchar(2), *stack_top++ = (d)
#define CANNY_POP(d) (d) = *--stack_top
int prev_flag = 0;
bool canny_push = false;
for (int j = 0; j < src.cols; j++)
{
#define CANNY_SHIFT 15
const int TG22 = (int)(0.4142135623730950488016887242097*(1<<CANNY_SHIFT) + 0.5);
int m = _mag[j];
if (m > low)
{
int xs = _x[j];
int ys = _y[j];
int x = std::abs(xs);
int y = std::abs(ys) << CANNY_SHIFT;
int tg22x = x * TG22;
if (y < tg22x)
{
if (m > _mag[j-1] && m >= _mag[j+1]) canny_push = true;
}
else
{
int tg67x = tg22x + (x << (CANNY_SHIFT+1));
if (y > tg67x)
{
if (m > _mag[j+magstep2] && m >= _mag[j+magstep1]) canny_push = true;
}
else
{
int s = (xs ^ ys) < 0 ? -1 : 1;
if (m > _mag[j+magstep2-s] && m > _mag[j+magstep1+s]) canny_push = true;
}
}
}
if (!canny_push)
{
prev_flag = 0;
_map[j] = uchar(1);
continue;
}
else
{
// _map[j-mapstep] is short-circuited at the start because previous thread is
// responsible for initializing it.
if (!prev_flag && m > high && (i <= boundaries.start+1 || _map[j-mapstep] != 2) )
{
CANNY_PUSH(_map + j);
prev_flag = 1;
}
else
_map[j] = 0;
canny_push = false;
}
}
// scroll the ring buffer
_mag = mag_buf[0];
mag_buf[0] = mag_buf[1];
mag_buf[1] = mag_buf[2];
mag_buf[2] = _mag;
}
// now track the edges (hysteresis thresholding)
while (stack_top > stack_bottom)
{
if ((stack_top - stack_bottom) + 8 > maxsize)
{
int sz = (int)(stack_top - stack_bottom);
maxsize = maxsize * 3/2;
stack.resize(maxsize);
stack_bottom = &stack[0];
stack_top = stack_bottom + sz;
}
uchar* m;
CANNY_POP(m);
// Stops thresholding from expanding to other slices by sending pixels in the borders of each
// slice in a queue to be serially processed later.
if ( (m < map + (boundaries.start + 2) * mapstep) || (m >= map + boundaries.end * mapstep) )
{
borderPeaks.push(m);
continue;
}
if (!m[-1]) CANNY_PUSH(m - 1);
if (!m[1]) CANNY_PUSH(m + 1);
if (!m[-mapstep-1]) CANNY_PUSH(m - mapstep - 1);
if (!m[-mapstep]) CANNY_PUSH(m - mapstep);
if (!m[-mapstep+1]) CANNY_PUSH(m - mapstep + 1);
if (!m[mapstep-1]) CANNY_PUSH(m + mapstep - 1);
if (!m[mapstep]) CANNY_PUSH(m + mapstep);
if (!m[mapstep+1]) CANNY_PUSH(m + mapstep + 1);
}
}