Извлечение размеров прямоугольника

Question

Итак, я работаю над проектом обработки изображений с помощью сонарных изображений. Чтобы быть более точным, я пытаюсь извлечь размеры изображения пула, снятого сонарным сканером. Мне удалось извлечь прямоугольную область пула, но я не могу понять, как получить размеры каждого края в пикселях. Я работаю с OpenCV в Python.

Буду признателен, если кто-нибудь подскажет, как это можно сделать.

Я уже пытался найти пересечения линий грубых линий, но это не дало многообещающих результатов.

Код пока.

import cv2 
import numpy as np 
from scipy import ndimage as ndi
from scipy.ndimage.measurements import label

def largest_component(indices):
    #this function takes a list of indices denoting
    #the white regions of the image and returns the largest 
    #white object of connected indices 

    return_arr = np.zeros((512,512), dtype=np.uint8)
    for index in indeces:
        return_arr[index[0]][index[1]] = 255

    return return_arr

image = cv2.imread('sonar_dataset/useful/sonarXY_5.bmp', 0)
image_gaussian = ndi.gaussian_filter(image, 4)
image_gaussian_inv = cv2.bitwise_not(image_gaussian)
kernel = np.ones((3,3),np.uint8)

# double thresholding extracting the sides of the rectangle
ret1, img1 = cv2.threshold(image_gaussian_inv, 120, 255, cv2.THRESH_BINARY)
ret2, img2 = cv2.threshold(image_gaussian_inv, 150, 255, cv2.THRESH_BINARY)

double_threshold = img1 - img2
closing = cv2.morphologyEx(double_threshold, cv2.MORPH_CLOSE, kernel1)

labeled, ncomponents = label(closing, kernel)
indices = np.indices(closing.shape).T[:,:,[1, 0]]
twos = indices[labeled == 2]
area =[np.sum(labeled==val) for val in range(ncomponents+1)]

rectangle = largest_component(twos)

cv2.imshow('rectangle', rectangle)
cv2.waitKey(0)

Исходное изображение и извлеченный объект находятся ниже.

Answer 1

Итак, вот что я придумала - это немного трудоемко, но в конечном итоге оно дает нам правильный ответ. Я буду напрямую использовать вывод подключенных компонентов, который вы показали с последним изображением.

1. Используйте морфологическое скелетирование изображения , чтобы мы получили скелет капли. Таким образом, это даст нам самое минимальное представление контура, так что мы получим границу шириной в один пиксель, которая проходит через середину каждого толстого края. Этого можно добиться с помощью метода skeletonize Scikit-image.

2. Используйте Преобразование Хафа , которое является методом обнаружения линий на скелетонизированном изображении. Таким образом, он параметризует линии в полярной области, и на выходе будет набор rho и theta, которые сообщают нам, какие линии обнаружены на скелетонизированном изображении. Для этого мы можем использовать OpenCV cv2.HoughLines. Очень важно, чтобы вы делали это на скелетонизированном изображении, иначе у нас будет много линий-кандидатов, параллельных истинному контуру ограничивающего прямоугольника, и вы не сможете различить их.

3. Возьмите каждую пару линий и найдите их точку пересечения. Мы ожидаем, что со всеми парами линий будет 4 преобладающих группы пересечений, которые дают нам угол каждого прямоугольника.

4. Из-за шума в контурах мы можем получить более четырех точек пересечения. Мы можем использовать выпуклую оболочку , чтобы наконец получить 4 точки пересечения для прямоугольника. Таким образом, алгоритм выпуклой оболочки работает со списком точек, где он определяет подмножество точек, которые могут минимально охватывать список точек. Мы можем использовать cv2.convexHull.

5. Наконец, из-за квантования преобразования Хафа, может быть несколько точек, которые находятся в непосредственной близости от каждого угла. Поэтому примените K-Means кластеризацию , чтобы найти 4 группы точек и, таким образом, найти их центроиды. Для этого мы можем использовать cv2.kmeans.

6. Как только мы найдем центроиды, мы можем просто циклически перебирать каждую пару точек, чтобы, наконец, найти расстояния до каждого угла и, таким образом, найти расстояния, которые вас интересуют.

