Commit f5e8a972 authored by Tomaš Vank's avatar Tomaš Vank
Browse files

Upload New File

parent 50b74cb2
import matplotlib.pyplot as plt
import numpy as np
from skimage import color, draw, feature, filters, io, morphology
from skimage.color import rgb2gray
from skimage.draw import circle_perimeter
import scipy.ndimage as ndi
from skimage.measure import regionprops,label
def maximum_Filter(img):
threshold = np.max(img) / 1.5
img_max = ndi.maximum_filter1d(img, size=2, axis=0,mode='nearest')
img_max = ndi.maximum_filter1d(img_max, size=2, axis=1,mode='nearest')
img_t = img > threshold
return img_t,img_max,threshold
def hough_Peaks(hspaces, radii, total_num_peaks):
radius, cx, cy, acc = [], [], [], []
for rad, hough_peaks in zip(radii, hspaces):
img = hough_peaks.copy()
img_t, img_max, threshold = maximum_Filter(img)
img_peaks, y_peaks, x_peaks = [], [], []
label_img = label(img_t)
pointsWithMaxValues = regionprops(label_img, img_max)
pointsWithMaxValues = np.flipud(sorted(pointsWithMaxValues, key=lambda x: x.max_intensity)) # s najdenych max hodnot zo sortujeme od zaciatku
coords = np.array([np.round(p.centroid) for p in pointsWithMaxValues], dtype=int) # zo sortovanych bodov najdeme centroidi
for y_idx, x_idx in coords:
accum = img_max[y_idx, x_idx] # bod koordinatov s odfiltrovaneho obrazka
if accum > threshold:
img_peaks.append(accum)
y_peaks.append(y_idx)
x_peaks.append(x_idx)
img_peaks = np.array(img_peaks)
y_peaks = np.array(y_peaks)
x_peaks = np.array(x_peaks)
radius[len(radius):] = (rad,) * len(img_peaks)
cx[len(cx):] = x_peaks
cy[len(cy):] = y_peaks
acc[len(acc):] = img_peaks
radius = np.array(radius)
cx = np.array(cx)
cy = np.array(cy)
acc = np.array(acc)
s = [x for x, y in sorted(enumerate(acc), key=lambda x: x[1])]
cx_sorted, cy_sorted, r_sorted = np.flipud(cx[s]), np.flipud(cy[s]), np.flipud(radius[s])
return (cx_sorted[:total_num_peaks],cy_sorted[:total_num_peaks],r_sorted[:total_num_peaks])
def hough(img, radius):
x, y = np.where(img == 1)#hodnoty ktore maju nenulove hodnoty
hough_space = np.zeros((radius.size, img.shape[0], img.shape[1]), dtype=np.double)
i = 0
for rad in radius:
circle_x, circle_y = draw.circle_perimeter(0, 0, rad) #kruznica okolo polomeru rad
num_circle_pixels = circle_x.size
if num_circle_pixels > 800: # kruh ktory ma viac ako 80 px zmensim na polovicu vzniknu nepresnoti ale z rývhli sa kod
num_circle_pixels = circle_x.size//2
increase = 1.0 / num_circle_pixels
for p in range(x.size):#velkost nenulovej casti
for c in range(num_circle_pixels):#velkost kruznice v pixeloch
tx,ty = circle_x[c] + x[p],circle_y[c] + y[p]#urcovanie bodov pre nakreslenie kruznice
if img.shape[0] > tx >= 0 and img.shape[1] > ty >= 0:
hough_space[i, tx, ty] += increase
i = i + 1
return hough_space
def findBigCircles(cy, cx, radii):
for center_y, center_x, radius in zip(cy, cx, radii):
already_there = False
circy, circx = circle_perimeter(center_y, center_x, radius,
shape=image.shape)
for x, y in middle_coor_Big:
if ((x - center_x) ** 2 + (y - center_y) ** 2) ** 0.5 < 100:
already_there = True
if not already_there:
image[circy, circx] = (255, 255, 255)
image[center_y-2:center_y+2, center_x:center_x + radius] = (255, 255, 255)
middle_coor_Big.append([center_x, center_y])
return center_x, center_y, radii
def findCircles(cy, cx, radii,BigCenter_x,BigCenter_y,BigRadii):
for center_y, center_x, radius in zip(cy, cx, radii):
already_there = False
circy, circx = circle_perimeter(center_y, center_x, radius,
shape=image.shape)
point = BigRadii ** 2 - ((BigCenter_x - center_x) ** 2 + (BigCenter_y - center_y) ** 2)
for x, y in middle_coor:
if ((x - center_x) ** 2 + (y - center_y) ** 2) ** 0.5 < 50:#kedy ze sa mi prekryvaju kruznice
already_there = True
if not already_there and point[0] > 0:
image[circy, circx] = (20, 254, 53)
image[center_y-2:center_y+2, center_x:center_x + radius] = (20, 254, 53)
middle_coor.append([center_x, center_y])
return len(middle_coor)
I = rgb2gray(io.imread("WAC_GL000.tif")[8200:9200, 15800:16800])
#plt.imshow(I, cmap="gray")
#plt.figure()
filtered_moon = filters.rank.mean(I, selem=morphology.selem.disk(80))/255
#plt.imshow(filtered_moon, cmap="gray")
#plt.figure()
corrected_moon = I - filtered_moon
#plt.imshow(corrected_moon, cmap="gray")
#plt.figure()
prepared_moon = morphology.erosion(corrected_moon, selem=morphology.selem.disk(3))
plt.title("prepared")
plt.imshow(prepared_moon, cmap="gray")
plt.figure()
canny_moon = feature.canny(prepared_moon, sigma=10.5,low_threshold=2)
plt.imshow(canny_moon, cmap="gray")
plt.title("canny")
plt.figure()
#hough hough_circle
hough_radii = np.arange(10, 15, 1)
hough_radii_Big = np.arange(354, 355, 1)
hough_res = hough(canny_moon, hough_radii)
hough_res_Big = hough(canny_moon, hough_radii_Big)
total_num_peaks = 20
cx, cy, radii = hough_Peaks(hough_res, hough_radii, total_num_peaks=total_num_peaks)
cxBig, cyBig, radiiBig = hough_Peaks(hough_res_Big, hough_radii_Big, total_num_peaks=3)
middle_coor = []
middle_coor_Big = []
image = color.gray2rgb(I)
BigCenter_x,BigCenter_y,BigRadii = findBigCircles(cxBig, cyBig, radiiBig)
for BigCircle in middle_coor_Big:
smallCircles = findCircles(cy, cx, radii, BigCenter_x, BigCenter_y, BigRadii)
if smallCircles >= 5:
middle_coor_Big.clear()
findBigCircles(cxBig, cyBig, radiiBig)
break
middle_coor.clear()
image = color.gray2rgb(I)
print("Velke kratery")
number = 1
labeled = []
for idx in range(len(middle_coor_Big)):
x = cxBig[idx] + 8200
y = cyBig[idx] + 15800
if [cxBig[idx],cyBig[idx]] not in labeled:
print(str(number) + " " + str(x) + " / " + str(y))
plt.text(cxBig[idx], cyBig[idx], str(number))
number = 1 + number
labeled.append([cxBig[idx], cyBig[idx]])
number = 1
labeled = []
print("Male kratery")
for center_y, center_x in zip(cy, cx):
already_there = False
for i, j in middle_coor:
if ((i - center_x) ** 2 + (j - center_y) ** 2) ** 0.5 ==0: # kedy ze sa mi prekryvaju kruznice
already_there = True
point = BigRadii ** 2 - ((BigCenter_x - center_x) ** 2 + (BigCenter_y - center_y) ** 2)
if point[0] > 0 and already_there and [center_x,center_y] not in labeled:
x = center_x + 8200
y = center_y + 15800
print(str(number) + " " + str(x) + " / " + str(y))
plt.text(center_x, center_y, str(number))
number = 1 + number
labeled.append([center_x, center_y])
plt.imshow(image, cmap="gray")
plt.show()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment