Wie man den Abstand zwischen realer Position und rekonstruierter Position Radio Tomoghraphie Bildgebung verringert ⇐ Python

[Guest]

Ich hatte den Python -Code im Zusammenhang mit der Bildgebung von Radio Tomoghraphie. Wenn ich diesen Code ausführe, scheint der reale Punkt (gezeigt durch blaue Farbe) und rekonstruierter Punkt (durch rotes Kreuz gezeigt) getrennt zu sein. Ich möchte also den Fehler zwischen dem wahren Standort und der geschätzten Position des Objekts für die Wärmekarte reduzieren.

Code: Select all

import os

import numpy as np

import matplotlib.pyplot as plt

class RTI:
'''this class is used to run RTI algorithm'''

def __init__(self, specific_links, empty_room_rssi_mean, number_of_nodes=12, alpha=100, lamda=0.01, side_pixels=20, area_width=6.5,
area_length=4, thershold=1, sigma=0.05, delta=0.5, kernel=False, reset=False, use_exp=False, dr=0.02):
self.specific_links = specific_links
self.empty_room_rssi_mean = empty_room_rssi_mean
self.init_other_params(number_of_nodes, area_width, area_length, kernel)
self.init_rti_params(alpha, lamda, side_pixels, thershold, sigma, delta)
self.init_others()
self.reset = reset
self.dr = dr
self.use_exp = use_exp

def init_other_params(self, number_of_nodes, area_width, area_length, kernel):
'''other parameters related to RTI but not directly'''
self.number_of_nodes = number_of_nodes
self.number_of_links = len(self.specific_links)  # Dynamic based on specific_links
self.area_width = area_width
self.area_length = area_length
self.is_kernel = kernel
self.position = np.zeros(shape=(2, 1))
self.pixels_positions = []
self.links_distnaces = {}

def init_rti_params(self, alpha, lamda, side_pixels, thershold, sigma, delta):
'''RTI's params'''
self.alpha = alpha
self.lamda = lamda
self.side_pixels = side_pixels
self.thershold = thershold
self.sigma = sigma
self.delta = delta

def init_others(self):
"""init variables needed by RTI"""
self.pixel_width = self.area_width / self.side_pixels
self.pixel_length = self.area_length / self.side_pixels
self.image = np.zeros(shape=(self.side_pixels, self.side_pixels))
self.weights = []
self.prev_pos = np.array([0, 0])

def create_pixel_positions_matrix(self):
"""create matrix that contains pixels' positions"""
x = self.pixel_length / 2
y = self.pixel_width / 2
for _ in range(0, self.side_pixels):
for _ in range(0, self.side_pixels):
self.pixels_positions.append([x, y])
x = x + self.pixel_length
x = self.pixel_length / 2
y = y + self.pixel_width

def create_link_distance_matrix(self):
"""create link distance matrix based on specific_links"""
for node_a, node_b in self.specific_links:
pos_a = np.array(self.nodes_positions[node_a])
pos_b = np.array(self.nodes_positions[node_b])
self.links_distnaces[f"{node_a}{node_b}"] = {
"node_a": pos_a,
"node_b": pos_b,
"distance": np.linalg.norm(pos_a - pos_b)
}

def use_ellipical_model(self, d):
return 1 / np.sqrt(d)

def use_exponential_model(self, l):
return np.exp((-1 / 2) * (l / self.dr))

def create_weights_matrix(self):
columns = []
for _, link in self.links_distnaces.items():
rows = []
d = link["distance"]
for pixel_pos in self.pixels_positions:
pixel_pos = np.array(pixel_pos)
d1 = np.linalg.norm(link["node_a"] - pixel_pos)
d2 = np.linalg.norm(link["node_b"] - pixel_pos)
if self.use_exp:
l = d1 + d2 - d
rows.append(self.use_exponential_model(l) if l > 0 else 0)
else:
rows.append(self.use_ellipical_model(d) if (d1 + d2 < d + self.lamda) else 0)
columns.append(rows)
self.weights = np.array(columns)

def create_PI_matrix(self):
if self.use_exp:
cov = np.linalg.inv(self.load_or_create_img_cov_matrix())
self.PI = np.linalg.inv(self.weights.T @ self.weights + self.alpha * cov) @ self.weights.T
else:
self.PI = np.linalg.inv(self.weights.T @ self.weights + self.alpha * np.identity(self.side_pixels ** 2)) @ self.weights.T

def load_or_create_img_cov_matrix(self) ->  np.array:
cov_file = f"params\\img_cov_{self.number_of_nodes}_{self.side_pixels}_{self.delta}.dat"
if os.path.exists(cov_file) and not self.reset:
self.cov = np.fromfile(cov_file).reshape(self.side_pixels ** 2, self.side_pixels ** 2)
return self.cov
else:
self.cov = np.zeros(shape=(self.side_pixels ** 2, self.side_pixels ** 2))
for i in range(0, self.side_pixels ** 2):
for j in range(0, self.side_pixels ** 2):
self.cov[i, j] = np.exp(-1 * np.linalg.norm(np.array(self.pixels_positions[i]) - np.array(self.pixels_positions[j])) / self.delta)
self.cov.tofile(cov_file)
return self.cov

def run_rti(self, z):
z = np.array(z)
y = self.empty_room_rssi_mean - z
y[y < self.thershold] = 0
x = self.PI @ y
self.position = self.pixels_positions[x.argmax()]
self.image = x.reshape(self.side_pixels, self.side_pixels)
self.image = (self.image - np.min(self.image)) / np.ptp(self.image)
self.image = np.flip(self.image, axis=0)

def start_rti(self):
self.create_pixel_positions_matrix()
self.create_link_distance_matrix()
self.create_weights_matrix()
self.create_PI_matrix()

def on_close(self, event):
plt.ioff()
plt.close()
exit(0)

def config_drawing(self):
plt.ion()
self.fig, self.ax = plt.subplots()
self.fig.canvas.mpl_connect('close_event', self.on_close)
plt.xlabel('x in meters', fontsize=18)
plt.ylabel('y in meters', fontsize=16)
self.fig.suptitle("RTI Image")
self.aimx = self.ax.imshow(np.random.rand(self.side_pixels, self.side_pixels), interpolation='quadric', extent=[0, self.area_length, 0, self.area_width], vmin=0, vmax=5e-12, cmap="jet")
self.fig.colorbar(mappable=self.aimx)

def draw(self):
self.aimx.set_clim(vmin=0, vmax=self.image.max())
self.aimx.set_data(self.image)
self.fig.canvas.flush_events()

def save_fig(self, path):
'''save rti fig'''
self.fig.savefig(path)

def create_links_columns(self) ->  list:
'''create links names'''
columns = []
for i in range(1, self.number_of_nodes):
for j in range(i + 1, self.number_of_nodes + 1):
columns.append(f"{i}{j}")
return columns

@property
def get_position(self):
'''rti position'''
return self.position

# Define node positions (same for both views)
nodes_positions = {1:[3.52, 6.05],2:[3.65, 0.34],3:[0.46, 0.50],4:[0.48, 6.05],5:[2.03, 5.90],6:[3.50, 4.64],7:[3.48, 3.12],8:[3.54, 1.70],9:[1.88, 0.51],10:[0.46, 1.76],11:[0.50, 3.21],12:[0.56, 4.59]}

top_links = [(1,2),(1, 3), (1, 4), (2, 3), (2, 4), (3,4), (5, 6), (5, 7),(5, 8),(5,9),(5,10),(5,11),(5,12), (6, 7),(6,8),(6,9),(6,10),(6,11),(6,12), (7,8),(7,9),(7,10),(7,11),(7,12),(8,9),(8,10),(8,11),(8,12),(9,10),(9,11),(9,12),(10,11),(10,12),(11,12)]

top_empty_rssi = np.array([-22.41, -16.43, -40.45, -22.22, -23.41, -15.89, -16.22, -10.86, -21.21, -17.6, -17.3, -8.11, -14.94, -9.38, -16.64, -23.38, -23.94, -13.63, -14.81, -11.32, -17.58, -21.98, -27.37, -12.53, -19.03, -19.23, -15.52, -14.88, -10.24, -21.97, -32.28, -10.98, -17.36, -21.71])

side_links =[(1,5),(1, 6), (1, 7), (1, 8),(1,9),(1,10),(1,11),(1,12), (2, 5), (2,6), (2, 7), (2, 8),(2,9),(2,10),(2,11),(2,12),(3, 5), (3, 6),(3,7), (3,8), (3,9),(3,10),(3,11),(3,12),(4, 5), (4, 6),(4,7),(4,8),(4,9),(4,10),(4,11),(4,12)]
side_empty_rssi = np.array([-15.79, -12.07, -24.71, -20.56, -26.76, -28.42, -17.32, -19.68, -29.4, -30.24, -24.94, -24.71, -21.7, -27.22, -28.63, -26.77, -13.34, -13.17, -22.54, -16.2, -12.02, -13.23, -15.33, -22.65, -24.75, -9.59, -15.18, -23.67, -19.89, -13.66, -12.95, -14.29])

# Create RTI instances with correct configurations
rti_top = RTI(
specific_links=top_links,
empty_room_rssi_mean=top_empty_rssi,
number_of_nodes=12,
alpha=1000,
lamda=0.02
)
rti_side = RTI(
specific_links=side_links,
empty_room_rssi_mean=side_empty_rssi,
number_of_nodes=12,
alpha=1000,
lamda=0.02
)

# Set node positions and initialize
rti_top.nodes_positions = nodes_positions
rti_side.nodes_positions = nodes_positions
rti_top.start_rti()
rti_side.start_rti()

# Run RTI with measurement data
#rti_top.run_rti(np.array([-22.4, -15.94, -19.2, -22.41, -28.42, -17.63, -9.98, -10.86, -21.44, -17.9, -14.14, -11.52, -24.61, -10.21, -18.39, -23.82, -23.46, -15.33, -12.16, -12.47, -18.77, -21.99, -17.42, -12.59, -19.21, -20.08, -15.47, -17.92, -9.75, -17.34, -44.52, -10.44, -20.56, -15.45]))
#rti_side.run_rti(np.array([-16.7, -9.94, -29.74, -18.89, -27.48, -22.97, -13.25, -27.92, -33.98, -29.96, -23.69, -23.97, -21.24, -28.26, -35.19, -28.72, -15.63, -13.29, -21.84, -16.66, -11.96, -11.66, -14.98, -21.76, -26.54, -17.59, -19.56, -19.96, -21.78, -15.22, -12.84, -11.08]))
#rti_top.run_rti(np.array([-23.45, -13.79, -21.76, -24.18, -24.55, -17.56, -11.7, -15.35, -21.01, -20.33, -20.47, -8.91, -14.35, -8.42, -15.64, -25.42, -23.33, -10.84, -20.78, -12.77, -25.06, -25.73, -18.97, -14.1, -19.24, -19.16, -14.54, -16.21, -9.72, -21.69, -25.93, -11.82, -15.82, -18.49]))  # Dataset 2 (Top View)
#rti_side.run_rti(np.array([-20.48, -11.84, -21.65, -16.08, -26.59, -18.26, -28.9, -12.97, -29.1, -41.37, -24.72, -24.67, -22.42, -25.72, -36.02, -26.87, -15.9, -13.4, -23.08, -17.37, -12.96, -14.19, -17.39, -19.36, -24.94, -15.29, -17.7, -20.05, -19.85, -14.68, -14.04, -10.84]))  # Dataset 2 (Side View)

# Plot results before optimization
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.imshow(rti_top.image, interpolation='quadric', extent=[0, rti_top.area_length, 0, rti_top.area_width],
vmin=0, vmax=1, cmap="jet")
plt.colorbar()
plt.title("Top View Before Optimization")
plt.xlabel('x (m)')
plt.ylabel('y (m)')

plt.subplot(1, 2, 2)
plt.imshow(rti_side.image, interpolation='quadric', extent=[0, rti_side.area_length, 0, rti_side.area_width],
vmin=0, vmax=1, cmap="jet")
plt.colorbar()
plt.title("Side View Before Optimization")
plt.xlabel('x (m)')
plt.ylabel('y (m)')

plt.tight_layout()
plt.show()
print("Top Position Before Optimization:", rti_top.get_position)
print("Side Position Before Optimization:", rti_side.get_position)