Давайте пройдемся по каждой точке по очереди:

Шаг № 1 - скелетизация морфологического изображения

Используя Scikit-image skeletonize, мы можем скелетировать изображение подключенных компонентов, которое вы показали выше. Обратите внимание, что вам нужно преобразовать изображение в двоичный файл, прежде чем продолжить. После того, как вы вызовете метод, нам нужно будет преобразовать обратно в 8-разрядное целое число без знака после для остальной части процесса. Я скачал изображение выше и сохранил его локально. Мы можем запустить метод skeletonize после:

from skimage.morphology import skeletonize

im = cv2.imread('K7ELI.png', 0)
out = skeletonize(im > 0)

# Convert to uint8
out = 255*(out.astype(np.uint8))

Мы получаем это изображение:

Шаг № 2 - Используйте преобразование Хафа

Используя преобразование Хафа, мы можем обнаружить наиболее заметные линии на этом изображении:

lines = cv2.HoughLines(out,1,np.pi/180,60)

Здесь мы указываем пространство поиска, так что мы ищем строки, где размер ячейки имеет длину 1, а углы имеют ячейку 1 градус, или pi / 180 радиан. Таким образом, преобразование Хафа просматривает каждую точку края и перебирает диапазон углов theta, которые проходят от начала координат до каждой точки края, и вычисляет соответствующее значение rho с учетом размера ячейки. Эта пара регистрируется в 2D гистограмме, и мы регистрируем голосование. Мы порождаем эту двумерную гистограмму так, чтобы любые ячейки за пределами определенного значения были кандидатами в строки. В приведенной выше строке кода установите пороговое значение для количества бинов равное 60.

Этот код не является обязательным, но я хотел показать вам, как выглядят визуализированные строки:

img_colour = np.dstack([im, im, im])
lines = cv2.HoughLines(edges,1,np.pi/180,60)
for rho,theta in lines[:,0]:
    a = np.cos(theta)
    b = np.sin(theta)
    x0 = a*rho
    y0 = b*rho
    x1 = int(x0 + 1000*(-b))
    y1 = int(y0 + 1000*(a))
    x2 = int(x0 - 1000*(-b))
    y2 = int(y0 - 1000*(a))
    cv2.line(img_colour,(x1,y1),(x2,y2),(0,0,255),2)

Этот код я вытащил из следующего урока . Он рисует обнаруженные линии преобразования Хафа на изображении красным цветом. Я получаю следующее изображение:

Как мы видим, на изображении четыре точки пересечения.Наша задача - найти эти точки пересечения.

Шаг № 3 - Найти точки пересечения

В преобразовании Хафа мы можем связать длину линии от начала координат до точки(x, y) в изображении, представленном под углом theta с помощью:

rho = x*cos(theta) + y*sin(theta)

Мы также можем сформировать уравнение прямой y = m*x + c в декартовой форме.Мы можем преобразовать между ними, разделив обе стороны уравнения rho на sin(theta), а затем переместив соответствующие термины в каждую сторону:

Следовательно, мы должны выполнить циклвсе уникальные пары линий и, используя приведенное выше уравнение, мы можем найти их точку пересечения, установив их декартовы формы равными друг другу.Это я не буду выводить для вас в целях экономии места, а просто установите две линии в декартовой форме равными друг другу и решите для x координаты пересечения.Как только это будет сделано, подставьте эту точку в любую из двух строк, чтобы найти координату y.Очевидно, что мы должны пропустить точки пересечения, которые выходят за пределы изображения в случае двух почти параллельных линий или если мы выбираем две пары линий, которые идут в одном направлении и не пересекаются.

pts = []
for i in range(lines.shape[0]):
    (rho1, theta1) = lines[i,0]
    m1 = -1/np.tan(theta1)
    c1 = rho1 / np.sin(theta1)
    for j in range(i+1,lines.shape[0]):
        (rho2, theta2) = lines[j,0]
        m2 = -1 / np.tan(theta2)
        c2 = rho2 / np.sin(theta2)
        if np.abs(m1 - m2) <= 1e-8:
            continue
        x = (c2 - c1) / (m1 - m2)
        y = m1*x + c1
        if 0 <= x < img.shape[1] and 0 <= y < img.shape[0]:
            pts.append((int(x), int(y)))

pts - это список кортежей, в который мы добавляем все точки пересечения, которые находятся внутри изображения, которые не выходят за границы.

Шаг # 4 - Использование выпуклой оболочки

Мы можем использоватьэтот список кортежей и использовать выпуклую оболочку, чтобы мы нашли список точек, которые определяют внешний периметр прямоугольника.Обратите внимание, что порядок точек, определяющих прямоугольник, идет против часовой стрелки.Это не имеет значения для этого шага, но будет иметь значение позже:

pts = np.array(pts)
pts = pts[:,None] # We need to convert to a 3D numpy array with a singleton 2nd dimension
hull = cv2.convexHull(pts)

hull содержит массив 3D NumPy, который является подмножеством исходных точек пересечения, которые создают внешнюю границу изображения.Мы можем использовать эти точки, чтобы нарисовать, где они расположены на изображении для иллюстрации

out2 = np.dstack([im, im, im])
for pt in hull[:,0]:
    cv2.circle(out2, tuple(pt), 2, (0, 255, 0), 2)

Я взял исходное изображение и нарисовал угловые точки зеленым цветом.Мы получаем это изображение:

Шаг № 5 - Применение кластеризации K-Means

Как вы можете видеть на изображении выше,Есть несколько точек, которые отображаются в каждом углу.Было бы хорошо, если бы мы могли объединить несколько точек на каждом углу в одну точку.Одним из способов является усреднение всех точек в каждом углу, и самый простой способ сделать это из коробки - использовать кластеризацию K-Means.Нам нужны центроиды, чтобы, таким образом, дать нам последние угловые точки прямоугольника.Нам нужно убедиться, что мы указываем 4 кластера для поиска.

Из учебника по кластеризации K-Means из документации OpenCV мы можем использовать этот код:

# Define criteria = ( type, max_iter = 10 , epsilon = 1.0 )
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)

# Set flags (Just to avoid line break in the code)
flags = cv2.KMEANS_RANDOM_CENTERS

# Apply KMeans
# The convex hull points need to be float32
z = hull.copy().astype(np.float32)
compactness,labels,centers = cv2.kmeans(z,4,None,criteria,10,flags)

Первый параметр - это выпуклая оболочка точек, которые должны быть в float32, как того требует алгоритм.Второй параметр указывает количество кластеров, которые мы хотим найти, поэтому 4 в нашем случае.Третий параметр вы можете пропустить.Это заполнитель для лучшего идентификатора кластера, которому назначена каждая точка, но нам не нужно его использовать.criteria - это параметры K-средних, используемые для механики алгоритма, а пятый параметр говорит нам, сколько попыток нам нужно предпринять, чтобы найти лучшие кластеры.Мы выбираем 10, что означает, что мы запускаем K-Means 10 раз и выбираем конфигурацию кластеризации, которая имеет наименьшее количество ошибок.Ошибка сохраняется в переменной compactness, которая выводится из алгоритма.Наконец, последняя переменная - необязательные флаги, и мы установили ее так, чтобы начальные центроиды алгоритма просто выбирались случайным образом из точек.

labels указывает, какой идентификатор кластера назначен для каждой точки, а centers - это ключевая переменная, которая нам нужна и которая возвращает:

array([[338.5   , 152.5   ],
       [302.6667, 368.6667],
       [139.    , 340.    ],
       [178.5   , 127.    ]], dtype=float32)

Это четыре угловые точки прямоугольника.,Мы можем видеть, где они располагаются, рисуя их непосредственно на исходном изображении, и мы также получаем это изображение:

out3 = np.dstack([im, im, im])
for pt in centers:
    cv2.circle(out3, tuple(pt), 2, (0, 255, 0), 2)

Шаг № 6 - Измерьте длину сейчас

