Как уменьшить расстояние между реальным положением и реконструированным положением радиомогарафий ⇐ Python

[Anonymous]

у меня был код Python, связанный с радиоможрафийской визуализацией. Когда я запускаю этот код, то реальная точка (показанная синим цветом) и реконструированная точка (показанная Красным Крестом), по -видимому, является отдельной. Таким образом, я хочу уменьшить ошибку между истинным местоположением и предполагаемым местоположением объекта для тепловой карты. < /P>
Любая помощь в этом отношении будет высоко оценена.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) - 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()


Подробнее здесь: https://stackoverflow.com/questions/794 ... ition-radi