class PSO:
def __init__(self, rti_top, rti_side, num_particles=80, max_iter=50, w=0.9, c1=2,  c2=2):
self.rti_top = rti_top
self.rti_side = rti_side
self.num_particles = num_particles
self.max_iter = max_iter
self.w = w
self.c1 = c1
self.c2 = c2
self.gbest_position = None
self.gbest_value = float('inf')
self.particles = []
self.initialize_particles()
self.rmse_history = []  # Add this line to track RMSE history

def initialize_particles(self):
for _ in range(self.num_particles):
alpha = np.random.uniform(10, 1000)
lamda = np.random.uniform(0.01, 1)
side_pixels = np.random.randint(10, 50)  # Ensure side_pixels is within a valid range
particle = {
'position': np.array([alpha, lamda, side_pixels]),
'velocity': np.random.uniform(-1, 1, 3),
'pbest_position': np.array([alpha, lamda, side_pixels]),
'pbest_value': float('inf')
}
self.particles.append(particle)

def evaluate(self, particle):
alpha, lamda, side_pixels = particle['position']
side_pixels = int(max(1, side_pixels))  # Ensure side_pixels is an integer and positive

# Create new RTI instances for each particle
rti_top = RTI(
specific_links=top_links,
empty_room_rssi_mean=top_empty_rssi,
number_of_nodes=12,
alpha=alpha,
lamda=lamda,
side_pixels=side_pixels
)
rti_side = RTI(
specific_links=side_links,
empty_room_rssi_mean=side_empty_rssi,
number_of_nodes=12,
alpha=alpha,
lamda=lamda,
side_pixels=side_pixels
)
rti_top.nodes_positions = nodes_positions
rti_side.nodes_positions = nodes_positions
rti_top.start_rti()
rti_side.start_rti()

# Calculate RMSE for all datasets
rmse_total = 0
for rss_data_top, rss_data_side, ground_truth_position in zip(rss_datasets_top, rss_datasets_side, ground_truth_positions):
rti_top.run_rti(rss_data_top)
rti_side.run_rti(rss_data_side)
estimated_position_top = rti_top.get_position
estimated_position_side = rti_side.get_position
rmse_top = np.sqrt(np.mean((estimated_position_top - ground_truth_position)**2))
rmse_side = np.sqrt(np.mean((estimated_position_side - ground_truth_position)**2))
rmse_total += (rmse_top + rmse_side) / 2  # Average RMSE for top and side views

return rmse_total / len(rss_datasets_top)  # Return average RMSE

def optimize(self):
self.rmse_history = []  # Reset RMSE history for new optimization
for iteration in range(self.max_iter):
for particle in self.particles:
fitness_candidate = self.evaluate(particle)
if fitness_candidate < particle['pbest_value']:
particle['pbest_value'] = fitness_candidate
particle['pbest_position'] = particle['position']
if fitness_candidate < self.gbest_value:
self.gbest_value = fitness_candidate
self.gbest_position = particle['position']
self.rmse_history.append(self.gbest_value)  # Record best RMSE
for particle in self.particles:
inertia = self.w * particle['velocity']
cognitive = self.c1 * np.random.rand() * (particle['pbest_position'] - particle['position'])
social = self.c2 * np.random.rand() * (self.gbest_position - particle['position'])
particle['velocity'] = inertia + cognitive + social
particle['position'] = particle['position'] + particle['velocity']
print(f"Iteration {iteration + 1}, Best RMSE: {self.gbest_value}")
if self.gbest_value <  0.1:
break
return self.gbest_position, self.gbest_value, self.rmse_history