Наконец, мы можем перебрать каждую пару линий и найти соответствующие измерения. Обратите внимание, что поскольку K-Means имеет центроиды в случайном порядке из-за случайной природы алгоритма, мы можем запустить выпуклую оболочку на этих центроидах, чтобы убедиться, что порядок круговой.

centers = cv2.convexHull(centers)[:,0]
for (i, j) in zip(range(4), [1, 2, 3, 0]):
    length = np.sqrt(np.sum((centers[i] - centers[j])**2.0))
    print('Length of side {}: {}'.format(i+1, length))

Таким образом, мы получаем:

Length of side 1: 219.11654663085938
Length of side 2: 166.1582489013672
Length of side 3: 216.63160705566406
Length of side 4: 162.019287109375

Если вы хотите, чтобы перспектива видела, как выровняется ограничивающая рамка, давайте нарисуем эти линии на изображении, которые определены в этих центрах:

out4 = np.dstack([im, im, im])
for (i, j) in zip(range(4), [1, 2, 3, 0]):
    cv2.line(out4, tuple(centers[i]), tuple(centers[j]), (0, 0, 255), 2)

Получаем:

Чтобы увидеть, где это совпадает с исходным изображением, давайте просто повторим приведенный выше код, но рисуем линии на исходном изображении. Я скачал копию исходного изображения, чтобы сделать это:

out5 = cv2.imread('no8BP.png') # Note - grayscale image read in as colour
for (i, j) in zip(range(4), [1, 2, 3, 0]):
    cv2.line(out5, tuple(centers[i]), tuple(centers[j]), (0, 0, 255), 2)

---

Для полноты, вот весь код от начала до конца без всех выходных данных отладки - мы переходим от чтения изображения к рисованию линий в исходном изображении с печатью длин каждой стороны в обнаруженном прямоугольнике.

from skimage.morphology import skeletonize
import cv2
import numpy as np

# Step #1 - Skeletonize
im = cv2.imread('K7ELI.png', 0)
out = skeletonize(im > 0)

# Convert to uint8
out = 255*(out.astype(np.uint8))

# Step #2 - Hough Transform
lines = cv2.HoughLines(out,1,np.pi/180,60)

# Step #3 - Find points of intersection
pts = []
for i in range(lines.shape[0]):
    (rho1, theta1) = lines[i,0]
    m1 = -1/np.tan(theta1)
    c1 = rho1 / np.sin(theta1)
    for j in range(i+1,lines.shape[0]):
        (rho2, theta2) = lines[j,0]
        m2 = -1 / np.tan(theta2)
        c2 = rho2 / np.sin(theta2)
        if np.abs(m1 - m2) <= 1e-8:
            continue
        x = (c2 - c1) / (m1 - m2)
        y = m1*x + c1
        if 0 <= x < img.shape[1] and 0 <= y < img.shape[0]:
            pts.append((int(x), int(y)))

# Step #4 - Find convex hull
pts = np.array(pts)
pts = pts[:,None] # We need to convert to a 3D numpy array with a singleton 2nd dimension
hull = cv2.convexHull(pts)

# Step #5 - K-Means clustering
# Define criteria = ( type, max_iter = 10 , epsilon = 1.0 )
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)

# Set flags (Just to avoid line break in the code)
flags = cv2.KMEANS_RANDOM_CENTERS

# Apply KMeans
# The convex hull points need to be float32
z = hull.copy().astype(np.float32)
compactness,labels,centers = cv2.kmeans(z,4,None,criteria,10,flags)

# Step #6 - Find the lengths of each side
centers = cv2.convexHull(centers)[:,0]
for (i, j) in zip(range(4), [1, 2, 3, 0]):
    length = np.sqrt(np.sum((centers[i] - centers[j])**2.0))
    print('Length of side {}: {}'.format(i+1, length))

# Draw the sides of each rectangle in the original image
out5 = cv2.imread('no8BP.png') # Note - grayscale image read in as colour
for (i, j) in zip(range(4), [1, 2, 3, 0]):
    cv2.line(out5, tuple(centers[i]), tuple(centers[j]), (0, 0, 255), 2)

