import math
import numpy as np
import pandas as pd
from typing import Tuple
def find_baseline_carti(mask: np.ndarray) -> Tuple[Tuple[int, int], Tuple[int, int]]:
"""
Funkcja znajdująca baseline, który ma najwięcej punktów pokrywających się
z górnym obrysem bony roof. Algorytm iteruje po pikselach wysuniętych
najbardziej na prawo punktów bony roof, które są oddalone o 30 pikseli
w górę i 10 pikseli w dół w celu zmniejszenia ilości obliczeń.
Funkcja przyjmuje:
- maskę w formacie png z zaznaczonymi strukturami
Funkcja zwraca:
- współrzędne prawego i lewego punktu baseline
"""
height, width = mask.shape
label_upper_contour_bony = 5
label_upper_contour_carti = 4
binary_mask_bony = (mask == label_upper_contour_bony).astype(np.uint8)
binary_mask_carti = (mask == label_upper_contour_carti).astype(np.uint8)
upper_contour_bony = []
# Znalezienie punktów górnego obrysu bony roof
for x in range(width):
column_bony = binary_mask_bony[:, x]
if np.any(column_bony):
y_bony = np.where(column_bony)[0][0]
upper_contour_bony.append((x, y_bony))
# Znalezienie najbardziej na prawo wysuniętego punktu carti roof
x_max_carti = -1
for x in range(width):
column_carti = binary_mask_carti[:, x]
if np.any(column_carti):
x_max_carti = x # Ustalamy x_max_carti na najbardziej wysunięty na prawo punkt
y_max_carti = np.where(column_carti)[0][0] # Używamy tylko y, ale nie potrzebujemy
# Ustalanie prawego punktu bony roof (+1 piksel w prawo)
rightmost_point_x = x_max_carti + 1 if x_max_carti != -1 else None
# Znalezienie punktu bony roof dla współrzędnej x = x_max_carti + 1
bony_point_y = None
if rightmost_point_x is not None and rightmost_point_x < width:
column_bony = binary_mask_bony[:, rightmost_point_x]
if np.any(column_bony):
bony_point_y = np.where(column_bony)[0][0] # Znajdź punkt na bony roof w kolumnie x_max_carti + 1
# Definiowanie punktów startowych (iteracja 30 pikseli w górę i 10 w dół od punktu bony roof)
start_points = []
if bony_point_y is not None:
start_points = [(rightmost_point_x, y) for y in range(max(0, bony_point_y - 5),
min(height, bony_point_y + 5))]
best_line = []
max_overlap = 0
rightmost_point = (rightmost_point_x, bony_point_y) if rightmost_point_x is not None else (None, None)
for x0, y0 in start_points:
for angle in range(360):
line_points = []
angle_rad = math.radians(angle)
sin_angle = math.sin(angle_rad)
cos_angle = math.cos(angle_rad)
# Definiujemy maksymalną liczbę kroków
n_max = 2 * (width + height)
n = np.arange(n_max)
x = x0 + n * cos_angle
y = y0 - n * sin_angle
x_rounded = np.round(x).astype(int)
y_rounded = np.round(y).astype(int)
# Filtrujemy punkty znajdujące się w granicach obrazu
valid_mask = (x_rounded >= 0) & (x_rounded < width) & \
(y_rounded >= 0) & (y_rounded < height)
x_valid = x_rounded[valid_mask]
y_valid = y_rounded[valid_mask]
line_points = list(zip(x_valid, y_valid))
# Obliczenie overlapu między linią a konturem
overlap = len(set(line_points) & set((p for p in upper_contour_bony if p[0] <= rightmost_point_x)))
if overlap > max_overlap:
max_overlap = overlap
best_line = line_points
# Znalezienie najbardziej na lewo wysuniętego punktu linii
leftmost_point = None
if best_line:
for point in best_line:
x, y = point
if mask[y, x] == label_upper_contour_bony:
if leftmost_point is None or x < leftmost_point[0]:
leftmost_point = (x, y)
if leftmost_point is None:
leftmost_point = (None, None)
return rightmost_point, leftmost_point
def find_baseline(mask: np.ndarray) -> Tuple[Tuple[int, int], Tuple[int, int]]:
"""
Funkcja znajdująca baseline, który ma najwięcej punktów pokrywających się
z górnym obrysem bony roof. Algorytm iteruje po pikselach wysuniętych
najbardziej na prawo punktów bony roof, które są oddalone o 10 pikseli
w górę i w dół w celu zmniejszenia ilości obliczeń.
Funkcja przyjmuje:
- maskę w formacie png z zaznaczonymi strukturami
Funkcja zwraca:
- współrzędne prawego i lewego punktu baseline
"""
height, width = mask.shape
label_upper_contour = 5
binary_mask = (mask == label_upper_contour).astype(np.uint8)
upper_contour = []
# Znalezienie punktów górnego obrysu (konturu)
for x in range(width):
column = binary_mask[:, x]
if np.any(column):
y = np.where(column)[0][0]
upper_contour.append((x, y))
# Znalezienie najbardziej na prawo wysuniętego punktu
x_max = max(p[0] for p in upper_contour)
y_max = max(p[1] for p in upper_contour if p[0] == x_max)
# Iterowanie po pikselach znajdujących się 30 pikseli ponad i 10
# pikseli podnajbardziej wysuniętymi na prawo punktami bony roof
start_points = [(x, y) for x, y in upper_contour if x == x_max] +\
[(x_max, y) for y in range(max(0, y_max - 30),
min(height, y_max + 10))]
best_line = []
max_overlap = 0
rightmost_point = None
for x0, y0 in start_points:
for angle in range(360):
line_points = []
angle_rad = math.radians(angle)
sin_angle = math.sin(angle_rad)
cos_angle = math.cos(angle_rad)
# Definiujemy maksymalną liczbę kroków
n_max = 2 * (width + height)
n = np.arange(n_max)
x = x0 + n * cos_angle
y = y0 - n * sin_angle
x_rounded = np.round(x).astype(int)
y_rounded = np.round(y).astype(int)
# Filtrujemy punkty znajdujące się w granicach obrazu
valid_mask = (x_rounded >= 0) & (x_rounded < width) & \
(y_rounded >= 0) & (y_rounded < height)
x_valid = x_rounded[valid_mask]
y_valid = y_rounded[valid_mask]
line_points = list(zip(x_valid, y_valid))
# Obliczenie overlapu między linią a konturem
overlap = len(set(line_points) & set(upper_contour))
if overlap > max_overlap:
max_overlap = overlap
best_line = line_points
rightmost_point = (x0, y0)
# Znalezienie najbardziej na lewo wysuniętego punktu linii
leftmost_point = None
if best_line:
for point in best_line:
x, y = point
if mask[y, x] == label_upper_contour:
if leftmost_point is None or x < leftmost_point[0]:
leftmost_point = (x, y)
if leftmost_point is None:
leftmost_point = (None, None)
return rightmost_point, leftmost_point
def old_find_baseline(class_mask: np.ndarray) -> Tuple[
Tuple[int, int] | None, Tuple[int, int] | None]:
"""Bazując na położeniu bony i carti roof określa punkt początkowy
oraz końcowy baseline dla maksymalnej odległości do 5 pikseli.
Algorytm przeszukuje obszar jednego piksela dookoła (czyli 3x3),
a jeżeli nie znajdzie punktów, to zwiększa obszar poszukiwań aż
do progu jakim jest 5 pikseli. Dodatkowo na początku algorytm
zawęża obszar poszukiwań tylko do obszarów bony roof i carti
roof oraz pikselom im sąsiadującym
Args:
class_mask (np.ndarray): Maska wygenerowana z txt
Returns:
Współrzędne punktów baseline
"""
mask_4_indices = np.argwhere(class_mask == 4)
mask_5_indices = np.argwhere(class_mask == 5)
if mask_4_indices.size == 0 or mask_5_indices.size == 0:
return None, None
min_y = max(0, min(np.min(mask_4_indices[:, 0]), np.min(mask_5_indices[
:, 0])))
max_y = min(class_mask.shape[0], max(np.max(mask_4_indices[:, 0]), np.max(
mask_5_indices[:, 0])))
min_x = max(0, min(np.min(mask_4_indices[:, 1]), np.min(
mask_5_indices[:, 1])))
max_x = min(class_mask.shape[1], max(np.max(mask_4_indices[:, 1]), np.max(
mask_5_indices[:, 1])))
for distance in range(1, 6):
right_point = None
left_point = None
for y in range(min_y, max_y + 1):
for x in range(min_x, max_x + 1):
surroundings = class_mask[max(0, y - distance):min(
class_mask.shape[0], y + distance + 1),
max(0, x - distance):min(
class_mask.shape[1], x + distance + 1)]
if (surroundings == 5).any() and (surroundings == 4).any():
if right_point is None or (x > right_point[0] or (
x == right_point[0] and y > right_point[1])):
right_point = (x, y)
if left_point is None or (x < left_point[0] or (
x == left_point[0] and y > left_point[1])):
left_point = (x, y)
if right_point is not None and left_point is not None:
if distance == 5 and right_point == left_point:
left_point = (left_point[0] - 1, left_point[1])
if right_point == left_point:
continue
return right_point, left_point
return None, None
import os
import cv2
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
# Słownik kolorów dla poszczególnych klas
CLASS_COLORS = {
0: {'color_rgb': (0, 0, 0), 'label': 'background'},
1: {'color_rgb': (255, 0, 0), 'label': 'chondroosseous border'},
2: {'color_rgb': (0, 255, 0), 'label': 'femoral head'},
3: {'color_rgb': (0, 0, 255), 'label': 'labrum'},
4: {'color_rgb': (255, 255, 0), 'label': 'cartilagineous roof'},
5: {'color_rgb': (0, 255, 255), 'label': 'bony roof'},
6: {'color_rgb': (159, 2, 250), 'label': 'bony rim'},
7: {'color_rgb': (255, 132, 0), 'label': 'lower limb'},
8: {'color_rgb': (255, 0, 255), 'label': 'baseline'},
9: {'color_rgb': (66, 135, 245), 'label': 'lower limb template'},
10: {'color_rgb': (255, 69, 0), 'label': 'lower limb - alg. v3'}
}
def visualize_masks_with_baseline(masks_dir: str, output_directory: str):
# Upewnij się, że katalog wyjściowy istnieje
if not os.path.exists(output_directory):
os.makedirs(output_directory)
# Iterowanie po katalogu z maskami
for mask_filename in os.listdir(masks_dir):
if mask_filename.endswith('.png'):
mask_path = os.path.join(masks_dir, mask_filename)
mask = cv2.imread(mask_path, cv2.IMREAD_GRAYSCALE) # Wczytaj maskę w skali szarości
# Znajdź punkty baseline za pomocą nowego algorytmu
rightmost_point, leftmost_point = find_baseline_carti(mask)
# Debug: Sprawdzenie punktów
print(f"Maska: {mask_filename}, Punkt prawy: {rightmost_point}, Punkt lewy: {leftmost_point}")
# Stwórz obraz RGB, na który będziemy rysować
mask_rgb = np.zeros((mask.shape[0], mask.shape[1], 3), dtype=np.uint8)
# Iteracja po klasach w masce i kolorowanie
for class_value, class_info in CLASS_COLORS.items():
mask_rgb[mask == class_value] = class_info['color_rgb']
# Rysowanie białej linii między punktami zwróconymi przez find_baseline
if rightmost_point and leftmost_point and rightmost_point != (None, None) and leftmost_point != (None, None):
cv2.line(mask_rgb, rightmost_point, leftmost_point, (255, 255, 255), 1) # Biała linia
cv2.circle(mask_rgb, rightmost_point, 1, (255, 255, 255), -1) # Prawy punkt
cv2.circle(mask_rgb, leftmost_point, 1, (255, 255, 255), -1) # Lewy punkt
# Zapisz wizualizację do katalogu
output_path = os.path.join(output_directory, f'visualization_{mask_filename}')
cv2.imwrite(output_path, mask_rgb)
# Wyświetlenie obrazu
plt.figure(figsize=(10, 10))
plt.imshow(cv2.cvtColor(mask_rgb, cv2.COLOR_BGR2RGB)) # Zamiana na RGB dla wyświetlenia
plt.title(f'Visualization for {mask_filename}')
plt.axis('off')
# Przykładowe wywołanie funkcji
masks_directory = './Angles/dane/masks_baseline'
output_directory = './Angles/dane/visualizations_carti'
#visualize_masks_with_baseline(masks_directory, output_directory)
import os
import numpy as np
import matplotlib.pyplot as plt
from typing import Tuple, List
import math
def slope(x1: int, y1: int, x2: int, y2: int) -> float:
"""Calculates the slope of the line through two points."""
# Sprawdzanie czy punkty są None
if None in (x1, y1, x2, y2):
print(f"Invalid slope points: ({x1}, {y1}), ({x2}, {y2}). Skipping.")
return None # Zwracamy None, aby pominąć obliczenia
if x2 - x1 == 0: # Unikamy dzielenia przez 0
return float('inf') # Linia pionowa
return (y2 - y1) / (x2 - x1)
def find_angle(line1_point1: Tuple[int, int],
line1_point2: Tuple[int, int],
line2_point1: Tuple[int, int],
line2_point2: Tuple[int, int],
return_absolute_value: bool = True) -> float:
"""Finds the angle between two lines."""
if len(line1_point1) != 2 or len(line1_point2) != 2 or \
len(line2_point1) != 2 or len(line2_point2) != 2:
print('Point with wrong number of coordinates')
return np.NaN # Zwracamy NaN zamiast zgłaszać wyjątek
M1 = slope(line1_point1[0], line1_point1[1], line1_point2[0], line1_point2[1])
M2 = slope(line2_point1[0], line2_point1[1], line2_point2[0], line2_point2[1])
PI = 3.14159265
# Sprawdzamy, czy M1 i M2 nie są None
if M1 is None or M2 is None:
print(f"Invalid slopes: M1={M1}, M2={M2}. Skipping angle calculation.")
return None
if M2 != float('inf') and M1 * M2 != -1:
if return_absolute_value:
angle = abs((M2 - M1) / (1 + M1 * M2))
else:
angle = (M2 - M1) / (1 + M1 * M2)
ret = math.atan(angle)
val = (ret * 180) / PI
return round(val, 1)
elif M1 * M2 == -1:
return 90.0
else:
return np.NaN
def read_doctor_coordinates(file_path: str) -> Tuple[Tuple[int, int], Tuple[int, int]]:
"""Reads doctor coordinates from a text file."""
with open(file_path, 'r') as file:
lines = file.readlines()
if len(lines) < 2:
return None, None
try:
point1 = tuple(map(int, lines[-2].strip().split()))
point2 = lines[-1].strip()
if point2 == "0":
return None, None
point2 = tuple(map(int, point2.split()))
except ValueError as e:
print(f"Error reading coordinates from {file_path}: {e}")
return None, None
return point1, point2
def calculate_angle(pt1: Tuple[int, int], pt2: Tuple[int, int], baseline: Tuple[Tuple[int, int], Tuple[int, int]]) -> float:
"""Calculates the angle between the line formed by pt1 and pt2 and the baseline."""
# Check if points are None
if pt1 is None or pt2 is None or baseline[0] is None or baseline[1] is None:
print(f"Invalid points: pt1={pt1}, pt2={pt2}, baseline={baseline}. Skipping.")
return None # Return None to skip the calculation
angle_doctor = find_angle(pt1, pt2, baseline[0], baseline[1], return_absolute_value=False)
if angle_doctor is None:
return None
angle_baseline = find_angle(baseline[0], baseline[1], baseline[0], baseline[1], return_absolute_value=True)
angle_diff = angle_baseline - angle_doctor
return angle_diff
def compare_baselines_to_doctor(masks_dir: str, coords_dir: str, method: int) -> np.ndarray:
"""Compares baseline methods with doctor's coordinates and generates a histogram."""
angle_diffs = []
for filename in os.listdir(masks_dir):
if not filename.endswith('.png'):
continue
mask_path = os.path.join(masks_dir, filename)
coords_path = os.path.join(coords_dir, filename.replace('.png', '.txt'))
doctor_coords = read_doctor_coordinates(coords_path)
if doctor_coords == (None, None):
print(f"Skipping file {coords_path}, invalid doctor coordinates.")
continue
mask = np.array(Image.open(mask_path))
if method == 1:
rightmost, leftmost = find_baseline_carti(mask)
elif method == 2:
rightmost, leftmost = find_baseline(mask)
elif method == 3:
rightmost, leftmost = old_find_baseline(mask)
else:
raise ValueError("Invalid method specified.")
if rightmost is None or leftmost is None:
print(f"Skipping file {mask_path}, invalid baseline points.")
continue
baseline = (leftmost, rightmost)
angle_diff = calculate_angle(doctor_coords[0], doctor_coords[1], baseline)
if angle_diff is not None:
angle_diffs.append(angle_diff)
return np.array(angle_diffs)
def plot_histogram(angle_diffs: np.ndarray):
"""Plots histogram of angle differences."""
plt.figure(figsize=(10, 6))
if angle_diffs.size == 0:
print("No angle differences to plot.")
return
mean_angle = np.mean(angle_diffs)
std_angle = np.std(angle_diffs)
plt.hist(angle_diffs, bins=np.linspace(-15, 15, 30), edgecolor='black')
plt.title(f'Histogram of Angle Differences - Wersja najstarsza\nMean: {mean_angle:.2f}°, Std: {std_angle:.2f}°')
plt.xlabel('Angle Difference (degrees)')
plt.ylabel('Frequency')
plt.xlim(-15, 15)
plt.grid(axis='y')
plt.axvline(0, color='red', linestyle='dashed', linewidth=1)
plt.show()
# Example usage:
masks_directory = './Angles/dane/masks_baseline'
coords_directory = './Angles/dane/templates_baseline'
method_type = 3 # Choose the method (1, 2, or 3)
angle_differences = compare_baselines_to_doctor(masks_directory, coords_directory, method_type)
plot_histogram(angle_differences)