# Define RSS datasets and ground truth positions (example data)
rss_datasets_top = [
np.array([-22.4, -15.94, -19.2, -22.41, -28.42, -17.63, -9.98, -10.86, -21.44, -17.9, -14.14, -11.52, -24.61, -10.21, -18.39, -23.82, -23.46, -15.33, -12.16, -12.47, -18.77, -21.99, -17.42, -12.59, -19.21, -20.08, -15.47, -17.92, -9.75, -17.34, -44.52, -10.44, -20.56, -15.45]),  # Dataset 1 (Top View)
np.array([-23.45, -13.79, -21.76, -24.18, -24.55, -17.56, -11.7, -15.35, -21.01, -20.33, -20.47, -8.91, -14.35, -8.42, -15.64, -25.42, -23.33, -10.84, -20.78, -12.77, -25.06, -25.73, -18.97, -14.1, -19.24, -19.16, -14.54, -16.21, -9.72, -21.69, -25.93, -11.82, -15.82, -18.49]),  # Dataset 2 (Top View)
np.array([-21.31, -14.85, -21.72, -21.08, -23.56, -16.22, -16.45, -12.89, -21.94, -18.61, -21.4, -7.55, -24.87, -8.5, -16.19, -34.29, -20.31, -9.13, -12.96, -10.48, -13.8, -26.57, -24.89, -12.01, -16.93, -19.18, -14.58, -16.03, -9.42, -21.75, -32.24, -11.27, -17.79, -16.38]),
np.array([-22.36, -18.5, -36.37, -21.32, -24.91, -16.5, -18.29, -11.54, -26, -17.67, -17.8, -8.01, -15.56, -5.32, -30.28, -29.48, -26.24, -12.7, -15.48, -17.46, -53.61, -16.18, -20.36, -12.31, -19.11, -18.85, -13.18, -19.96, -10.25, -21.17, -32.6, -10.92, -17.43, -23.66]),

]

rss_datasets_side = [
np.array([-16.7, -9.94, -29.74, -18.89, -27.48, -22.97, -13.25, -27.92, -33.98, -29.96, -23.69, -23.97, -21.24, -28.26, -35.19, -28.72, -15.63, -13.29, -21.84, -16.66, -11.96, -11.66, -14.98, -21.76, -26.54, -17.59, -19.56, -19.96, -21.78, -15.22, -12.84, -11.08]),  # Dataset 1 (Side View)
np.array([-20.48, -11.84, -21.65, -16.08, -26.59, -18.26, -28.9, -12.97, -29.1, -41.37, -24.72, -24.67, -22.42, -25.72, -36.02, -26.87, -15.9, -13.4, -23.08, -17.37, -12.96, -14.19, -17.39, -19.36, -24.94, -15.29, -17.7, -20.05, -19.85, -14.68, -14.04, -10.84]),  # Dataset 2 (Side View)
np.array([-18.59, -13.54, -22.36, -16.14, -33.09, -22.31, -31.15, -20.93, -30.52, -27.54, -24.46, -25.06, -20.54, -29.35, -31.05, -26.44, -12.96, -13.42, -23.58, -15.5, -12.44, -12.99, -15.13, -23.25, -18.7, -9.58, -16.09, -23.46, -19.93, -13.35, -12.96, -13.7]),
np.array([-16.57, -11.54, -27.56, -24.54, -28.98, -29.41, -16.7, -19.89, -32.82, -25.03, -33.11, -32.54, -20.8, -24.16, -35.31, -28.56, -12.97, -12.97, -24.55, -16.74, -11.85, -12.99, -15.57, -24.7, -25.35, -9.56, -18.83, -23.55, -18.86, -14.57, -13.69, -15.26]),

]

ground_truth_positions = [
np.array([0.98, 5.52]),  # Real position for Dataset 1
np.array([2.06, 5.05]),
np.array([2.95, 5.98]),# Real position for Dataset 2
np.array([3.44, 2.59]),

]

for i, (rss_top, rss_side, real_pos) in enumerate(zip(rss_datasets_top, rss_datasets_side, ground_truth_positions)):
# Run RTI with original parameters
rti_top.run_rti(rss_top)
rti_side.run_rti(rss_side)

# Get reconstructed positions
recon_top = rti_top.get_position
recon_side = rti_side.get_position

# Plot
plt.figure(figsize=(12, 6))
plt.suptitle(f"Dataset {i+1} - Before Optimization", fontsize=14)

# Top View
plt.subplot(1, 2, 1)
plt.imshow(rti_top.image, interpolation='quadric', extent=[0, rti_top.area_length, 0, rti_top.area_width],
vmin=0, vmax=1, cmap="jet")
plt.scatter(real_pos[0], real_pos[1], color='blue', marker='o', s=100, label='Real Position')
plt.scatter(recon_top[0], recon_top[1], color='red', marker='x', s=100, label='Reconstructed')
plt.colorbar()
plt.title(f"Top View (Error: {np.linalg.norm(real_pos - recon_top):.2f}m)")
plt.legend()