# Show the image
cv2.imshow('Output', out5); cv2.waitKey(0); cv2.destroyAllWindows()

Answer 2

Уже есть два хороших решения, я хотел предложить более простое, основанное на другом способе обнаружения прямоугольника.

(здесь я использую MATLAB с DIPimage , потому чтодля меня быстрее собрать вместе концепцию, чем в Python, но точно такая же функциональность доступна в Python, см. в конце поста. Отказ от ответственности: я автор DIPimage.)

Потому чтопрямоугольник - это яркая фигура на темном фоне, и (я полагаю) гарантированно окружит центр изображения, мы можем создать семя в центре изображения и семя на периферии и использовать водораздел, чтобы найтипрямоугольникВ этом случае водораздел гарантируется для получения одного замкнутого контура толщиной 1 пиксель.

img = readim('https://i.stack.imgur.com/no8BP.png');
seeds = clone(img);
seeds(rr(seeds)<50) = 1;
seeds(rr(seeds)>250) = 2;
rect = waterseed(seeds,gaussf(img));
overlay(img,rect) % for display only

Обратите внимание, что я сглаживалнемного входное изображение.Но прямоугольник все еще довольно шумный, что повлияет на измерение размера, которое мы сделаем позже.Мы можем сгладить его, используя морфологическое отверстие с большим круговым структурирующим элементом.Эта операция также обрезает углы, но закругленные углы не повлияют на результат измерения.

rect = opening(fillholes(rect),35);
overlay(img,rect-berosion(rect)) % for display only

Теперь у нас хорошая формаэто подходит для измерения.Диаметры Ферета - это длины выступов формы.Мы измеряем длину самой короткой проекции (равной ширине прямоугольника) и длину проекции, перпендикулярной самой короткой (равной длине прямоугольника).См. мой пост в блоге для подробного описания алгоритма, который вычисляет эти длины.

msr = measure(rect,[],'feret');
sz = msr(1).feret(2:3)

Это возвращает sz = [162.7506, 215.0775].

---

ВотPython-эквивалент кода выше (запускаются точно такие же алгоритмы).PyDIP, привязки Python для библиотеки DIPlib, не так совершенен, как набор инструментов DIPimage, который я использую выше, и часть синтаксиса немного более многословна (хотя в основном специально).Коллега работает над упаковкой бинарного дистрибутива PyDIP, до этого вам нужно было бы собрать его из источников , что, как мы надеемся, довольно просто, если вы будете следовать инструкциям .

import PyDIP as dip

img = dip.ImageRead('no8BP.png')

seeds = img.Similar()
seeds.Fill(0)
rr = dip.CreateRadiusCoordinate(seeds.Sizes())
seeds[rr<50] = 1
seeds[rr>250] = 2

rect = dip.SeededWatershed(dip.Gauss(img), seeds)
dip.viewer.Show(dip.Overlay(img,rect))
dip.viewer.Spin()

rect = dip.Opening(dip.FillHoles(rect),35)
dip.viewer.Show(dip.Overlay(img,rect-dip.BinaryErosion(rect,1,1)))
dip.viewer.Spin()

msr = dip.MeasurementTool.Measure(dip.Label(rect),features=['Feret'])
sz = (msr[1]['Feret'][1],msr[1]['Feret'][2])
print(sz)

---

Скорее всего, это можно реализовать и в OpenCV, но это может быть немного сложнее.Например, две меры Ферета, которые мы здесь вычисляем, эквивалентны тому, что возвращает OpenCV minAreaRect, а отсеянный водораздел включен в OpenCV watershed.

Answer 3

Это не идеально, но этот простой подход должен быть хорошей отправной точкой для вас:

import cv2, math
import numpy as np

img = cv2.imread(R'D:\dev\projects\stackoverflow\dimensions_of_rectangle\img1.png')
print(img.shape)
img_moments=cv2.moments(img[:,:,0]) #use only one channel here (cv2.moments operates only on single channels images)
print(img_moments)
# print(dir(img_moments))

# calculate centroid (center of mass of image)
x = img_moments['m10'] / img_moments['m00']
y = img_moments['m01'] / img_moments['m00']

# calculate orientation of image intensity (it corresponds to the image intensity axis)
u00 = img_moments['m00']
u20 = img_moments['m20'] - x*img_moments['m10']
u02 = img_moments['m02'] - y*img_moments['m01']
u11 = img_moments['m11'] - x*img_moments['m01']

u20_prim = u20/u00
u02_prim = u02/u00
u11_prim = u11/u00

angle = 0.5 * math.atan(2*u11_prim / (u20_prim - u02_prim))
print('The image should be rotated by: ', math.degrees(angle) / 2.0, ' degrees')

cols,rows = img.shape[:2]
# rotate the image by half of this angle
rotation_matrix = cv2.getRotationMatrix2D((cols/2,rows/2), math.degrees(angle / 2.0), 1)
img_rotated = cv2.warpAffine(img, rotation_matrix ,(cols,rows))
# print(img_rotated.shape, img_rotated.dtype)

cv2.imwrite(R'D:\dev\projects\stackoverflow\dimensions_of_rectangle\img1_rotated.png', img_rotated)

img_rotated_clone = np.copy(img_rotated)
img_rotated_clone2 = np.copy(img_rotated)

# first method - just calculate bounding rect
bounding_rect = cv2.boundingRect(img_rotated[:, :, 0])
cv2.rectangle(img_rotated_clone, (bounding_rect[0], bounding_rect[1]), 
    (bounding_rect[0] + bounding_rect[2], bounding_rect[1] + bounding_rect[3]), (255,0,0), 2)

# second method - find columns and rows with biggest sums
def nlargest_cols(a, n):
    col_sums = [(np.sum(col), idx) for idx, col in enumerate(a.T)]
    return sorted(col_sums, key=lambda a: a[0])[-n:]

def nlargest_rows(a, n):
    col_sums = [(np.sum(col), idx) for idx, col in enumerate(a[:,])]
    return sorted(col_sums, key=lambda a: a[0])[-n:]

top15_cols_indices = nlargest_cols(img_rotated[:,:,0], 15)
top15_rows_indices = nlargest_rows(img_rotated[:,:,0], 15)

for a in top15_cols_indices:
    cv2.line(img_rotated_clone, (a[1], 0), (a[1], rows), (0, 255, 0), 1)

for a in top15_rows_indices:
    cv2.line(img_rotated_clone, (0, a[1]), (cols, a[1]), (0, 0, 255), 1)

cv2.imwrite(R'D:\dev\projects\stackoverflow\dimensions_of_rectangle\img2.png', img_rotated_clone)

Конечно, вам необходимо скорректировать пути.img1.png - это второе изображение из вашего вопроса, img1_rotated - результат поворота изображения:

, а img2 - окончательный результат:

Синий прямоугольник - это метод1 (просто ограничивающий прямоугольник), а зеленые и красные линии (15 красных и 15 зеленых - все шириной в 1 пиксель) - второй метод.

Алгоритм довольно прост:

1. Рассчитайте моменты изображения, чтобы определить основную ось интенсивности изображения (я не знаю, как это хорошо описать - проверьте страницу вики https://en.wikipedia.org/wiki/Image_moment#Examples_2).В основном это угол, на который вам нужно повернуть изображение, чтобы белые пиксели распределились по горизонтали или вертикали.
2. Как только вы узнаете угол - поверните изображение (и сохраните результат).
3. Метод 1 - вычислить и нарисовать повернутый прямоугольник всех пикселей.
4. Метод 2 - найти 15 строк и 15 столбцов с наибольшей суммой (== наибольшее количество белых пикселей) и нарисовать горизонтальные / вертикальные линии в этих строках / столбцах,Обратите внимание, что число 15 было выбрано методом проб и ошибок, но должно быть легко выбрать 2 столбца (и строки) с большой суммой, которые не близки друг к другу.Эти столбцы / строки являются хорошими кандидатами на границы прямоугольника.

Надеюсь, это то, что вы искали, дайте мне знать, что у вас возникнут вопросы.