# Side View
plt.subplot(1, 2, 2)
plt.imshow(rti_side.image, interpolation='quadric', extent=[0, rti_side.area_length, 0, rti_side.area_width],
vmin=0, vmax=1, cmap="jet")
plt.scatter(real_pos[0], real_pos[1], color='blue', marker='o', s=100, label='Real Position')
plt.scatter(recon_side[0], recon_side[1], color='red', marker='x', s=100, label='Reconstructed')
plt.colorbar()
plt.title(f"Side View (Error:  {np.linalg.norm(real_pos - recon_side):.2f}m)")
plt.legend()

plt.tight_layout()
plt.show()

# =============================================
# Run PSO Optimization with All Datasets
# =============================================
pso = PSO(rti_top, rti_side)
best_position, best_value, rmse_history = pso.optimize()

# Create optimized RTI instances
alpha_opt, lamda_opt, side_pixels_opt = best_position
rti_top_opt = RTI(
specific_links=top_links,
empty_room_rssi_mean=top_empty_rssi,
number_of_nodes=12,
alpha=alpha_opt,
lamda=lamda_opt,
side_pixels=int(side_pixels_opt)
)
rti_side_opt = RTI(
specific_links=side_links,
empty_room_rssi_mean=side_empty_rssi,
number_of_nodes=12,
alpha=alpha_opt,
lamda=lamda_opt,
side_pixels=int(side_pixels_opt)
)
rti_top_opt.nodes_positions = nodes_positions
rti_side_opt.nodes_positions = nodes_positions
rti_top_opt.start_rti()
rti_side_opt.start_rti()

# =============================================
# Display Results After Optimization
# =============================================
for i, (rss_top, rss_side, real_pos) in enumerate(zip(rss_datasets_top, rss_datasets_side, ground_truth_positions)):
# Run RTI with optimized parameters
rti_top_opt.run_rti(rss_top)
rti_side_opt.run_rti(rss_side)

# Get reconstructed positions
recon_top_opt = rti_top_opt.get_position
recon_side_opt = rti_side_opt.get_position

# Plot
plt.figure(figsize=(12, 6))
plt.suptitle(f"Dataset {i+1} - After Optimization", fontsize=14)

# Top View
plt.subplot(1, 2, 1)
plt.imshow(rti_top_opt.image, interpolation='quadric', extent=[0, rti_top_opt.area_length, 0, rti_top_opt.area_width],
vmin=0, vmax=1, cmap="jet")
plt.scatter(real_pos[0], real_pos[1], color='blue', marker='o', s=100, label='Real Position')
plt.scatter(recon_top_opt[0], recon_top_opt[1], color='red', marker='x', s=100, label='Reconstructed')
plt.colorbar()
plt.title(f"Top View (Error: {np.linalg.norm(real_pos - recon_top_opt):.2f}m)")
plt.legend()

# Side View
plt.subplot(1, 2, 2)
plt.imshow(rti_side_opt.image, interpolation='quadric', extent=[0, rti_side_opt.area_length, 0, rti_side_opt.area_width],
vmin=0, vmax=1, cmap="jet")
plt.scatter(real_pos[0], real_pos[1], color='blue', marker='o', s=100, label='Real Position')
plt.scatter(recon_side_opt[0], recon_side_opt[1], color='red', marker='x', s=100, label='Reconstructed')
plt.colorbar()
plt.title(f"Side View (Error: {np.linalg.norm(real_pos - recon_side_opt):.2f}m)")
plt.legend()

plt.tight_layout()
plt.show()

# =============================================
# Print Final Results
# =============================================
print("\nOptimized Parameters:")
print(f"Alpha: {alpha_opt:.2f}, Lambda: {lamda_opt:.2f}, Side Pixels: {int(side_pixels_opt)}")
print(f"Final Average RMSE: {best_value:.4f} meters")

# Plot RMSE convergence
plt.figure()
plt.plot(rmse_history)
plt.title("PSO Convergence")
plt.xlabel("Iteration")
plt.ylabel("RMSE (meters)")
plt.grid(True)
plt.show()